Учебное пособие_Численные методы -26 -дек_2014 (97-2003)

}

void trap(){

double integral=0, t; for(t=a+h; t<=b-h; t=t+h){

integral=integral+f(t);

}

integral=h*((f(a)+f(b))/2.0+integral); cout<<"integral trap="<<integral<<endl;

}

void simps(){

double s1=0,s2=0,integral,t; for(t=a+h; t<=b-h; t=t+2.0*h)

s1=s1+f(t);

for(t=a+2.0*h; t<=b-2.0*h; t=t+2.0*h) s2=s2+f(t);

integral=h/3.0*(f(a)+f(b)+4.0*s1+2.0*s2); cout<<"integral simps="<<integral<<endl;

}

void centr(){

double integral=0,tc; for(tc=a+h/2.0; tc<b; tc=tc+h){

integral=integral+f(tc);

}

integral=integral*h;

cout<<"integral centr="<<integral<<endl;

}

void main(){ prav(); lev(); trap(); simps(); centr();

}

Результат выполнения программы представлен на рис. 6.11.

Рис. 6.11. Результат выполнения программы

Пример 6.8

Производство оборудования некоторого вида характеризуется

темпом роста его выпуска K y 1 , где y - прирост выпуска этого

t y

оборудования за промежуток времени t , а y – уровень его производства за единицу времени на момент времени t. Найти общее количество оборудования, произведенного за 10 лет, полагая, что К – известная постоянная величина, единицей времени является год, а начальный момент времени t=0.

191

Решение. Перейдем к пределу при t 0, полагая, что он существует. Будем также полагать, что y является непрерывной функцией от времени

t. Получаем

	y 1			y	'
K lim
t o t		y		y

ln

y '

.

Интегрируя это	равенство в пределах	от 0	до	t, получаем
Kt ln y y e Kt .
Суммарное количество оборудования, например при K=0,05 (5%
ежегодного темпа	роста), выпущенного за	10	лет,	вычисляется

определенным интегралом

Y

10
e	0,05t	dt
0

.

Вычислим данный определенный интеграл с числом разбиений 100,

1000 и 10000.

Решение в Microsoft Excel при n=100 представлено на рис. 6.12.

Рис. 6.12. Решение в Microsoft Excel при n=100

Решение в Microsoft Visual Studio

Код программы С++ взят из предыдущего примера.

Результаты выполнения программ при различных n представлен на рис. 6.13 – рис. 6.15.

192

Рис. 6.13. Результат выполнения программы при n=100

Рис. 6.14. Результат выполнения программы при n=1000

Рис. 6.15. Результат выполнения программы при n=10000

При определении экономической эффективности капитальных вложений встречаются так называемые задачи дисконтирования: определение начальной суммы S через время t по ее конечной величине St при процентной ставке p. При непрерывном начислении процента

конечная сумма вычисляется по формуле	St Se	rt	, где r 0,01p .
Если сумма St также является	функцией времени f(t), то

дисконтированная сумма к моменту времени t составит S f (t)e rt . Полная дисконтированная сумма за время t вычисляется по

		t
формуле Sd			f (t)e	rt	dt.
		0

Пример 6.9
Определим	дисконтированную сумму S d в течении 10
f (t) S0 1 kt , где		S 0 =3000 – начальное капиталовложения,
ежегодная доля их увеличения, r=0,009 процентная ставка.
		10
Решение. Sd	3000(1 0,05t)e 0,009t dt .
		0

лет при k=0,05 –

Вычислим данный определенный интеграл с числом разбиений 10, 100, 1000 и 10000.

193

Решение в Microsoft Visual Studio

Код программы С++ взят из предыдущего примера.

Результаты выполнения программ при различных n представлен на рис. 6.16рис. 6.19.

Рис. 6.16. Результат выполнения программы при n=10

Рис. 6.17. Результат выполнения программы при n=100

Рис. 6.18. Результат выполнения программы при n=1000

Рис. 6.19. Результат выполнения программы при n=10000

194

Контрольные вопросы

1.Какую формулу называют квадратурной?

2.Как определяется остаточный член квадратурной формулы?

3.Какой вид имеет квадратурная формула Ньютона – Котеса?

4.Какой вид имеет интерполяционный многочлен Лагранжа n-й степени?

5.В чём заключается метод прямоугольников?

6.В чём отличие метода левых, правых и центральных прямоугольников?

7.Как вычисляется интеграл по методу центральных прямоугольников?

