kontrolnye / Обчислювальна математика / Лекции / r9

Приближенное решение не совпадает, вообще говоря с ui, поэтому невязка ψi в i-ой точке не равна нулю. Вводят следующее определение: численный метод аппроксимирует исходное дифференциальное уравнение, если ψi → 0 при h → 0 и имеет p-й порядок точности, если ψi = O(hp).

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

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

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

	vn+1 = vn + hnvn0	+	1	hn2 vn00	+ . . . ,	(9.6)
			2
	где обозначено vn = v(xn), vn0 = v0(xn), hn = xn+1	− xn .Отметим, что согласно (9.1) vn0				= ϕ(xn, vn),
	vn00 = ϕx0 (xn, vn)ϕv0 (xn, vn).

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

v00 = (v0 )0 = ϕ(˜x, v˜) − ϕ(xn, vn),

n n x

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

x˜ = xn + γh, v˜ = vn + δh.

Обозначим приближенное значение решения в узле с номером n через yn (именно это решение будет получаться после того, как мы ограничим ряд членами с порядком не выше второго).

Имеем

•Назад •Первая •Предыдущая •Следующая •Последняя •Перейти •Предметный указатель

yn+1 = yn + hnϕ(xn, yn) +	1	h2	[	ϕ(xn + γhn, yn + δhn) − ϕ(xn, yn)	] =
				x
	2 n

= yn + hn[βϕ(xn, yn) + αϕ(xn + γhn, yn + δhn)].

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

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

yn+1 = yn + hn { β · ϕ(xn, yn) + α · [ϕ(xn, yn)+

+ϕ0x(xn, yn)γhn + ϕ0y(xn, yn)δhn] } =

= yn + (α + β)hnϕ(xn, yn) + αh2

[γϕ0 (xn, yn) + δϕ0 (xn, yn)].

(9.7)

n

x

y

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

1

1

α + β = 1, αγ =

, αδ =

ϕ(xn, yn).

2

2

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

1

1

β = 1 − α, γ =

, δ =

ϕ(xn, yn)

2α

2α

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

1

hn

(9.8)

yn+1 = yn + hn {(1 − α) ϕ(xn, yn) + αϕ(xn +

hn, yn +

ϕ0(xn, yn)] } .

2α

2α

•Назад •Первая •Предыдущая •Следующая •Последняя •Перейти •Предметный указатель

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

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

Впрактике расчетов используют формулы (9.8) при значениях параметра α = 12, α = 1.

Случай α =	1	.
	2
Из (9.8) выводим
		yn+1	= yn +	h	[ϕ (xn, yn) + ϕ (xn + hn, yn + hnϕ (xn, yn))] .	(9.9)
			2

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

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

y¯n+1 = yn + hnϕ (xn, yn) .

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

y¯n0 +1 = ϕ (xn+1, y¯n+1) .

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

yn0 + 12 = 12 ϕ (xn, yn) + y¯n0 +1 ,

•Назад •Первая •Предыдущая •Следующая •Последняя •Перейти •Предметный указатель

4. Рассчитывается значение функции в (n + 1)-м узле

yn+1 = yn + hyn0 + 12 .

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

Согласно (9.8) получаем

yn+1 = yn + hnϕ[xn + h2n , yn + h2n ϕ(xn, yn)].

Задача решается посредством следующих шагов:

1. Вычисляется значение функции в половинном узле

hn

yn+ 12 = yn + 2 ϕ(xn, yn).

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

y0	1	= ϕ(xn +	hn	, y		1 ).
					n+
n+	2	2				2

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

yn+1 = yn + hnyn0 + 12 .

•Назад •Первая •Предыдущая •Следующая •Последняя •Перейти •Предметный указатель

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

yn+1 = yn + (k1 + 2k2 + 2k3 + k4) /6,	(9.10)
k1 = hnϕ (xn, yn) , k2 = hnϕ (xn + hn/2, yn + k1/2) ,
k3 = hnϕ (xn + hn/2, yn + k2/2) , k4 = hnϕ (xn + hn, yn + k3) .

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

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

Оценку погрешностей легко получить для частного случая правой части дифференциального уравнения

ϕ(x, v) ≡ ϕ(x).

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

yn+1 = yn + h2n [ϕ (xn) + ϕ (xn+1)] ,

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

yn+1 = yn + h2n [ϕ (xn) + 4ϕ (xn + hn/2) + ϕ (xn + hn)] ,

•Назад •Первая •Предыдущая •Следующая •Последняя •Перейти •Предметный указатель

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

Мажорантные оценки погрешности формул трапеций и Симпсона известны (см. раздел 7). Из (7.4) и (7.5) видно, что точность схем Рунге-Кутта достаточно высока.

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

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

можно описать следующей системой уравнений

a0yn + a1yn−1 + . . . + amyn−m

= b

ϕ

+ b

ϕ

n−1

+ . . . + b

ϕ

n−m

,

(9.11)

h

0

n

1

m

где n = m, m + 1, . . ., а ak, bk — числовые коэффициенты, k = 0, 1, 2, . . . , m, a0 6= 0.

Согласно данному уравнению расчет начинается с n = m. В этом случае получается соотношение вида

a0ym + a1ym−1 + . . . + amy0 = b0ϕm + b1ϕm−1 + . . . + bmϕ0, h

