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

,

где

~ ~

.

x, v пока неизвестные величины. Пусть

Обозначим приближенное значение решения в узле с номером n

через yn (именно это решение

будет

получаться

после

того, как мы

ограничим ряд членами с порядком не выше второго).

Имеем

y

y

h (x

, y

)

1

h

2

(x

h

, y

h ) (x

, y

)

n

n

n

n

n

n

n 1

n

n

n

n

2

n

x

y

n

h

(x

n

, y

) (x

n

h

, y

n

h ) .

n

n

n

n

Введенные здесь параметры α, β, γ и δ подлежат определению.

Раскладывая правую часть в ряд Тейлора и приводя подобные члены,

получим последовательно

yn 1 yn hn

(xn , yn ) (xn , yn ) x (xn , yn ) hn y (xn , yn ) hn

y ( )h (x , y ) h2 (x , y ) (x , y ) .

(7.7)

n

n

n

n

n

x

n n

y

n

n

Условием выбора параметров α, β, γ и δ поставим близость выражения (7.7) ряду (7.6), тогда

.

Один параметр остается свободным. Пусть это будет α, тогда

и окончательно из (7.7) с учетом найденных отношений для β, γ и δ получим

.

(7.8)

Соотношение (7.8) описывает однопараметрическое семейство двучленных формул Рунге-Кутта.

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

значениях параметра , α=1.

Рассмотрим случай

y	n 1

. Из (7.8) выводим

	yn	h	(xn	, yn ) xn	hn , yn(xn , yn ) .
		2

(7.9)

Применение формулы (7.9) сводится к следующей последовательности шагов:

1. Вычисляется грубо значение функции yn 1 (по схеме ломаных)

.

201

2. Определяется наклон интегральной кривой в точке (xn 1 , yn 1 )

.

3. Находится среднее значение производной функции на шаге

h	n

	y			1	(x		, y		)
		1		2		n		n
	n
		2

4. Рассчитывается значение функции в

y
n 1
	.

(n+1)-м узле

y	n 1

y	n

	hy	1
	n
		2

.

Данная схема имеет специальное название «предиктор-корректор». Рассмотрим случай α=1. Согласно (7.8) получаем

y

y

h

x

h

, y

h

(x

, y

)

n

n

n 1

n

n

n

2

n

2

n

n

Задача решается посредством следующих шагов: 1. Вычисляется значение функции в половинном

	.

узле

y			y			h	(x
						n
	n	1		n		2	n
		2

2.Определяется значение производной в

, y	n	)
		.

узле n+(1/2)

					h
y			x		n	, y
n	1		n		2		n	1
	2							2		.

3.Находится значение функции в(n+1) - м узле

y	n 1	x	n	h	n	y	1
						n
							2	.

Помимо рассмотренных выше двучленных схем, широкое распространение в практике расчетов имеют схемы Рунге-Кутта

четвертого порядка точности. Приведем соответствующие формулы:

y

n 1

y

k

2k

2k

k

/ 6

,

n

1

2

3

4

k

h x

, y

,

k

h x

h

/ 2, y

k

/ 2

,

(7.10)

1

n

n

n

1

n

n

n

n

1

k

h x

h

/ 2, y

k

/ 2

2

n

n

n

n

1

,

k h x h / 2, y k / 2

,

k

h x

h

, y

k

.

3

n

n

n

n

2

4

n

n

n

n

3

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

Погрешности приведенных схем Рунге-Кутта определяются максимальными значениями соответствующих производных.

Оценку погрешностей легко получить для частного случая правой части дифференциального уравнения x, v x . В этом случае решение уравнения может быть сведено к квадратуре, и все схемы разностного решения переходят в формулы численного интегрирования. Например, схема (7.9) принимает вид

202

y	n 1

y	n

hn

2

(x	)
n

(x	)
n 1

,

то есть имеет вид формулы трапеций, а схема

y		y			h	(x	) 4 (x	h	/
					n
	n 1		n		2	n	n	n

(7.10) переходит в схему

2) (x	h )	,
n	n