8.Как вычисляется интеграл по методу трапеций?

9.В чём заключается метод Монте–Карло?

10.Чем отличается метод Монте–Карло с заданной точностью?

195

11.

7. Решение однородных дифференциальных уравнений

Приближённые методы решения дифференциального уравнения – это методы получения аналитических выражений, либо численных значений,

	приближающих		искомое		частное решение
	дифференциального уравнения.
	Для получения приближённого решения в
	виде аналитического			выражения используют
	метод Эйлера (метод ломаных),основанный на
	приближённом вычислении квадратуры. Этот
	метод был предложен в 1768 г. Метод Эйлера
	относится к семейству методов Рунге-Кутты:
	это одношаговый метод численного решения
Леонард Эйлер	задачи	Коши	для	системы обыкновенных
(1707 – 1783)	дифференциальных уравнений.
Основная идея	метода	применительно		к	дифференциальным

уравнениям первого порядка была предложена Карлом Рунге, немецким

	физиком и математиком, в
	1885 г.
	Дальнейшее		развитие
	этот метод получил в работах
	Вильгельма Кутта в 1901 г.
	Французский математик
	Эмиль	Пикар	развил
	ещё один метод после-
	довательных приближений,
Вильгельм Кутта	который был назван в его			Эмиль Пикар
	честь.			(1856 – 1941)
(1867 – 1944)

Общие замечания

Обыкновенными дифференциальными уравнениями называются уравнения с одной независимой переменной. Обыкновенное дифференциальное уравнение любого порядка при помощи замены переменных может быть сведено к системе уравнений первого порядка.

В общем виде преобразование выглядит следующим образом: дифференциальное уравнение n-го порядка

196

заменой переменных v(k)=vk сводятся к системе n уравнений первого порядка

v	v	k 1	,	0 k n 2,
k

где обозначено v0≡v.

В соответствии с изложенным, далее будут рассматриваться системы уравнений первого порядка:

Решение системы n-го порядка зависит от n параметров c1 , c2 ,..., cn .

Единственное решение получается при использовании дополнительных условий для искомой функции. В зависимости от того, каким образом ставятся данные условия, различают три типа задач для обыкновенных дифференциальных уравнений: задача Коши, краевая задача и задача на собственные значения. В задаче Коши все дополнительные условия ставятся в одной точке: vk(x0)=vk,0, 1≤k≤n. Решение отыскивается в некотором интервале x0≤x≤xl.

Если правые части φk	уравнений	непрерывны в некоторой
окрестности начальной точки	x0 , v1,0 , v2,0 ,..., vn,0	, и удовлетворяют условию

Липшица по переменным vk, то решение задачи Коши существует, единственно и непрерывно зависит от координат начальной точки, т.е. задача является корректной. φ

Условие Липшица формулируется следующим образом:

для любых точек x, v1,l , v2,l ,..., vn,l , x, v1,m , v2,m ,..., vn,m .

Методы решения

Можно выделить три типа методов решения обыкновенных дифференциальных уравнений: точные, приближенные и численные. Точные методы предусматривают получение решения в виде комбинации элементарных функций или в виде квадратур от последних. Возможности точных методов ограничены. Приближенные методы сводятся к построению последовательности функций wn(x), имеющих пределом искомую функцию vn(x). Обрывая эту последовательность на каком-то k, получают приближенное решение. Наиболее универсальными методами решения являются численные. Их основной недостаток – возможность получения только частного решения. Следует иметь в виду следующее обстоятельство. Успех от применения численного метода сильно зависит от обусловленности задачи, т.е. задача должна быть хорошо обусловлена,

197

а именно, малые изменения начальных условий должны приводить к малому изменению решения. В противном случае (слабой устойчивости) малые погрешности в начальных данных или погрешности численного метода могут приводить к большим погрешностям в решении.

Далее будут рассматриваться алгоритмы решения задачи Коши на примере одного уравнения первого порядка . Обобщение на случай системы n уравнений осуществляется заменой на и

на , где

v(x) v	v	2	... v	n	,
1

.

7.1. Метод Пикара

Данный метод является представителем класса приближенных методов решения.

Идея метода чрезвычайно проста и сводится к процедуре последовательных приближений для решения интегрального уравнения, к которому приводится исходное дифференциальное уравнение.

Пусть поставлена задача Коши

,

v(x	0	)

Проинтегрируем выписанное уравнение

