Методы вычислительной математики

Видно, что с уменьшением шага интегрирования h численное решение приближается к точному. Противоположную картину можно наблюдать при увеличении шага интегрирования.

На рис. 8.4 приведены результаты расчетов той же задачи с большими шагами интегрирования. При шаге h = 1 величины yk уже на втором шаге интегрирования обращаются в 0 и в дальнейшем не изменяются (см. рис. 8.4, а).

Рис. 8.4. Результаты численного решения разностного уравнения с большими шагами интегрирования: h = 1 (а), h = 1,5 (б), h = 2 (в), h = 2,1 (г); — – точное решение, —о— – численное решение

Шаг интегрирования h = 1,5, дает отрицательные значения yk для некоторых значений xk, что противоречит точному решению задачи (см. рис. 8.4, б).

171

Следует отметить, что даже в этом случае можно наблюдать сходимость численного решения к точному при большом числе шагов. Рис. 8.4, в показывает, что при h = 2 численное решение дает значения yk, осциллирующие возле действительного значения. При дальнейшем увеличении шага получается расходящееся решение (см. рис. 8.4, г). Представленные результаты говорят о необходимости проведения анализа результатов расчетов, полученных с использованием ЭВМ, с целью проверки сходимости численных решений к точным.

8.4. Метод Рунге–Кутты1 второго порядка аппроксимации

Для построения разностной схемы используется разложение решения обыкновенного дифференциального ураынения (8.1) в ряд Тейлора:

y(xk+1 )= y(xk )+ y′x (xk )h + y′xx′ (xk )h2 2 +….

Вторая производная в этом разложении заменяется выражением

наибольшей точности записанного выражения. Для дальнейших выкладок производится замена величины ~y разложением в ряд Тейлора:

Вводятся обозначения

α = h2 x, β =1 − h2 x, γ = xh, δ = f (xk , yk ) xh ,

которые позволяют записать предыдущее выражение в форме yk+1 = yk + h[βf (xk , yk )+ αf (xk + γh, yk + δh)].

Введенные коэффициенты зависят от x и могут быть определены через коэффициент α, который в этом случае играет роль параметра,

1 Кутта Мартин Вильгельм [3.11.1867 – 25.12.1944] – немецкий физик и математик, был профессором Высшей технической школы Штудгарта.

172

β =1 − α, γ =12α, δ = f (xk , yk )2α.

Окончательно схема Рунге–Кутты принимает форму выражения

Та же схема в форме разностного аналога уравнения (8.1) имеет вид

(yk +1 − yk )h = (1 − α)f (xk , yk ) + αf (xk + h2α, yk + hf (xk , yk )2α).

При α = 0 получается, как частный случай, известная схема Эйлера: yk+1 = yk + hf (xk , yk ).

При α = 1 выражение (8.11) записывается в форме

yk+1 = yk + hf (xk + h2, yk + hf (xk , yk )2).

В этом случае проведение расчетов на очередном шаге интегрирования можно рассматривать как последовательность операций (рис. 8.5).

1. Рассчитывается значение

K1 = f (xk , yk )

производной y′x искомого решения в узле xk разностной сетки. 2. Вычисляется выражение

ykЭ+1/ 2 = yk + hf (xk , yk )2 = yk + hK1 2 ,

представляющее собой полушаг интегрирования по схеме Эйлера, то есть определяется приближенное значение искомой функции в точке xk + h2.

y

Рис. 8.5. Схема интегрирования обыкновенного дифференциального уравнения методом Рунге–Кутты второго порядка с параметром α = 1

3. Определяется приближенное значение

173

производной y′x для той же промежуточной точки xk + h2.

4. Вычисляется уточненное значение искомой функции в конечной точке xk+1 = xk + h всего шага по схеме Эйлера с вычисленным на предыдущем ша-

ге значением производной K2,

yk+1 = yk + hK2 .

Геометрические построения (см. рис. 8.5) показывают, что получаемая с помощью такой последовательности точка yk+1 лежит ближе к истинному решению, чем вычисляемая по схеме Эйлера, то есть решение, получаемое методом Рунге–Кутты, оказывается более точным.

На рис. 8.6 показана геометрическая интерпретация разностной схемы1, получающейся из выражения (8.11) при α =12 ,

yk+1 = yk + h[f (xk , yk )+ f (xk + h, yk + hf (xk , yk ))]2 .

y

yk+1 K2

производной y′x искомого решения в узле xk разностной сетки.

2. Выполняется полный шаг метода Эйлера для определения приближенного значения искомой функции на конце отрезка интегрирования xk +1 = xk + h:

ykЭ+1/ 2 = yk + hf (xk , yk )= yk + hK1 .

1 Получающееся выражение называется схемой (методом) Эйлера–Коши.

174

5. Вычисляется значение искомой функции в конечной точке xk+1 всего шага по схеме Эйлера с усредненным значением производной:

yk+1 = yk + h(K1 + K2 )2 .

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

Пример 8.5. Решить уравнение y′ = −y, y(0) =1 методом Рунге–Кутты вто-

K2 = f (xk + h, ykЭ+1 ) = −yk (1 − h);; (K1 + K2 )/ 2 = −yk (2 − h)/ 2;

yk+1 = yk − hyk (2 − h)2 = yk [(h −1)2 +1]2 .

Строится последовательность значений искомой функции, y0 = y(0) =1,

y1 = y0 [(h −1)2 +1]2= [(h −1)2 +1]2 , y2 = y1 [(h −1)2 +1]2 = {[(h −1)2 +1]2}2 , y3 = y2 [(h −1)2 +1]2 = {[(h −1)2 +1]2}3 , ...,

175

ym = {[(h −1)2 +1]2}m .

Рзультаты получаемого численного решения для значения аргумента x = 10 при различных шагах интегрирования приведены в табл. 8.2. Три верные значащие цифры получены теперь для шага h = 0,01.

Таблица 8.2 Результаты численного решения ym методом Рунге–Кутты второго порядка

дифференциального уравнения y′ = −y с начальным условием y(0) =1

Для оценки погрешности аппроксимации уравнения (8.1) разностной схемой метода Рунге–Кутты точное решение подставляется в разностный аналог исходного дифференциального уравнения и вычисляется невязка:

ψk = [y(xk +1 )− y(xk )]h − (1 − α)f (xk , y(xk ))−

−αf (xk + h2α, y(xk )+ hf (xk , y(xk ))2α).

Разложения функций в ряды Тейлора

y(xk+1 )= y(xk )+ y′(xk )h + y′′(xk )h2 2 + O(h3 ),

f(xk + h2α, y(xk ) + hf (xk , y(xk ))2α)=

=f (xk , y(xk ))+ hfx′(xk , y(xk ))2α + O(h2 )

подставляются в полученное выражение:

ψk = [y(xk ) + y′x (xk )h + y′xx′ (xk )h2 2 +O(h3 )− y(xk )]h − (1 − α)f (xk , y(xk ))−

−α{f (xk , y(xk ))+ hfx′(xk , y(xk ))2α + O(h2 )}=

=y′(xk )− f (xk , y(xk ))+h[y′′(xk )− fx′(xk , y(xk ))]2 +O(h2 ).

Учитывая уравнение (8.1), а также выражение для производной,

y′′(xk ) = [y′(xk )]′ = f x′(xk , y(xk )),

можно получить, что ψk = O(h2 ), то есть метод Рунге–Кутты, независимо от значения параметра α имеет погрешность аппроксимации второго порядка.

176

8.5. Методы Рунге–Кутты третьего и четвертого порядков аппроксимации

Ниже приводятся две схемы Рунге–Кутты, предназначенные для численного решения обыкновенных дифференциальных уравнений первого порядка и имеющие погрешность аппроксимации третьего порядка [22]:

и две схемы Рунге–Кутты, имеющие погрешность аппроксимации четвертого порядка:

Пример 8.6. Решить методом Рунге–Кутты четвертого порядка уравнение y′ = −y, y(0) =1.

Правая часть дифференциального уравнения имеет вид f (x, y) = −y , поэтому схема метода Рунге–Кутты четвертого порядка представляется следующим образом:

K1 = −yk ;

K 2 = −(yk + hK 1 2) = − yk (1 − h2);

K 3 = −(yk + hK 2 2) = − yk [1 − h(1 − h2)2]; K 4 = −(yk + hK 3 ) = − yk {1 − h[1 − h(1 − h2)2]};

yk+1 = yk + h(K1 + 2K2 + 2K3 + K4 )6 = yk [1 + (h4 − 4h3 +12h2 − 24h)24].

Строится последовательность значений искомой функции,

177

y0 = y(0) =1,

y1 = y0 [1 + (h4 − 4h3 +12h2 − 24h)24]= [1 + (h4 − 4h3 +12h2 − 24h)24], y2 = y1[1 + (h4 − 4h3 +12h2 − 24h)24]= [1 + (h4 − 4h3 +12h2 − 24h)24]2 , y3 = y2 [1 + (h4 − 4h3 +12h2 − 24h)24]= [1 + (h4 − 4h3 +12h2 − 24h)24]3,