представляющую собой формулу Симпсона с шагом	n	.
	h / 2

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

следующими соображениями. Если функция

x, v

в правой части

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

Помимо схем Рунге-Кутта, практический интерес представляют многошаговые методы, которые можно описать следующей системой уравнений:

a y		a y	n 1	a	y	n m b			b		b
0	n	1		m				n		n 1
				h		0			1		m
где n m, m 1, ; ak , bk		- числовые коэффициенты;									k 0, m ,

Согласно данному уравнению расчет начинается случае получается соотношение вида

n m

a	0

с

,

0 . n

m

(7.11)

. В этом

a y

a y

m 1

a

y

b

b

b

0

m

1

m

0

m 1

,

h

0

m

1

m

0

т.е. для начала счета надо иметь m начальных значений

yi , i 0, m 1 . Эти

значения

yi

приходится

вычислять

каким-либо

другим

методом,

например, методом Рунге-Кутта.

Среди многошаговых методов наиболее распространен метод

Адамса, схема реализации которого следует из (7.11)

при a0

a1 1 и

ak 0 , k

:

2, m

y

y

m

n

n 1

b

n k

h

k

.

k 0

При b0

0 метод Адамса оказывается явным, а при b0

0 - неявным.

203

7.3. Неявные методы

Введем понятие устойчивости разностного метода. Для этого рассмотрим уже упоминавшееся разностное уравнение многошагового метода

m		m
	ak	yn k bk xn k	, yn k
				(7.12)
k 0 h		k 0	, n m, m 1, .
Однородное разностное уравнение, соответствующее (7.12), имеет вид
		m
		ak yn k	0	(7.13)
		k 0	.

Говорят, что уравнение (7.13) устойчиво по начальным данным, если существует постоянная M, не зависящая от n, такая, что при любых

начальных данных	y0 , y1, , ym 1			имеет место неравенство
	y	n	M			max	y	j	n m, m 1,
						0 j m 1
								,
Вопрос устойчивости				по		начальным данным решается путем

рассмотрения корней так называемого характеристического уравнения,

получаемого из (7.13), если

решение

этого

уравнения искать в виде

yn k q

n k

. Подставляя данное

yn k в (7.13)

и сокращая на

q

n m

получим

характеристическое уравнение для нахождения q:

a q

m

a q

m 1

a

q a

0

.

(7.14)

0

1

m 1

m

Справедлива следующая теорема. Для устойчивости уравнения (7.13) по начальным данным необходимо и достаточно, чтобы выполнялось

условие корней, а именно: все корни

q1, q2 , , qm

характеристического

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

Доказывается следующее				утверждение.			Пусть	0 nh T , условие
корней выполнено,	i	i	0		при	h 0 ,	i 0, m 1, и разностное
	y	x

уравнение (7.12) аппроксимирует исходное дифференциальное уравнение (7.1). Тогда решение разностной задачи (7.12) сходится при h 0 к решению исходной задачи (7.1). Говоря другими словами, из аппроксимации и устойчивости по начальным данным следует сходимость на ограниченном отрезке 0,T .

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

уравнения (7.1).

Пусть x, x x , 0 , т.е.

204

d dx

x

.

(7.15)

Решение этого уравнения асимптотически устойчиво, т.е. при любых

x 0	справедлива оценка
	x h x	(7.16)
	.

Потребуем, чтобы и разностное уравнение давало решение, обладающее свойством (7.16). Используя явный метод Эйлера первого порядка аппроксимации, получим разностный аналог (7.15)

y	n 1	y	n
	h

y	n

,

n

0,1,2,

.

(7.17)

или

Оценка (7.16) будет q 1 , так как тогда y

y	n 1	1 h y	n , т.е.	q 1 h	.

выполнена для (7.17) только в том случае, если n 1 yn . Из q 1 следует ограничение на шаг

h:

0 h

2

.

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

