book1989
.pdfчае одного уравнения, а именно, метод имеет квадратичную сходи мость, если начальное приближение выбрано достаточно хорошо.
|
Приведем без доказательства одну из теорем |
о |
сходимости |
метода |
Нью |
||||
тона. |
|
|
вещественных |
векторов с |
нормой |
1|х|| = |
|||
/ |
Пусть Ет — множество m-мерных |
||||||||
т |
у/, |
|
|
|
|
|
|
|
|
= [ |
|
,|И!| — норма матрицы А, |
подчиненная |
|
данной |
норме вектора. |
|||
V i= 1 . |
|
|
|
|
|
|
|
||
Обозначим |
|
|
|
|
|
|
|
||
|
|
Ur (x°) = {.<<=Ет: ||х—х°\\ < |
г} |
|
|
|
|
||
и предположим, что в шаре UT(x°) функции fi(x), |
i = 1, 2, |
. . . , |
т, непрерывно |
||||||
дифференцируемы. |
что в |
UT(x°) матрица F'(x) |
удовлетворяет ус |
||||||
|
Т е о р е м а 2. Предположим, |
||||||||
ловию Липшица с постоянной L, |
т. е. |
|
|
|
|
|
|
||
НЕ' (х1)—F' (л?) || L\\x'—дг2||
для любых х’, x2e U r(x°). Пусть в UT(xA) матрица (Е '(х))-1 существует, причем элементы ее непрерывны и
1!(F’(x))-4]^ M .
Если начальное приближение х° таково, что ||.F(xo) И=^т) и
причем
оо
Мц 2 92*-1 < г, k—{'
то система уравнений (2 ) имеет решение x . e U r(x0), к которому сходится метод
Ньютона (14). Оценка погрешности дается неравенством
, |
— -X» | 5С |
(А 1 |
II к |
Мт\ ----- -k . |
|
|
|
1-<?2 |
Доказательство теоремы 2 можно найти в [42].
П р и м е р 4. Модифицированный метод Ньютона имеет вид
F'(x°) (xk+i—xh)+F(xk) =0 |
(15) |
и обладает линейной сходимостью. Упрощение в численной реали зации по сравнению с обычным методом Ньютона состоит в том, что матрицу F'(х) надо обращать не на каждой итерации, а лишь один раз. Возможно циклическое применение модифицированного метода Ньютона, когда F'(x) обращается через определенное чис ло итераций.
П р и м е р 5. Метод Ньютона с параметром имеет вид
F' (.хk) |
+ F (xk) = 0. |
(16) |
THI
Рассмотренные до сих пор методы являлись линейными отно сительно новой итерации я*41. Возможны и нелинейные методы, ког
211
да для вычисления xh+i приходится решать нелинейные системы уравнений. Приведем примеры таких методов.
П р и м е р 6. Нелинейный метод Якоби для системы |
(1) имеет |
||
вид |
|
|
|
k k+l |
k |
40 = 0, |
(17) |
/; (4, 4, ■ , Xj—i, Xi |
t Xij-iy . |
||
i = 1, 2 , . . . , m.
Здесь для отыскания xk+l необходимо решить m независимых скалярных уравнений. Для решения скалярного уравнения можно применить какой-либо из итерационных методов, рассмотренных в § 1, причем не обязательно применять один и тот же метод для всех уравнений.
П р и м е р |
7. Нелинейный метод Зейделя |
состоит в последова |
|||||
тельном решении уравнений |
|
|
|
|
|
||
|
f.(yk+l |
xk+l |
4 +l уk |
...» |
4л n |
(18) |
|
|
11 |
, л2 , . . . , л1 , л4+1, |
Лт)— |
и |
|||
|
|
гаЬ. 1 |
i= 1, 2, . . . , m. |
|
|
|
|
относительно переменной х,- , |
|
методы, |
когда |
||||
Большое |
распространение |
получили |
гибридные |
||||
внешние итерации осуществляются одним методом, а внутренние — другим. При этом число внутренних итераций может быть фикси рованным и не очень большим, так что внутренние итерации не до водятся до сходимости. В результате получается некоторый новый метод, сочетающий свойства исходных методов. Приведем приме ры таких методов.
П р и м е р 8. Внешние итерации — по Зейделю |
и внутренние — |
по Ньютону. Здесь в качестве основной (внешней) итерации выби |
|
рается нелинейный метод Зейделя (18), а для |
нахождения 4 +l |
используется метод Ньютона. Обозначим г/,= х*+1. Тогда итерации определяются следующим образом:
4 + 1 |
, |
r * + l |
y .S |
4 ч / Л - 1 |
S\ . |
|
|
(,-П , Х2 |
. . . , Л,-!, |
уи л4+1, |
. . . , Хщ) [tji |
— Уц -f- |
|
||
|
|
|
+ П4+\ |
„s-и |
|
||
|
|
|
, 44,4,4».........4.) = о, |
||||
|
|
5 = 0 , 1 , . . . , / , |
4 = 4 , |
yl4 = x Y \ |
(19) |
||
|
|
|
i'=l, 2........m. |
|
|
|
|
Здесь индексом s обозначен номер внутренней итерации.
Иногда в (19) делают всего одну внутреннюю итерацию, пола
гая 1 = 0, |
yai= x t,y \ |
|
— 4 fl. Тогда |
приходят |
к следующему |
итера |
||||||
ционному методу: |
|
|
|
|
|
|
|
|
|
|||
4 li. (vk+l |
, |
rk+l |
• • • |
t |
4+i |
4 |
xk |
, x m) (X; +1 — |
X;) |
|
||
dxi |
Л2 » |
|
Li Xi, |
Л^-i, . . • |
|
|||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
4- f- (xkH |
• • t |
v - +1 |
4 4 |
|
0, |
k = 0, 1 , . . . |
(20) |
|||||
1 |
/ l \Л 1 |
* • |
|
1» |
X/ > X l + l t . . . , 4 . ) = |
|||||||
212 |
|
|
|
|
|
|
|
|
|
|
|
|
В частности, при т = 2 метод (20) принимает вид
|
dfl(xJ' ^ (4Л- 4) + h(4,4) = 0, |
||
|
dXl |
|
|
|
|
|
(21) |
|
(*Г\ 4 |
(х!н - х к) + h (4 +1, **) = |
0. |
|
иХ% |
|
|
П р и м е р |
9. Внешние итерации— по Ньютону и внутренние — |
||
по Зейделю. Запишем метод Ньютона для системы (2) в виде |
|||
|
F' (xft) (x*+1—хк) +F {хк) =0, |
(22) |
|
где F'(xk) = |
(йц), ау = |
dfi (*к) |
m. Для решения |
---------, i, / = 1, 2, .. ., |
|||
|
|
dxj |
|
системы линейных уравнений (22) воспользуемся методом Зейделя. Напомним (см. § 1 гл. 2), что для линейной системы
Aw+F = Q |
(23) |
метод Зейделя строится следующим образом. Матрица А представ ляется в виде суммы А = Л_+/)-|-Л+, где матрицы Л_, А+, D соот ветственно нижняя треугольная, верхняя треугольная и диагональ ная. Итерации метода Зейделя строятся по правилу
{A_+D)wa+l+ A +wa+F = 0, s = 0, |
1, . . . , / , |
(24) |
|
и система (24) решается |
путем обращения |
нижней |
треугольной |
матрицы A-A-D. |
надо положить A = F'(xk), вычислить по |
||
В случае системы (22) |
|||
следовательно векторы ws согласно (24), начиная с w° = 0, и поло жить wl+l= xk+l—хк, так что хк+>= xk-\-wI+l.
Заметим, что итерации по Зейделю можно осуществлять и от
носительно вектора х',+1. |
|
|||
Пусть в (24) |
совершается только одна итерация, т. е. 1 = 0. Тог |
|||
да, учитывая, что ш°= 0, wl= xk+l—хк, получим метод |
|
|||
|
|
|
(А_+£>) (xk+l—хк) + F (**) = 0, |
(25) |
где Л _ + 0 — «нижняя треугольная» часть матрицы Якоби |
(11), |
|||
вычисленной при х= хк. |
|
|||
В частности, при т = 2 метод (25) принимает вид |
|
|||
|
|
и |
(хки 4 ) (хк+1- 4 ) + h (хк, Хк) = о, |
|
|
|
|
(26) |
|
|
|
|
|
|
дх\ |
(А4)(4п-А +-&-(4 , А ( 4 м - А +и А А = |
о. |
||
|
|
дх2 |
|
|
Сопоставление (21) и (26) показывает, что методы, рассмот ренные в двух последних примерах, не совпадают.
213
Г Л А В А 6
ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ ЗАДАЧИ КОШИ ДЛЯ ОБЫКНОВЕННЫХ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ
§1. Исходная задача и примеры численных методов ее решения
1.Постановка исходной задачи. Будем рассматривать задачу
Коши для системы |
обыкновенных |
дифференциальных уравнений |
|||||
at |
= |
|
t > О, |
u(0) = |
u<°> |
(1) |
|
|
|
|
|
|
|
|
|
или, подробнее, |
|
|
|
|
|
|
|
= fi if, |
«1, «г, |
• • • |
, Um), |
t > |
0, i = |
1, 2, . . . , m, |
(2) |
at |
|
|
|
|
|
|
|
щ (0) = |
u f\ |
i — 1,2, |
... , m. |
|
(3) |
||
Хорошо известны условия, гарантирующие существование и единственность решения задачи Коши (см. [39, с. 49]). Предполо жим, что функции fu t =l , 2, . . . , т, непрерывны по всем аргумен там в замкнутой области
D = { \t \^ a , |
| ut— ы10>| sg b, |
i = 1,2, . . . , т). |
|||
Из непрерывности |
функций |
следует их ограниченность, т. е. |
|||
существование константы |
0 такой, что всюду в D выполняются |
||||
неравенства |Д| |
i= 1, 2, . . . , т. |
|
удовлетворяют |
||
Предположим, кроме того, что в D функции |
|||||
условию Липшица по аргументам и{, и2, ... , ит, т. е. |
|
||||
fi (t, Щ, U2, . .. , и'т) — fi (t, Uu и2, |
, |
ит) I ^ |
|
||
|
L {| и[ —и[ I + I и2 — и, I -f ... + I ит— ит|} |
||||
для любых точек (/, ии . . . , ит) |
и (t, ии и2, . . . , ит) |
области D. |
|||
Если выполнены сформулированные выше предположения, то |
|||||
существует единственное решение |
|
|
|||
ui= ul(t), u2 = u2(t), . |
. . , um= u m(t) |
|
|||
системы (2), определенное |
при |
|/| ^ / 0 = min(a, b/М) и принимаю |
|||
щее при ^ = 0 заданные начальные значения (3).
При исследовании численных методов для задачи Коши будем заранее предполагать, что ее решение существует, единственно и обладает необходимыми свойствами гладкости.
2. Примеры численных методов. Существуют две группы чис ленных методов решения задачи Коши: многошаговые разностные методы и методы Рунге — Кутта. Приведем примеры и поясним ос новные понятия, возникающие при использовании численных мето дов. Для простоты изложения будем рассматривать сейчас одно
214
уравнение |
|
|
|
- ± = f(t,u), |
/ > 0, |
и (0) — ии. |
(4) |
a t |
|
|
|
Введем по переменному t равномерную сетку с шагом т> 0 , т. е. |
|||
рассмотрим множество точек |
|
|
|
Шт= { t n = tix, п = 0, |
1,2,...}. |
|
|
Будем обозначать через u(t) точное решение задачи |
(4), а че |
||
рез y n = y ( t n) — приближенное |
решение. Заметим, что приближен |
||
ное решение является сеточной функцией, т. е. определено только в точках сетки шг.
П р и м е р 1. Метод Эйлера. Уравнение (4) заменяется разност
ным уравнением |
|
■П¥1~ Уп — f (t,i, уф = 0, я = 0, 1, 2, . .. , у0 = и,. |
(5) |
т |
|
Решение этого уравнения находится явным образом по рекур рентной формуле
уп+1 = Уп+т/(t n, у„), п = 0, 1, . . . , у0—ий.
При использовании приближенных методов основным является вопрос о сходимости. Понятие сходимости приближенного метода можно сформулировать по-разному. Применительно к разностным методам, к которым относится и метод Эйлера (5), наибольшее рас пространение получило понятие сходимости при т—>-0. Оно означает
следующее. Фиксируем |
точку |
t |
и |
построим |
последовательность |
|
сеток шт таких, что т—>-0 и tn = m = t |
(тогда необходимо п->-оо). Го |
|||||
ворят, что метод (5) сходится |
в |
точке t, если |
\уп— |
при |
||
т^О, tn = t. |
отрезке |
(0, Т\, если он сходится |
в каждой |
|||
Метод сходится на |
||||||
точке 1<= (0, Т]. |
|
|
|
|
|
|
Говорят, что метод имеет р-й порядок точности, если существует число р > 0 такое, что | уп—u(tn) | = 0 (%“) при т—^0.
Получггм уравнение, которому удовлетворяет погрешность ме
тода z„= уп—u(tn). Подставляя y n = zn-\-un в (5), получим |
|
|||||||
Znt-L ~~гп |
|
Г,, |
| |
ч |
' |
игч-1~ип |
• |
(6) |
т |
—/ viii Мц |
“р %п) |
Т |
|||||
|
|
|
|
|
|
|
||
Правую часть уравнения |
(6) можно представить в виде суммы |
|
||||||
где |
|
|
|
|
|
|
|
|
.1 (!) |
|
Un+l |
Utl |
, |
г /, |
ч |
|
|
'ри — ---------------- Ь f (tn, ип), |
|
|
||||||
|
|
т |
|
|
|
|
|
|
фг?’ = |
|
f (tn, Un + |
гп) — f (tn, Un). |
|
|
|||
Функция фгё1) называется невязкой или погрешностью аппрок симации разностного уравнения (5) на решении исходного урав нения (4). Видно, что невязка представляет собой результат подста-
215
новки точного решения u = u(t) в левую часть разностного уравне ния (5). Если бы приближенное решение уп совпадало с точным u(tn), то невязка равнялась бы нулю. Говорят, что разностный ме тод аппроксимирует исходное дифференциальное уравнение, если
при т-»-0. Разностный метод имеет р-й порядок аппрокси
мации, если г|)^ = 0 (т р). В дальнейшем будет показано, что при очень общих предположениях порядок точности разностного мето да совпадает с порядком аппроксимации.
Функция
фп’= f (In, ип + zn) — f (in, ип)
обращается в нуль, если правая часть f не зависит от решения и.
В общем случае |
> пропорциональна погрешности zn, так как по |
формуле конечных приращений имеем |
|
ty[n |
= ^ - ( tn ,un + ten)Zn, | 9 | < 1 . |
|
ди |
Порядок аппроксимации метода Эйлера (5) нетрудно найти, используя разложение по формуле Тейлора. Поскольку
= |
+ |
т ) , |
т |
|
|
то в силу уравнения (4) |
|
|
фя1>= — и’ (tn) + / (tn, u„) -f 0 (т) = 0 (т),
т. е. метод Эйлера имеет первый порядок аппроксимации. При вы воде предполагалась ограниченность и" (t).
П р и м е р |
2. Симметричная схема. Уравнение (4) |
заменяется |
разностным уравнением |
|
|
Уп+,~ Уп _ ± |
у (/я> уп) + f {in+i> yn+i)) = Q> n = 0, i t ... , |
у0 = и0ш (7) |
т2
Данный метод более сложен в реализации, чем метод Эйлера (5), так как новое значение yn+i определяется по найденному ра нее Уп путем решения уравнения
Уп+1 0Дт/(^п+1) Уп+1) Fn,
где F„= y„-Г0,5т/(<п, уп). По этой причине метод называется неяв ным. Преимуществом метода (7) по сравнению с (5) является бо лее высокий порядок точности.
Для невязки
ФД = - |
+ ~ ( f |
( U Уп) + |
/ (in,и УпЛ) |
справедливо разложение |
|
|
|
W = - ип ~ 7 и’п + 0 И + |
\ К + |
Un J = |
|
= - |
ип - |
f 'и"п+ \ |
К + ип + хи'п + 0 (**)). |
216
сказывается с невысокой точностью 0 (т), а на втором этапе (9) это предсказанное значение исправляется, так что результирующая погрешность имеет второй порядок по т.
Тот же самый метод (10) можно реализовать несколько иначе. А именно, сначала вычислим последовательно функции
ki = f (tn, у„), k2 = j(tn+ 0,5т, уп+ 0,5т*!),
а затем найдем уп+1 из уравнения (уп+1—yn)li = k2.
Такая форма реализации метода (10) называется методом Рунге — Кутта. Поскольку требуется вычислить две промежуточные функции, /г, и k2, данный метод относится к двухэтапным методам. В следующем параграфе будут рассмотрены более общие т-этап- ные методы Рунге —Кутта, позволяющие получить большую точ ность.
§2. Методы Рунге — Кутта
1.Общая формулировка методов. Семейство методов второго порядка. По-прежнему рассматриваем задачу Коши для одного уравнения
|
^ = f(t,u), |
i > 0, |
и (0) = ы0. |
|
(1) |
|
|
at |
|
|
|
|
|
Явный |
m-этапный метод Рунге —Кутта состоит |
в |
следующем. |
|||
Пусть решение yn = y{tu) |
уже |
известно. Задаются числовые коэф |
||||
фициенты |
а„ b,j, 1= 2, 3, |
..., |
m, \ = 1, |
2, .... т—1, |
о„ |
£=1, 2, ... |
..., т, и последовательно вычисляются функции |
|
|
||||
ki = / (tn, Уп), /?2 = / (tn + а2Т, у,г + Ь21Т&,), К = / (tn + а3х, уп+ b3]Tkt + b32т/е,),
К. = / (tn + а„Д, Уп +bm]%kl + bnl2lk2+ . . . + b.n^xkm-i).
Затем из формулы
^Яг1~ - = |
Oikt |
(2) |
ТI—1
находится новое значение i/n+1 = i/(£n+1).
Коэффициенты аи b,h сн выбираются из соображений точности.
Например, для того чтобы уравнение (2) аппроксимировало исход- |
|||||
|
|
|
/П |
|
|
мое уравнение |
(1), необходимо потребовать |
Отметим, |
|||
|
|
|
1=1 |
|
|
что методы Рунге —Кутта при т > 5 не используются. |
|
т = |
|||
Остановимся |
более подробно на отдельных методах. При |
||||
= 1 получаем схему Эйлера, |
рассмотренную в примере 1 из |
пре |
|||
дыдущего параграфа. При т = 2 получаем |
семейство |
методов |
|
||
k, |
= f(t„, ул), |
kz = f(tn+a2т, |
уп+Ь21т/г,), |
|
(3) |
|
|
|
|
|
|
//„-и =yn +l(o 1k,+Olk2).
Исследуем погрешность аппроксимации методов (3) в зависи мости от выбора параметров. Исключая из последнего уравнения
218
функции k, |
и k2, получаем |
|
|
|
|
|
|
|
|
|
|||
—^--—=■o j (tn, Уп) + Ozf (tn 4~ a.2X, Уп 4" ^21Tf (tn, Уп))- |
(4) |
||||||||||||
|
T |
|
|
|
|
|
|
|
|
|
|
|
|
По определению |
погрешностью |
аппроксимации |
или невязкой |
||||||||||
метода (3) называется выражение |
|
|
|
|
|
|
|
|
|||||
фпг) = — — ---- - + |
Gif (-«. “г) + |
|
o j (tn+ |
a2x, ua4- b2J (tn, un)), |
(5) |
||||||||
|
T |
|
|
|
|
|
|
|
|
|
|
|
|
полученное заменой |
в |
(4) |
приближенного решения уп точным ре |
||||||||||
шением un = u(tn). |
|
|
|
|
|
|
|
|
|
|
|
||
Найдем порядок погрешности аппроксимации в предположении |
|||||||||||||
достаточной гладкости решения u(t) |
и функции f(t, |
и). Для этого |
|||||||||||
разложим все величины, входящие в выражение (5), |
по формуле |
||||||||||||
Тейлора в точке tn. Имеем |
|
|
|
|
|
|
|
|
|
||||
|
|
-n- l ~ Un= и’ (tn) 4 -1 и |
(tn) |
4- О(т2), |
|
|
|
||||||
|
|
|
т |
|
|
|
2 |
|
|
|
|
|
|
f (tn + |
а.2х, tin + |
b2lrfn) = f n + |
a2x |
dt |
+ |
b2lxfn |
+ |
О(т2), |
|
||||
|
|
|
|
|
|
|
|
|
да |
|
|
|
|
где f„= f(tn, |
un), |
— — — |
(tn, un). Далее, согласно уравнению |
(1), |
|||||||||
получим |
|
ди |
|
да |
|
|
|
|
|
|
|
|
|
|
|
,, |
d f |
. d f |
, |
|
d f |
, |
d f t |
|
|
|
|
|
|
|
|
|
|
|
|||||||
|
|
u — — |
H— - u |
|
= — + — /• |
|
|
|
|||||
Поэтому |
|
|
|
d t |
d a |
|
|
d t |
|
d u |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Ф„f = — u’n4- (ay + cr2) fn 4- |
|
|
|
|
|
|
|
|
|
||||
|
+ T |
|
|
|
dfn |
|
|
|
|
4-0(т2). |
(6) |
||
|
(a2b2l — 0,5) fn ^ - + (оая2- ° , 5)^ ] |
||||||||||||
|
|
|
|
|
да |
|
|
|
|
|
|
|
|
Отсюда видно, что методы (3) имеют первый порядок аппрок симации, если Oi4-a2= 1.
Если же дополнительно потребовать o2a2 = o2b2I = 0,5, то полу чим методы второго порядка аппроксимации. Таким образом, име ется однопараметрическое семейство двухэтапных методов Рунге — Кутта второго порядка аппроксимации. Это семейство методов можно записать в виде
— — — = (1 — о) f (tn, Уп) + of (tn+ ах, у„. + axf (tn, уп)), (7) |
||
х |
|
|
где ста = 0,5. |
а = 0,5 получаем метод, |
рассмотренный |
В частности, при о=1, |
||
в примере 3 предыдущего |
параграфа. При с = 0,5, |
а—1 получаем |
другой метод второго порядка: |
|
|
ki = f(tn, Уг), k2= f(tn+ х, уп+xk,), yn+i= y,+Q,5x(kl+k2).
219
Двухэтапных методов третьего порядка аппроксимации не суще ствует. Чтобы убедиться в этом, достаточно рассмотреть уравне ние и' = и. Для него двухэтапный метод Рунге —Кутта (7) прини мает вид
|
----------= (14- хаа) уп |
и погрешность аппроксимации равна |
|
|
г^1’ = --- —— — 4- (1 4- таа) ип. |
|
т |
Разлагая |
по формуле Тейлора и учитывая, что и"' = и'' = и' = |
= и, получим |
|
= т (аа — 0,5) и — — и (1„4- 9т), 0 ^ 0 < 1.
6
Отсюда видно, что наивысший достижимый порядок аппроксима ции равен двум.
Приведем примеры методов Рунге —Кутта более высокого по рядка аппроксимации.
Метод третьего порядка:
K = f ( t n,y n), |
k2 = f(tn + J , Уп + |
, |
||||
К = f(tn 4- т, уп—- т/4 4- 2xk2), |
Упи — Уп = ± ( k l +4k2+ k3). |
|||||
Метод третьего порядка: |
|
|
О |
|
||
|
|
|
\kx |
|||
k i = f {tn, Уп), |
k2 = f ^ n + — , y n - |
|||||
з |
||||||
+ |
Уп + Щ |
, |
УП+1 ~ - = j - ( k l +3k3). |
|||
\ |
3 |
3 |
1 |
т |
4 |
|
Метод четвертого порядка: |
|
|
Xkj |
|||
h = f(tn,y,i), |
k2 = f[tn + — , уп . |
|||||
2 |
||||||
: / (tn 4- j |
, yn4- — 'j, |
k4 = f (tn 4- T, yn 4- тк3), |
||||
Уп¥у~ [,п= -i (ky + ?.k24- 2k3 + k,).
т6
Метод четвертого порядка: |
|
|
||
^1 — / (^n, */n)> ^2 —- / |
4~ — , |
Уп 4------1 |
||
h |
= |
f ( t n A - - , |
y«.-\- |
, |
fe4 = f(tn+ T, (/„ 4- xfe4 — 2xk24- 2xk3), |
||||
y.^ |
- |
y« = ± {kl + 4k3 + k4) |
||
|
x |
6 |
|
|
220