……………………………………………………………………………..., ym = [1 + (h4 − 4h3 +12h2 − 24h)24]m.

Результаты численного решения для значения аргумента x = 10 при различных шагах интегрирования приведены в табл. 8.3. Три верные значащие цифры получены для шага h = 0,25.

Таблица 8.3 Результаты численного решения ym методом Рунге–Кутты четвертого порядка

дифференциального уравнения y′x = −y с начальным условием y(0)=1

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

при том же шаге интегрирования или получить ту же точность при большем шаге интегрирования, и, следовательно, меньшем числе шагов, то есть приво-

дит к снижению требуемых ресурсов ЭВМ.

8.6. Метод Адамса1

Предполагается, что при решении задачи (8.1) для четырех последовательных точек xk−3 , xk−2 , xk −1, xk известны приближенные значения искомой функ-

1 Адамс Джон Кауч [5.6.1819 – 21.1.1892] – английский астроном и математик, член Лондонского королевского общества. В 1861 году стал директором Кембриджской астрономической обсерватории. В 1864 году был избран иностранным членом-корреспондентом Петербургской академии наук.

178

ции yk−3 , yk−2 , yk−1, yk , а значит, можно подсчитать соответствующие значения правой части дифференциального уравнения,

где

F(xk , xk−1 )= [F(xk ) − F(xk−1 )](xk − xk−1 ),

F(xk , xk −1, xk−2 )= [F(xk , xk−1 )− F(xk−1, xk−2 )](xk − xk−2 ),

F(xk , xk−1, xk−2 , xk−3 )= [F(xk , xk−1, xk −2 )− F(xk−1, xk −2 , xk−3 )](xk − xk−3 )

– разделенные разности, являющиеся разностными аналогами соответствующих производных функции F(x).

Уравнение (8.1) интегрируется на сегменте [xk ,xk+1 ],

xk +1

y(xk+1 )= y(xk )+ ∫ f (t, y(t))dt .

xk

С использованием аппроксимации (8.12) строится схема численного интегрирования дифференциального уравнения на этом сегменте:

xk +1

+ ∫F(xk , xk−1, xk −2 )(t − xk )(t − xk−1 )dt +

xk

xk +1

+ ∫F(xk , xk−1, xk−2 , xk−3 )(t − xk )(t − xk−1 )(t − xk−2 )dt .

xk

Полагая шаг интегрирования h постоянным для всей сеточной области Ωm , предыдущее выражение можно преобразовать к виду

179

yk+1 = yk + F(xk )h + F(xk , xk−1 )h2 2 + F(xk , xk−1, xk−2 )5h3 6 +

(8.13)

+ F(xk , xk−1, xk−2 , xk−3 )9h4 4 .

Погрешность аппроксимации уравнения (8.1) схемой (8.13) имеет четвертый порядок. Очевидный недостаток метода Адамса в том, что он не является самостартующим, поскольку для начала вычислений кроме начального условия y(0)= y0 необходимо дополнительно знать значения y1, y2 , y3 искомой функции. Недостающие значения могут быть определены методами Эйлера, Рунге–Кутты или другими способами. Кроме того, при интегрировании с переменным шагом h выражение (8.13) становится достаточно громоздким.

Вместе с тем метод Адамса имеет определенные преимущества, например, перед методом Рунге–Кутты. Действительно, для определения очередного значения искомой функции вычислять правую часть дифференциального уравнения приходится лишь один раз, тогда как в методе Рунге–Кутты четвертого порядка такие вычисления приходится выполнять четырежды. При достаточно сложном выражении, стоящем в правой части дифференциального уравнения, выигрыш во времени проведения расчетов может быть значительным.

В монографии [28] приведены схемы Адамса–Башфорта третьего yk+1 = yk + h[23 f (xk , yk ) −16 f (xk−1, yk−1 ) + 5 f (xk−2 , yk−2 )]12

и второго

yk+1 = yk + h[3 f (xk , yk ) − f (xk−1, yk −1 )]2

порядков погрешности аппроксимации.

8.7. Неявные схемы интегрирования

Рассмотренные методы решения задачи Коши относятся к явным схемам интегрирования, поскольку приводят на каждом шаге к решению одного линейного уравнения относительно очередного узлового значения искомой функции. Примером неявных схем являются разностные аналоги Адамса–Моултона [28] дифференциального уравнения (8.1) второго

yk+1 = yk + h[f (xk , yk ) + f (xk+1, yk+1 )]2

и четвертого

yk+1 = yk + h[9 f (xk +1, yk+1 ) +19 f (xk , yk )− 5 f (xk−1, yk−1 ) + f (xk−2 , yk−2 )]24

180