Деякі скінченно-різнецеві методи розвязування звичайних диференціальних рівнянь (стр. 3 из 6)
- Перетворюємо вектор Е на одне число є, що характеризує погрішність на цьому кроці.
Ця схема призводить до досить забавної поведінки "швидких" змінних. Поведінка виявляється абсолютно однаковою як у разі релаксаційного рівняння
, так і у разі осциляторного рівняння 
. Саме, при Ch>>1 "швидкі" змінні експоненциально вибуваєть, причому швидкість вибування абсолютно не залежить від величини C: x(t+h)=x(t)/3. Отже стійкість гарантована, але сподіватися на правильний опис еволюції "швидких" змінних не доводиться.Повний алгоритм рішення системи ОДУ з адаптивною зміною кроку будується на основі цього рецепту таким самим чином, як будувався повний алгоритм в звичайному (явному) методі Рунге-Кутта п'ятого по-рядка. Єдина модифікація полягає в тому, що міри "1/4" і "1 /5", відповідна п'ятому порядку точність явної схеми, належить замінити на мірі "1/3" і "1/4", відповідна четвертому порядку точність неявної схеми.
Необхідно також нагадати, що для "жорсткої" системи ОДР перетворити вектор Е в одне число 6, що оцінює погрішність, набагато важче, чим для звичайної системи ОДР.
6. Неявні інтерполяційні схеми
Алгоритм неявних інтерполяційних методів дослівно повторює алгоритм звичайних (явних) інтерполяційних методів. Єдина відмінність полягає в тому, що ми не використовуємо звичайну (явну) схему Рунге-Кутта другого порядку:
Ми заміняємо її на неявну схему Рунге-Кутта другого порядку:
Нагадаємо, що для цієї схеми можна використовувати два різних рецепта:
• шукати точне рішення цієї системи нелінійних рівнянь методом Ньютона;
• обмежитися лінеаризованим рівнянням, тобто зробити тільки перший крок методу Ньютона :
У будь-якому випадку доведеться вирішувати СЛР. В більшості випадків можна обійтись другим варіантом, але в дуже нестійких завданнях доведеться вдатися до першого (до речі, при цьому не доведеться сильно міняти текст програми).
Нагадаємо, що існують варіанти інтерполяційного методу з фіксованим кроком Н і змінним порядком m, з фіксованому порядку те і змінним кроком Н і, нарешті, із змінними m і Н. Тут можна привечти тільки простий алгоритм з фіксованим H.
Алгоритм:
Вибираємо постійний крок H за часом (як правило слід брати H близько 0.2). Рухаємося по траєкторії з цим кроком. Кожен окремий крок реалізується так:
Для k = 1,2,.. виконуємо наступне:
Для чергового Nk з набору {4, 6, 8, 12, 16, 24, 32, 48, 64, 96, 128} обчислюємо крок hk = H/Nk.
За допомогою цього кроку hk знаходимо
тут tn = t+nhk, n = 0.. Nk-1; при цьому
. Зрозуміло, тут не слід обчислювати зворотну матрицю, потрібно просто вирішувати СЛР.Для кожного хіпо k точок:
виконуємо поліноміальну інтерполяцію взалежності Y(X) в точку X = 0. Результат інтерполяції, отриманий на k -му етапі, означаємо
. Оцінюємо погрішність 
як 
Перетворюємо вектор
на число 
характеризуючи погрішність. Якщо 
менше замовленого 
те вважаємо 
і завершуємо даний крок t→t+H.Перша неприємна відмінність від явного методу полягає в тому, що кожен маленький крок на hk (великий крок H складається з Nk штук таких кроків) вимагає рішення СЛР.
Друга неприємна відмінність від явного методу полягає в тому, що тут набагато складніше перетворити вектор
на одне число 
, оцінюючи погрішність. На жаль, рецепт цього перетворення слідує підібрати залежно від конкретного завдання, тобто від конкретних властивостей "швидких" і "повільних" змінних.На перший погляд, цей алгоритм приведе не просто до багатьох помилок в описі динаміки "швидких" змінних, а до помилок якісних. Дійсно, якщо рівняння для "швидкої" змінної має вигляд
, то при Ch>>1 неявна схема Рунге-Кутта другого порядку замість експоненціального вибуває дає x(t+h)≈-x(t), тобто |x(t+H)|≈|x(t)|. На щастя, наступна інтерполяція виправляє ситуацію. Інтерполяційний алгоритм "помічає", що у міру зменшення hk абсолютна величина x(k)(t+H) зменшується. Причому це зменшення тим помітніше, чим менше hk - В результаті остаточна відповідь 
буде помітна менше за абсолютною величиною, чим x(t).7. Програма Рунге-Кутта на мові С#
Наведемо приклад пограми Рунге-Кутта на мові С#. В програмі використовується абстрактний клас TrungeKutta, в якому потрібно перекрити абстрактний метод F, який задає перші чаcтини рівняння.
using System;
using System.Collections.Generic;
namespace rwsh_rk4
{
abstract class TRungeKutta
{
public int N;
double t; // теперішній час
public double[] Y; // шукане число Y[0] – саме рішення, Y[i] - i-та змінна розвязку
double[] YY, Y1, Y2, Y3, Y4; // внутрішня змінна
public TRungeKutta(int N) // N – розмір системи
{
this.N = N; // зберегти розміри системи
if (N < 1)
{
this.N = -1; // якщо розмірність менше одиниці, то установити флаг помилки
return; // і вийти із конструктора
}
Y = new double[N]; // створення вектора розв’язку
YY = new double[N]; // внутрішній розв’язок
Y1 = new double[N];
Y2 = new double[N];
Y3 = new double[N];
Y4 = new double[N];
}
public void SetInit(double t0, double[] Y0) // встановлення початкових умов.
{ // t0 – початковий час, Y0 – початкова умова
t = t0;
int i;
for (i = 0; i < N; i++)
{
Y[i] = Y0[i];
}
}
public double GetCurrent() // повернути даний час
{
return t;
}
abstract public void F(double t, double[] Y, ref double[] FY); // перші частини с-ми.
public void NextStep(double dt) // наступний крок метода Рунге-Кутта, dt – крок по часу
{
if(dt<0)
{
return;
}
int i;
F(t, Y, ref Y1); // вирахувати Y1
for (i = 0; i < N; i++)
{
YY[i] = Y[i] + Y1[i] * (dt / 2.0);
}
F(t + dt / 2.0, YY, ref Y2);
for (i = 0; i < N; i++)
{
YY[i] = Y[i] + Y2[i] * (dt / 2.0);
}
F(t + dt / 2.0, YY, ref Y3); // вирахувати Y3
for (i = 0; i < N; i++)
{
YY[i] = Y[i] + Y3[i] * dt;
}
F(t + dt, YY, ref Y4); // вирахувати Y4
for (i = 0; i < N; i++)
{
Y[i] = Y[i] + dt / 6.0 * (Y1[i] + 2.0 * Y2[i] + 2.0 * Y3[i] + Y4[i]); // виразувати р-зок на новому кроці
}
t = t + dt; // збільшити крок
}
}
class TMyRK : TRungeKutta
{
public TMyRK(int aN) : base(aN) { }
public override void F(double t, double[] Y, ref double[] FY)
{
FY[0] = Y[1]; // приклад математичний маятник
FY[1] = -Y[0]; // y''(t)+y(t)=0
}
}
class Program
{
static void Main(string[] args)
{
TMyRK RK4 = new TMyRK(2);
double[] Y0 = {0, 1}; // задаємо початкові умови y(0)=0, y'(0)=1
RK4.SetInit(0, Y0);
while (RK4.GetCurrent() < 10) // розв’язуєм до 10
{
Console.WriteLine("{0}\t{1}\t{2}", RK4.GetCurrent(), RK4.Y[0], RK4.Y[1]); // вивести t, y, y'
RK4.NextStep(0.01); // розв’язати на наступному кроці, крок інтегрування dt=0.01
}
}
}
}
8. Програма Beeman
У програмі Вееman моделюється осцилятор Морза з допомогою алгоритму Бімана. Оскільки цей алгоритм не самостартуючий, то для всіх - числових значень х(, і використовується швидкісна форма алгоритму Верле.
PROGRAM Beeman І моделювання осцилятора Морза
CALL initialfx, v, aold, dt, dt2, nmax)
CALL energy(x, v, ecum, e2cum) 1 значення початкової енергії
CALL Verleg x, v, a, aold, dt, dt2) CALL energy{ x, v, ecum, e2cum) LET n = 1
DO whiie n < nmax
LET n = n + 1 1 число кроку
CALL Bceman{x, v, a, aold, dl, dl2)
І образування повної знергії після кожного кроку за часом CALL energY(x, v, ecum, e2cum) LOOP
CALL output{ ecum, e2cum, n) END
SUB initial( x, v, aold, dt, dt2, nmax) DECLARE DEF f LET х = 2 LET v = 0 LET aold = f(x)
INPUT prompt "крок по часу (c) = ": dt LET dt2 = dt" dt
INPUT prompt "тривалість = ": tmax LET nmax = tmax/dt END SUB
SUB Ver!et( x, v, a, aold, dt, dt2) DECLARE DEF f
LET x = x + v*dt + 0.5*ao!d*dt2 LET a = f(x)
LET v = v + 0.5"{a + ao!d)*dt END SUB
SUB Beeman(x, v, a, aold, dt, dt2) DECLARE DEF f
LET x = x + v*dt + (4*a - aold)*dt2/6
LET anew = f(x) І значення на (п+1) -му кроці LET v = v + (2*anew + 5*a - aold)*dt/6
LET aold = a значення на (n-1) -му кроці
LET а = anew значення на n-му кроці END SUB
DEF f(x) LET e = exp(- x) LET f = 2*e*(e - 1) END DEF
SUB energY(x, v, ecum, e2cuin) LET KE = 0.5*v" v LET e = exp(- x) LET PE = e*{e - 2) LET etot = KE + PE LET ecum = ecum + etot LET e2cum = e2cum + etot*etot END SUB
SUB output{ecum, e2cuiT!, n) LET n = n + 1 І вирахування початкового значення
LET ebar = ecum/n PRINT "середня енергія = ";ebar LET sigma2 = e2cum/n - ebar*ebar PRINT "sigma = "; sqr(sigma2) END SUB
Метод Адамса
Цей метод чисельного інтегрування розроблений Адамсом в 1855 році на прохання видомого англійського алтелериста Башфора, який займався внутрішньою балістикою. В подальшому цей метод був забутий і знову відкритий був норвезьким математиком Штермером. Популяризація метода Адамса і подальше його вдосконалення пов’язане із іменем Крилова.