logo
Деякі скінченно-різнецеві методи розв’язування звичайних диференціальних рівнянь

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} {1} {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 році на прохання видомого англійського алтелериста Башфора, який займався внутрішньою балістикою. В подальшому цей метод був забутий і знову відкритий був норвезьким математиком Штермером. Популяризація метода Адамса і подальше його вдосконалення повязане із іменем Крилова.

Запишемо рівняння першого порядку

З початковими умовами (1,2)

Нехай xi(i=0,1,2…)-система рівнозначних значень з кроком h i y(xi). Очевидно маємо

(3)

В силу другої інтерполяційної формули Ньютона з точністб до різниць четвертого порядку отримуємо:

(4)

де або (4а)

Підставляю вираз (4а) в формулу (3) і враховуючи те, що будемо мати

З відси отримуємо формулу експоляриціональну Адамса

(5)

Для початкового процессу потрібно чотири початкових значення y0, y1, y2, y3, - початковий відрізок, який приділяє, виходячи із початкових умов (2), яким-небуть чисельним методом. Мажна наприклад використати метод Рунге-Кутта або розкласти в ряд Тейлора

Де i=1,2,3 (або i=-1,1,2) із відповідною зміною нумерування. Знаючи ці значення, із рівнянь (1) можна знайти значення похідних і скласти таблицю

(6)

Подальше значення yi (i=4,5…) шуканого розвязку можна крок за кроком обчислювати за формулою Адамса, поповнюючи по мірі можливості таблицю різниць (6)

Вирахувавши перше наближення для по формулі

Визначити підрахувати кінцеві різниці

(7)

а потім знайти друге наближення для більш точній формулі

(8)

Якщо і відрізняються лишень на дкілька одиниць останнього зберігаючого десяткового розряду, то можна поставити а потім знайшовши перерахувавши кінцеві різниці (7). Після цього, потрібно знову знайти по формулі (8) Поту цей крок h повинен бути таким, щоб цей перерахунок був зміненим.

На практиці крок h вибирають малим, щоб можна було знехтувати членом в формулі (8)

Якщо за розбіжність величин і суттєва, то потрібно зменшити крок h.

Звичайно крок h зменшують рівно в 2 рази. Можна показати, як в цьому випадку, маючи до деякого значення і таблицю величин хj, yj, Yj=hyj (j<=i) з кроком , можна просто побудувати таблицю величин з кроком

На основі формули (4) будемо мати

(9)

Де Звідси, і і враховуючи, що заходимо

(10)

Аналогічно при із формули (9) отримаєм, що аргументу відповідає значення

(11)

Що стосується значень Yi-1 i Yi, то вони знаходяться в старій таблиці. Після цього складаємо початковий відрізок для нової таблиці:

і знаходимо кінцеві різниці:

Далі таблиця будується простим способом, подальшою модифікацією формули (5):

Для роботи на компютерах формулу Адамса (5) вигідно використовувати в розкритому виді. Враховуючи, що

Після цього маємо: причому

Метод Крилова

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

(1)

З початковими умовами

Введемо спочатку ряд допоміжних формул

В силу формули Адамса отримаємо

(2)

Введемо позначення

Формула (2) називається формулою похилого рядка, так як в ній використовуються різниці, які знаходяться на діагоналі таблиці різниць. Враховуючи, що

Із формули (2) будемо мати

Звідси отримуємо першу допоміжну формулу - яку ще можна назвати перша формула ламаного рядка

(3)

Далі враховуючи, що і із формули (3) виводимо другу формулу - друга формула ламаного рядка

(4)

Якщо отримаємо формулу горизонтального рядка

(5)

Підмітимо, що формулу (5) можна отримати безпосередньо за допомогою інтегрування, в межах від xi до xi+1 розкладанням за допомогою першої інтерполяційної формули Ньютона:

Перейдемо до опису метода Крилова послідовних наближень. Перше наближення полягає у тому, щоб знайти наближене значення

Після цього знайдемо і складає різницю , де .

Значення які знайшли заносимо в розділ (І) основного бланку (таблиця 1)

Схема обчислення відрізка методом послідовних наближень

№ наближення

і

x

y

І

0

1

x0

x1

ІІ

0

1

2

x0

x1

x2

ІІІ

0

1

2

3

x0

x1

x2

x3

Далі переходимо до другого наближення. Для того, використовуємо дані із знаходження ламаних рядків, обчислюємо значення і :

(7)

(8)

Двочленні формули отримуються відповідно із формули (5) при і=0 і із формули (2) при і=1 в результаті відкидання різниць порядка вищого ніж перший.

Таким чином, отримаємо можливість знайти

і ,

в результаті чого можна порахувати

і скласти різниці

Отримані результати записуємо у таблицю в розділ 2 основного бланка

Для знаходження третього наближення застосовуємо трьохчленні формули, які отримуються із формули (2) при і=2 після відкидання різниць третього порядку. Обчислюємо значення із трьохчленних формул:

(9)

(10)

(11)

Звідси можна знайти

і обчислити . Після цього можна заповнити розділ ІІІ в таблиці (І) знайшовши потрібні різниці звичайним порядком.