В качестве примера абсолютно устойчивого метода традиционно рассматривается неявный метод Эйлера, имеющий первый порядок аппроксимации

Из (7.18) следует

h

0

. т.е.

q

		y	n 1
			h
		1	1
		h
	1

y	n	y
			n 1

всегда,

.

при любых

h

0

.

(7.18)

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

Запишем разностное уравнение (7.12) для задачи (7.15):

m

ak bk yn k 0 , k 0

n

m, m

1,

,

(7.19)

где

h - в общем случае комплексный параметр. Характеристическое уравнение для (7.19) имеет вид

m
ak	bk q	m k	0 .
k 0

При малых µ корни (7.20) близки к корням (7.13).

(7.20)

205

Областью устойчивости метода (7.12) называется множество точек

комплексной плоскости

h , для которых данный метод, примененный

к уравнению специального вида (7.15), является устойчивым.
Для явного метода Эйлера условие				устойчивости	1 1	при
комплексном		0	i 1 0 Re , 1 Im выглядит следующим образом:
0 1	1 1,	т.е.	областью устойчивости	является круг	единичного
2	2

радиуса, центр которого находится в точке (-1;0) комплексной плоскости.

	Для неявного			метода		Эйлера условие	1	1	соответствует
							1
	неравенству	1 0		1	1, т.е.
						областью устойчивости является внешняя
			2	2

область круга единичного радиуса с центром в точке (1;0).

Разностный метод называется A − устойчивым, если область его

устойчивости включает левую полуплоскость

Re μ 0

(или

h Re μ 0 ).

Следует обратить внимание на то, что уравнение (7.15) асимптотически

устойчиво при	Re 0	. Следовательно A − устойчивый разностный метод

является абсолютно устойчивым (т.е. устойчивым при любых h>0), если устойчиво решение исходного дифференциального уравнения. Из приведенного рассмотрения видно, что неявный метод Эйлера обладает свойством A − устойчивости, а явный метод - нет.

Рассмотрим еще один неявный метод более высокого порядка аппроксимации (второго):

y

n 1

y

n

1

x

, y

x

, y

n 1

n

h

2

n 1

n

.

(7.21)

Этот метод получается заменой интеграла от правой части (7.1) на длине шага по формуле трапеций. Применительно к уравнению (7.15) метод (7.21) выглядит следующим образом:

т.е.

1 0,5μ 1 0,5μ

1

	yn 1		1 0,5μ		yn	,
			1	0,5μ

, если µ ≤ 0, т.е. метод (7.21) относится к A − устойчивым.

Существует доказательство следующих положений:

среди методов (7.12) не существует явных A − устойчивых методов;

среди неявных линейных многошаговых методов нет A − устойчивых методов, имеющих порядок точности выше второго.

A − устойчивые разностные схемы весьма эффективны при решении

так называемых жестких систем уравнений, так как эти методы не накладывают ограничений на шаг h. Рассмотрим подробнее это утверждение.

Система обыкновенных дифференциальных уравнений

206

	d	A
	dx

с независящей от x матрицей A(m×m) называется

	k=	и отношение	s	max Re	k	велико, где	k
				1 k m
				min Re	k
				1 k m

(7.22)

жесткой, если Reλk<0, - собственные числа

матрицы A. Величина зависит от x, то и k жесткости

s называется числом жесткости. Если матрица A -зависят от x, тогда вводится переменное число

	max Re	k
	s
	1 k m
	min Re	k
	1 k m
	и оперируют с величиной sup s x на отрезке интегрирования.

Отличительной особенностью жестких систем является наличие в их решении как быстро, так и медленно убывающих компонент. При x> 0 решение системы практически определяется медленно убывающей компонентой, однако, если воспользоваться явными разностными методами, то быстро убывающая составляющая будет отрицательно влиять на устойчивость, и в результате весь расчет необходимо вести с малым шагом интегрирования. При использовании же неявных методов ограничения на шаг сняты, и его величину определяют из условия достижения нужной точности, не заботясь особо об устойчивости.

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