т.е. для начала счета надо иметь m начальных значений yi, i = 0, 1, 2, . . . , m − 1. Эти значения yi приходится вычислять каким-либо другим методом, например, методом Рунге-Кутта.

Среди многошаговых методов наиболее распространен метод Адамса, схема реализации которого следует из (9.11) при a0 = −a1 = 1 и ak = 0 для k = 2, 3, . . . , m:

•Назад •Первая •Предыдущая •Следующая •Последняя •Перейти •Предметный указатель

m

yn − yn−1 = X bkϕn−k. h

k=0

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

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

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

m	ak		m
X			X	(9.12)
k=0	h	yn−k =	bkϕ(xn−k, yn−k), n = m, m + 1, . . . .
			k=0

Однородное разностное уравнение, соответствующее (9.12), имеет вид

m

X

akyn−k = 0.

(9.13)

k=0

Говорят, что уравнение (9.13) устойчиво по начальным данным, если существует постоянная M, не зависящая от n, такая, что при любых начальных данных y0, y1, . . . , ym−1 имеет место неравенство

|yn| 6 M · max |yj|, n = m, m + 1, . . . .

06j6m−1

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

•Назад •Первая •Предыдущая •Следующая •Последняя •Перейти •Предметный указатель

yn−k = qn−k . Подставляя данное yn−k в (9.13) и сокращая на qn−m, получим характеристическое уравнение для нахождения q:

a0qm + a1qm−1 + . . . + am−1q + am = 0.

(9.14)

Справедлива следующая теорема. Для устойчивости уравнения (9.13) по начальным данным необходимо и достаточно, чтобы выполнялось условие на корни, а именно: все корни q1, q2, . . . , qm характеристического уравнения должны располагаться внутри или на границе единичного круга комплексной плоскости, причем на границе не должно быть кратных корней.

Доказывается следующее утверждение.

Теорема 1 Пусть 0 6 nh 6 T , условие на корни выполнено, |yi − v(xi)| → 0 при h → 0, i = 0, 1, 2, . . . , m−1, и разностное уравнение (9.12) аппроксимирует исходное дифференциальное уравнение . Тогда решение

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

[0, T ].

Сформулированное условие устойчивости, базирующееся на анализе расположения корней характеристического уравнения (9.14), является весьма общим. Конкретизируем вопрос об устойчивости разностного уравнения применительно к асимптотически устойчивым решениям уравнения (9.1). Пусть ϕ(x, v(x)) = λ · v(x), λ < 0, т.е.

dv	= λ · v(x).	(9.15)
dx

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

|v(x + h)| 6 |v(x)| .

(9.16)

•Назад •Первая •Предыдущая •Следующая •Последняя •Перейти •Предметный указатель

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

	yn+1 − yn	= λ	·	y	n	, n = 0, 1, 2, . . . ,	(9.17)
	h

или yn+1 = (1 + hλ)yn , т.е. q = 1 + hλ .

Оценка (9.16) будет выполнена для (9.17) только в том случае, если |q| 6 1 , так как тогда |yn+1| 6 |yn|

. Из |q| 6 1 следует ограничение на шаг h: 0 6 h 6 |λ2 | .

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

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

Эйлера, имеющий первый порядок аппроксимации
yn+1 − yn		= λ			·	y	n+1	.	(9.18)
	h
Из (9.18) следует
yn+1 =	yn			=			yn		,
	1 − hλ				1 + h |λ |

1

т.е. |q| = 1 + h |λ| < 1 всегда, при любых h > 0.

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

•Назад •Первая •Предыдущая •Следующая •Последняя •Перейти •Предметный указатель

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

m
X	(9.19)
(ak − µ bk)yn−k = 0, n = m, m + 1, . . . ,
k=0
где µ = hλ — в общем случае комплексный параметр.
Характеристическое уравнение для (9.19) имеет вид
m
X	(9.20)
(ak − µbk)qm−k = 0.

k=0

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

Областью устойчивости метода (9.12) называется множество точек комплексной плоскости µ = hλ, для которых данный метод, примененный к уравнению специального вида (9.15), является устойчивым.

Для явного метода Эйлера условие устойчивости

+ µ

| 26

1 при комплексном µ = µ

+ iµ

(µ

=

|1 2

0

1

0

Reµ, µ1 = Im µ ) выглядит следующим образом: (µ0 + 1)

+ µ1 6 1, т.е. областью устойчивости является

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

1

6 1

соответствует неравенству (1 − µ0)2 + µ12 > 1 , т.е.

Для неявного метода Эйлера условие

|1 − µ|

областью устойчивости является внешняя область круга единичного радиуса с центром в точке (1; 0). Разностный метод называется A-устойчивым, если область его устойчивости включает левую полу-

плоскость Reµ < 0 (или h · Reλ < 0). Следует обратить внимание на то, что уравнение (9.15) асимптотически устойчиво при Reλ < 0. Следовательно A-устойчивый разностный метод является абсолютно устойчивым (т.е. устойчивым при любых h > 0), если устойчиво решение исходного дифференциального уравнения.

Из приведенного выше рассмотрения видно, что неявный метод Эйлера обладает свойством A- устойчивости, а явный метод — нет.

•Назад •Первая •Предыдущая •Следующая •Последняя •Перейти •Предметный указатель