v	0

.

(7.1)

. (7.2)

Процедура последовательных приближений метода Пикара реализуется согласно следующей схеме:

,

(7.3)

причем .

Докажем сходимость метода Пикара. Пусть в некоторой ограниченной области правая часть непрерывна и, кроме того, удовлетворяет условию Липшица по переменной v, т.е.

, где L - некоторая константа.

Всилу ограниченности области имеют место неравенства

xx0 l1 , v v0 V .

Вычтем из (7.3) формулу (7.2), получим выражения для модулей правой и левой частей:

,

или

198

.

Окончательно, используя условие непрерывности Липшица, получим

		x
zs (x)	L		zs 1	(t)dt ,
		x
		0

где zs (x) ys (x) v(x) − погрешность приближенного решения. Последовательное применение формулы (7.4) при

следующую цепочку соотношений при учете того, что

s 1,2,...

(7.4)

дает

z	0	(x) v		0	v(x)			V ,
z	1	(x) LV x x					0	,
z		(x)	1		2	V (x			x		)
	2		2	L						0
...
z		(x)	1		s	V (x			x		)
	s			L						0
			s!

2

s

,

.

Так как

Заменяя

погрешности

x x0 l1 , то имеем

по формуле Стирлинга, приближенного решения

.

окончательно получим оценку

		V	el L			s
z	(x)
				1
s		2 s		s		.
Из (7.4) следует, что при
		s		модуль погрешности

приближенное решение равномерно сходится к точному.

z	s	(x)

(7.5)

, т.е.

7.2. Методы Рунге-Кутта

Данные методы являются численными.

На практике применяются методы Рунге-Кутта, обеспечивающие построение разностных схем (методов) различного порядка точности. Наиболее употребительны схемы (методы) второго и четвертого порядков. Их мы и рассмотрим далее.

Предварительно введем некоторые понятия и определения. Сеткой на отрезке a, b называется фиксированное множество точек этого отрезка

. Функция, определенная в данных точках, называется сеточной функцией. Координаты точек xi удовлетворяют условиям

	a x0 x1 x2 ... xN 2	xN 1 xN	b
			.
Точки	являются узлами сетки.	Равномерной сеткой на a, b
называется множество точек			, где
шаг сетки.
	199

При решении дифференциальных уравнений приближенным методом основным является вопрос о сходимости. Применительно к разностным

методам традиционно более	употребительно понятие сходимости при
h 0. Обозначим значения	сеточной функции	yi	значения точного

решения дифференциального уравнения (7.1)		в узле i v(xi ) ( yi являются
приближенными значениями	v(xi ) ). Сходимость при h→0 означает
следующее. Фиксируем точку x и строим совокупность сеток			таким
образом, что h→0 и	(при этом	). Тогда считают, что
численный метод сходится в точке x, если		yi v(xi ) 0 при	h 0 xi x .

Метод сходится на отрезке a, b , если он сходится в каждой точке

x a,b .

Говорят, что метод имеет p-й порядок точности, если можно найти такое
число p 0 , что	yi v(xi ) O h	p	при h 0.

Введем далее понятие невязки, или погрешности аппроксимации разностного уравнения, заменяющего заданное дифференциальное уравнение, на решении исходного уравнения, т.е. невязка представляет

собой результат подстановки точного решения уравнения (7.1)

v(x)

в

разностное уравнение. Например, (7.1) можно заменить простейшим разностным уравнением

		.
Тогда невязка определится выражением		.
Приближенное решение не совпадает, вообще говоря, c Ui, поэтому
невязка	в i-й точке не равна нулю.	Вводят следующее определение:

численный метод аппроксимирует исходное дифференциальное

уравнение, если	при h 0, и имеет p-й порядок точности, если
.
Доказывается,	что порядок точности численного метода решения

дифференциального уравнения совпадает с порядком аппроксимации при достаточно общих предположениях.

Теперь перейдем к анализу схем Рунге-Кутта. Сначала обратимся к схемам второго порядка точности.

Используя формулу Тейлора, решение дифференциального уравнения (7.1) можно представить в виде

v

v

h v

1

h2 v ...

n 1

где обозначено vn

v(xn ),

n

n

n

2 n

n

,

(7.6)

vn v (xn ),

hn

xn 1

xn .

Отметим,

что

согласно

уравнению

(7.1)

,

.

Далее удерживаем только выписанные члены ряда. Представим вторую производную следующим образом:

200