ВУЗ: Не указан
Категория: Не указан
Дисциплина: Не указана
Добавлен: 02.08.2019
Просмотров: 4704
Скачиваний: 6
Численное решение дифференциальных уравнений
91
отклоняется от решения
y
k
системы разностных уравнений. Для погреш-
ности решения мы будем использовать величину
δ
h
= max
k
|y
k
− u(x
k
)
| =
||y
k
− u(x
k
)
||
∞
,
k = 0, 1, . . . , n.
Говорят, что численное решение имеет
p-ый порядок точности, если су-
ществует число
p > 0 такое, что δ
h
= O(h
p
) при h
→ 0.
Очевидно, в общем случае точное решение дифференциальной задачи
u(x) и решение разностной схемы y
k
,
k = 0, 1, . . . , n не совпадают. Поэто-
му, если решение дифференциального уравнения подставить в разностное
уравнение, последнее перестанет быть равенством из-за появления невяз-
ки. Данная невязка называется погрешностью аппроксимации и обозна-
чается
ψ. Если при уменьшении шага h
→ 0 выполняется ψ = O(h
p
) для
некоторого
p, то говорят, что погрешность аппроксимации имеет порядок
p. Если ψ
k
6→ 0 при h → 0, то говорят, что разностная схема не аппрокси-
мирует дифференциальное уравнение.
Найти аналитически погрешность решения — это часто очень сложная
или неразрешимая задача. Напротив, вычислить погрешность аппрокси-
мации довольно легко. На практике важно знать, какой порядок точности
имеет та или иная разностная схема. Оказывается, что ответить на этот
вопрос можно, зная порядок погрешности аппроксимации.
Теорема 11.1 (без доказательства). Пусть а) некоторая разностная схе-
ма аппроксимирует задачу (11.2) или задачу (11.3), причём
ψ
k
= O(h
p
);
б) решение разностной схемы
y
k
,
k = 0, 1, . . . , n устойчиво к погрешно-
стям входных данных (т.е. малые погрешности в функции
f из правой
части (11.2) или (11.3) приводят к незначительным изменениям реше-
ния
y
k
). Тогда решение
y
k
,
k = 0, 1, . . . , n разностной схемы сходится к
решению
u(x
k
) дифференциального уравнения при h
→ 0, и имеем место
следующая оценка погрешности
δ
h
= O(h
p
).
Мы не будем останавливаться на вопросе устойчивости решения
y
k
,
k = 0, 1, . . . , n разностных схем. Отметим лишь, что все рассматривае-
мые далее методы устойчивы.
Численное решение дифференциальных уравнений
92
Метод Эйлера.
Заменим в задаче Коши (11.2) производную правой
разделённой разностью
y
k+1
− y
k
h
= f (x
k
, y
k
), k = 0, 1, 2, . . . , n
− 1, y
0
= u
a
.
(11.5)
Решение этой системы уравнений находится явным образом по рекуррент-
ной формуле
y
k+1
= y
k
+ hf (x
k
, y
k
), k = 0, 1, . . . , n
− 1, y
0
= u
a
.
Геометрическая интерпретация метода состоит в замене точного реше-
ния ломаной (рис. 11.1). При этом угол наклона отрезка ломаной при
x
∈ [x
k
, x
k+1
] совпадает с углом наклона касательной к графику точного
решения
u(x) в точке (x
k
, u
k
).
x
0
= a x
1
x
2
b = x
n
u(x)
y
k+1
−y
k
h
= f (x
k
, y
k
)
y
k+1
−y
k
h
= f (x
k+1
, y
k+1
)
y
k+1
−y
k
h
=
f (x
k
,y
k
)+f (x
k+1
,y
k+1
)
2
y
0
y
1
y
2
y
n
Рис. 11.1: Геометрическая интерпретация метода Эйлера. Справа три мо-
дификации метода для разных способов выбора углов наклона отрезков
ломаной.
Найдём погрешность аппроксимации аппроксимации. Для этого нужно
в разностное уравнение (11.5) подставить точное решение
u(x). Применим
следующий метод. Зафиксируем узел
x
k
и выразим значение точного ре-
шения
u(x) в соседних узлах x
k
±1
, x
k
±2
, . . . в виде разложения в ряд Тейло-
ра в окрестности точки
x
k
. В нашем случае в разностное выражение (11.5)
входят только два соседних узла
x
k
и
x
k+1
. Представим
u
k+1
= u(x
k+1
) в
виде ряда Тейлора
u
k+1
= u
k
+ hu
0
k+1
+ O(h
2
).
Подставим точное решение
u
k
и
u
k+1
дифференциального уравнения (11.2)
Численное решение дифференциальных уравнений
93
в разностное уравнение (11.5) вместо
y
k
и
y
k+1
. Вычитая из правой части
левую, получим невязку
ψ
k
=
u
k+1
− u
k
h
− f(x
k
, y
k
) = u
0
k
+ O(h)
− f(x
k
, y
k
).
Из исходного уравнения (11.2) получаем, что
u
0
k
= f (x
k
, y
k
). Отсюда сле-
дует, что
ψ
n
= O(h).
Интуитивно ясно, что качество приближения метода Эйлера определя-
ется правильным выбором углов наклона отрезков ломаной. Можно рас-
смотреть модификацию метода с выбором угла наклона таким же, как у
касательной к графику
u(x) в точке (x
k+1
, u
k+1
). Ясно, что в этом случае
также будет
ψ
k
= O(h).
Лучшее приближение будет достигнуто, если угол наклона отрезков ло-
манной будет заключен между углом касательной к
u(x) в точках (x
k
, u
k
)
и
(x
k+1
, u
k+1
). Так мы приходим к симметричной схеме.
Симметричная схема.
y
k+1
− y
k
h
=
1
2
(f (x
k
, y
k
) + f (x
k+1
, y
k+1
)), k = 1, 2, . . . n
−1, y
0
= u
0
. (11.6)
Данный метод более сложен в реализации, чем метод Эйлера (11.5), так
как новое значение
y
k+1
определяется по найденному ранее
y
k
путём ре-
шения уравнения
y
k+1
− 0,5hf(x
k+1
, y
k+1
) = F
k
,
где
F
k
= y
n
+ 0,5hf (x
k
, y
k
). По этой причине метод называется неявным.
Преимуществом метода (11.6) по сравнению с (11.5) является более высо-
кий порядок точности.
Для невязки
ψ
k
=
u
k+1
− u
k
h
−
1
2
(u(x
k
, y
k
) + u(x
k+1
, y
k+1
))
справедливо разложение
ψ
k
= u
0
k
+
h
2
u
00
k
+ O(h
2
)
−
1
2
(u
0
k
+ u
00
k+1
) = u
0
k
+
h
2
u
00
k
−
1
2
(u
0
k
+ u
0
k
+ hu
00
k
+ O(h
2
)),
т.е.
ψ
k
= O(h
2
). Таким образом, метод (11.6) имеет второй порядок ап-
проксимации, а, следовательно, и второй порядок точности.
Численное решение дифференциальных уравнений
94
Метод предиктор–корректор.
Предположим, что приближённое зна-
чение
y
k
решения исходной задачи в точке
x = x
k
уже известно. Для
нахождения
y
k+1
= y(x
k+1
) поступим следующим образом. Сначала, ис-
пользуя схему Эйлера
y
k+1/2
− y
k
0,5h
= f (x
k
, y
k
),
(11.7)
вычислим промежуточное значение
y
k+1/2
, а затем воспользуемся разност-
ным уравнением
y
k+1
− y
k
h
= f (x
k
+ 0,5h, y
k+1/2
),
(11.8)
из которого явным образом найдём искомое значение
y
k+1
.
Для исследования невязки подставим промежуточное значение
y
k+1/2
=
y
k
+ 0,5hf
k
, где
f
n
= f (x
n
, y
n
), в уравнение (11.8). Тогда получим разност-
ное уравнение
y
k+1
− y
k
h
= f (x
k
+ 0,5h, y
k
+ 0,5hf
k
),
(11.9)
невязка которого равна
ψ
k
=
u
k+1
− u
k
h
− f(x
k
+ 0,5h, u
k
+ 0,5hf
k
).
Имеем
u
k+1
− u
k
h
= u
0
k
+ 0,5hu
00
k
+ O(h
2
).
Используя формулу Тейлора для функции нескольких переменных
, по-
лучим окрестности точки
(x
k
, u
k
)
f (x
k
+ 0,5h, u
k
+ 0,5hf (x
k
, u
k
)) = f (x
k
, u
k
)+
+ 0,5h
∂f (x
k
, u
k
)
∂x
+ f (x
k
, u
k
)
∂f (x
k
, u
k
)
∂u
+ O(h
2
) =
= f (x
k
, u
k
) + 0,5hu
00
k
+ O(h
2
),
8
Теорема. Пусть функция
u(x
1
, x
2
, . . . , x
m
) непрерывно дифференцируема (n
− 1)
раз в
ε-окрестности точки M
0
(
◦
x
1
,
◦
x
2
, . . . ,
◦
x
m
) и n раз дифференцируема в самой точке
M
0
. Тогда для любой точки
M из указанной ε-окрестности M
0
справедлива следующая
формула
u(M ) = u(M
0
) +
1
1!
du
M
0
+
1
2!
d
2
u
M
0
+ . . . +
1
n!
d
n
u
M
0
+ O(ρ
n+1
),
где
ρ = ρ(M, M
0
) — расстояние между точками M и M
0
.
Численное решение дифференциальных уравнений
95
так как в силу (11.2) справедливо равенство
u
0
= f (x, u) и, следовательно,
u
00
=
∂f
∂x
+ f
∂f
∂u
(производная функции двух переменных)
.
Таким образом, метод (11.9) имеет второй порядок погрешности аппрок-
симации,
ψ
k
= O(h
2
), и в отличии от (11.6) является явным.
Реализация метода (11.9) в виде двух этапов (11.7), (11.8) называется
методом предиктор–корректор (предсказывающе-исправляющим
), по-
скольку на первом этапе приближённое значение предсказывается с невы-
сокой точностью
O(h), а на втором этапе это предсказанное значение ис-
правляется, так что результирующая погрешность имеет второй порядок
по
h.
Тот же самый метод можно реализовать несколько иначе. А именно,
сначала вычислим последовательно функции
p
1
= f (x
k
, y
k
),
p
2
= f (x
n
+ 0,5h, y
k
+ 0,5hp
1
),
а затем найдём
y
k+1
из уравнения
(y
k+1
− y
k
)/h = p
2
.
Такая форма реализации метода предиктор–корректор называется ме-
тодом Рунге–Кутта. Поскольку требуется вычислить две промежуточные
функции
p
1
и
p
2
, данный метод относится в двухэтапным методам.
Метод Рунге–Кутта.
Явный
m-этапный метод Рунге-Кутта состоит в сле-
дующем. Пусть решение
y
k
= y(x
k
) уже известно. Задаются числовые коэффициенты
a
i
, b
ij
, i = 2, 3, . . . , m, j = 1, 2, . . . , m
− 1,
σ
i
,
i = 1, 2, . . . , m,
и последовательно вычисляются функции
p
1
= f (x
k
, y
k
),
p
2
= f (x
k
+ a
2
h, y
k
+ b
21
hp
1
),
p
3
= f (x
k
+ a
3
h, y
k
+ b
31
hp
1
+ b
32
hp
2
),
· · ·
p
m
= f (x
k
+ a
m
h, y
k
+ b
m1
hp
1
+ b
m2
hp
2
+ . . . + b
m,m−1
hp
m−1
).
Затем из формулы
y
k+1
− y
k
h
=
m
X
i=1
σ
i
p
i
(11.10)
находится новое значение
y
k+1
= y(x
k+1
).
9
От англ. predict – предсказывать, correct – исправлять
10
М´
артин Вильг´
ельм К´
утта — немецкий физик и математик (1867–1944).