m	h xn , yn
ak ym k
k 0	,

т.е. рассматривается частный случай метода (7.12), когда b1= b2=… =bm=0, а

b0=1.

При m = 1 и a0 = 1, a1 = −1 имеем

yn yn 1

h xn , yn , т.е. неявный метод

Эйлера. При m = 2 и m = 3 методы выглядят следующим образом:

3

y

2 y

1

y

h x

, y

,

(7.23)

n

n 1

n 2

2

2

n

n

11

y

3y

3

y

1

y

h x

, y

.

(7.24)

n

n 1

n 2

n 3

6

2

3

n

n

Разностное уравнение (7.23) имеет второй порядок точности, а (7.24) - третий. Чтобы найти область устойчивости метода, следует записать аналогичные уравнения для дифференциального уравнения (7.15). Например, (7.23) примет вид

32 yn 2 yn 1 12 yn 2 yn .

207

Соответствующее характеристическое уравнение запишется следующим образом:

3		2	2q	1	0
	q
2				2

.

(7.25)

Наша задача определить область комплексной плоскости µ=µ0+ iµ1, в точках которой оба корня (7.25) по модулю меньше единицы. Оказывается, что эта область целиком располагается в правой плоскости и метод (7.23) является A − устойчивым.

Метод (7.24) относится к так называемым A(α) – устойчивым методам.

7.4. Примеры

Пример 7.1

Дано уравнение y’=xy. Начальные условия x0 =1, y0 =1. Найти первые три точки y1, y2, y3, смещаясь с шагом h = 0,2 методом Эйлера.

Итерационная формула: yi+1=yi+f(xi, yi)(xi+1-xi). Решение. В данном случае yi+1=yi+xiyih.

y1=1+1∙1∙0,2=1,2; y2=1,2+1,2∙1,2∙0,2=1,488; y3=1,488+1,4∙1,488∙0,2=1,905.

Нашли приближенное решение, а теперь найдем точное – аналитически.

x3=1,6 y3*=2,181; δ=y*-y3=0,276.

Пример 7.2

Дано уравнение y’=xy, x0=1, h=0,2. Решить первой,второй и третьей модификацией метода Эйлера.

Решение.

Первая модификация метода Эйлера

1-й шаг:

y 1		1 1 1 0,1 1,1 , (этап 1)
	2
y1		1 1,1 1,1 0,2 1,242, (этап 2)

208

.

2-й шаг:

,

,

(этап 1)

(этап 2)

.

Вторая модификация метода Эйлера

,	(этап 1)
.	(этап 2)
1-й шаг:
,
.
2-й шаг:
	,
	.

Третья модификация метода Эйлера

1-й шаг:

,

.

2-й шаг:

,

.

Пример 7.3

Дано уравнение . Решить методом Рунге – Кутта.

Решение:

,

,

,

209

,

,

,

.

Пример 7.4

	v x x	3	v
Решить методом Пикара уравнение	v 0 0.

Решение этого уравнения не выражается функции.

3	,

через элементарные

			x							x	4
1		(x) 0			3	dt
y				t								,
										4
			0
			x						t		4		3					x	4						1
		(x) 0 t					3																				9
y2														dt						1				2		x			,
										4								4				4			13
			0
									t		4						1					3				x		4			x	13					3					3					1
y	3	(x) 0			t 3								1			2				t 9				dt							3			1		2		x9			2		x18			6		2		x 27
										4					4			13								4				4					4		22			4		13	31		4		13		40
																																13

ит.д.

7.5.Примеры применения дифференциальных уравнений

вэкономике

Дифференциальное исчисление – широко применяемый для экономического анализа математический аппарат [19, 20]. Базовой задачей экономического анализа является изучение связей экономических величин, записываемых в виде функций. Основная задача дифференциального исчисления заключается в следующем: дана функция F(x), требуется найти ее производную, т.е. f (x) F' (x) (например, найти предельные издержки, зная суммарные издержки).

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

210