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)
Звідси можна знайти
і обчислити . Після цього можна заповнити розділ ІІІ в таблиці (І) знайшовши потрібні різниці звичайним порядком.
- Розв’язування звичайних диференціальних рівнянь в системі mathcad
- Тема № 13. Елементи теорії звичайних диференціальних рівнянь
- 7.9. Чисельні методи розв’язування задачі Коші для систем звичайних диференціальних рівнянь першого порядку
- 7.10. Чисельні методи розв’язування задачі Коші для звичайних диференціальних рівнянь вищих порядків
- 7.8. Засоби середовища matlab розв’язування задачі Коші для звичайних диференціальних рівнянь першого порядку з використанням багатокрокових методів
- Лекція 7. Методи рішення звичайних диференціальних рівнянь
- Тема № 13. Елементи теорії звичайних диференціальних рівнянь
- 1.2 Методи розв’язування рчп