Метод Рунге-Кутта численного решения обыкновенного дифференциального уравнения.

В той же ситуации, что и выше, т.е. при поиске числа y(b) по числу y₀ = y(c) и равенству (4.8), существует еще один широко применяемый метод - метод Рунге-Кутта, который, как правило, быстрее приводит к числу y(b), чем метод Эйлера. Сформулируем действия по методу Рунге-Кутта:

1Й шаг. Фиксируем точность, с которой нужно найти значение .

Обозначим это число через ε. Поясним, что это означает, что числа, отличающиеся меньше, чем на ε, считаются одинаковыми.

2й шаг. Фиксируем произвольное n и разделим отрезок [a, b] на n равных частей: a = x₀ < x₁ < ... < x_n = b, где , x_i = a + i·h, i = 1,2, ..., n.

3й шаг. Построим последовательность чисел:

y_i = y_i-1 + λ_i-1, где

и

k₁ = h·f(x_i-1, y_i-1),

i = 1, 2, 3, ..., n.

в которой, напомним, y₀ = y(c). Обозначим y_n через U.

4й шаг. Заменим n на n+1и повторим шаги 2 и 3. Полученное число (т.е. последнее из вычисляемых на шаге 3) обозначим теперь через V.

5й шаг. Если окажется, что числа и отличаются друг от друга меньше, чем на , то число y(b) считается найденным и равным V. В противном случае переобозначим V через U и вернемся к шагу 4.

Можно доказать, что когда функция из (4.8) имеет непрерывные частные производные, описанный процесс обязательно конечен и ответ находится действительно с любой наперед заданной точностью.

#include <stdafx.h>

#include "rungekutta.h"

/*************************************************************************

Метод Рунге-Кутта четвертого порядка для решения уравнения первого порядка.

Решает дифференциальное уравнение y'=F(x,y) методом Рунге-Кутта 4 порядка.

Начальная точка имеет кординаты (x,y),

конечная - (x1, RungeKutt(x, x1, y, n)).

До конечной точки мы добираемся через n промежуточных

с постоянным шагом h=(x1-x)/n

*************************************************************************/

double solveoderungekutta(double x, double x1, double y, int n)

{

double result;

int i;

double h;

double y1;

double k1;

double k2;

double k3;

h = (x1-x)/n;

y1 = y;

i = 1;

do

{

k1 = h*f(x, y);

x = x+h/2;

y = y1+k1/2;

k2 = f(x, y)*h;

y = y1+k2/2;

k3 = f(x, y)*h;

x = x+h/2;

y = y1+k3;

y = y1+(k1+2*k2+2*k3+f(x, y)*h)/6;

y1 = y;

i = i+1;

}

while(i<=n);

result = y;

return result;

}

Численное решение системы обыкновенных дифференциальных уравнений.

Система обыкновенных дифференциальных уравнений - это система уравнений вида

f_i(x, y₁, ..., y_n, y₁', ..., y'_n, y₁^(S1), ..., y_n^(Sn)) = 0 (4.10)

Здесь x — аргумент, изменяющийся в некотором промежутке [a, b], а y₁, ..., y_n — функции этого аргумента, встречающиеся в равенствах системы (4.10) под знаками производных, порядок которых не выше S_i для y_i = 1, ..., n. Найти эти функции — пусть в виде таблиц значений, и означает решить систему (4.10). Существует класс задач, в которых надо найти совершенно определенные функции y₁, ..., y_n — функции, удовлетвряющие начальным условиям в некоторой точке :

y_i(c) = y_i⁽⁰⁾, y'_i(c) = y_i', y_i^(Si-1)(c) = y_iСуществует простой прием сведения произвольной системы уравнений типа (4.10) к системе уравнений, в которой производные не бывают порядка выше первого: если, допустим, S_i > 1, то полагают

y_i' = z_i, y_i'' = z_i', ..., y_i^(Si-1)

^(Si-1)(c) = z_i^(Si-2),

после чего в системе (4.10) появляется новая неизвестная функция z_i, но высший порядок производной, под которым встречается в новом варианте функция уже на единицу меньше, а высший порядок производной новой функции z_i, с которым она входит в новую систему, тоже не больше S_i – 1. Легко заметить, что таким приемом можно любую систему (4.10) свести к системе первого порядка:

f_i(x, y₁, ..., y_n, y₁', ..., y'_n) = 0 (4.11),

Так же, как и в случае одного уравнения, принято анализировать лишь ту ситуацию, при которой (4.11) можно разрешить относительно производных, т.е. свести всё к системе

y_i' = g_i(x, y₁, ..., y_n,) (4.12)

y_i(c) = y_i⁰, i = 1, 2, ..., n.

Легко, однако, для системы (4.12) воспроизвести, как формулы Эйлера, так и формулы Рунге-Кутта, причем утверждение о том, что с помощью этих формул можно приблизиться как угодно точно к значениям неизвестных функций в некоторой точке x = b > a при наличии у функций g_i определенных свойств, остается справедливым.

#include <stdafx.h>

#include "rungekuttasys.h"

/*************************************************************************

Метод Рунге-Кутта четвертого порядка для решения

системы дифферециальных уравнений.

Алгоритм решает систему дифференциальных уравнений y[i]'=F(i,x,y) для i=1..n

методом Рунге-Кутта 4 порядка.

Начальная точка имеет кординаты (x, y[1], ..., y[n])

До конечной точки мы добираемся через n промежуточных

с постоянным шагом h = (x1-x)/m

Результат помещается в переменную y

*************************************************************************/

void solvesystemrungekutta(const double& x,

const double& x1,

const int& m,

const int& n,

double& y[])

{

double h;

int i;

for(i = 0; i <= m-1; i++)

{

solvesystemrungekuttastep(x+i*(x1-x)/m, (x1-x)/m, n, y);

}

}

/*************************************************************************

Один шаг метода Рунге-Кутта четвертого порядка для решения

системы дифферециальных уравнений.

Алгоритм совершает один шаг метода для системы

диффуров y[i]'=F(i,x,y) для i=1..n

Начальная точка имеет кординаты (x,y[1], ..., y[n])

После выполнения алгоритма в переменной y содержится состояние

системы в точке x+h

*************************************************************************/

void solvesystemrungekuttastep(const double& x,

const double& h,

const int& n,

double& y)

{

int i;

double yt[n];

double k1[n];

double k2[n];

double k3[n];

double k4[n];

for(i = 1; i <= n; i++)

{

k1[i] = h*f(i, x, y);

}

for(i = 1; i <= n; i++)

{

yt[i] = y[i]+0.5*k1[i];

}

for(i = 1; i <= n; i++)

{

k2[i] = h*f(i, x+h*0.5, yt);

}

for(i = 1; i <= n; i++)

{

yt[i] = y[i]+0.5*k2[i];

}

for(i = 1; i <= n; i++)

{

k3[i] = h*f(i, x+h*0.5, yt);

}

for(i = 1; i <= n; i++)

{

yt[i] = y[i]+k3[i];

}

for(i = 1; i <= n; i++)

{

k4[i] = h*f(i, x+h, yt);

}

for(i = 1; i <= n; i++)

{

y[i] = y[i]+(k1[i]+2.0*k2[i]+2.0*k3[i]+k4[i])/6;

}

}

Метод Булирша-Штера с использованием рациональной экстраполяции для системы уравнений

Метод Булирша-Штера (Bulirsch-Stoer Method) - это метод решения системы обыкновенных дифференциальных уравнений первого порядка с гладкими правыми частями. Гладкость правых частей является необходимой для работы метода. Если правые части системы не являются гладкими или содержат разрывы, то лучше использовать метод Рунге-Кутта. В случае же гладкой системы метод Булирша-Штера позволяет добиться существенно большей точности, чем метод Рунге-Кутта.

Принцип работы метода

Основной идеей метода является вычисление состояния системы в точке x+h, как результата двух шагов длины h/2, четырех шагов длины h/4, восьми шагов длины h/8 и так далее с последующей экстраполяцией результатов. Метод строит рациональную интерполирующую функцию, которая в точке h/2 проходит через состояние системы после двух таких шагов, в точке h/4 проходит через состояние системы после четырех таких шагов, и т.д., а затем вычисляет значение этой функции в точке h = 0, проводя экстраполяцию.

Гладкость правых частей приводит к тому, что вычисленное при помощи экстраполяции состояние системы оказывается очень близко к действительному, а использование рациональной эстраполяции вместо полиномиальной позволяет ещё больше повысить точность.

Таким образом проводится один шаг метода, после чего принимается решение — следует ли изменять шаг, а если да — то в какую сторону. При этом используется оценка погрешности, которую мы получаем в качестве дополнительного результата при рациональной экстраполяции.

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

Описание программы

Для работы программы требуется объявить функцию F, принимающую массив (нумерация элементов от 1 до n), хранящий состояние системы, и заполняющую второй переданный ей массив (нумерация элементов от 1 до n) производными от переменных по времени.

#include <stdafx.h>

#include "bulirschstoer.h"

/*************************************************************************

Решение автономной системы диффуров первого порядка методом Булирша-Штера

(метод рациональной экстраполяции с переменным шагом).

Входные параметры:

x1 - начало временного отрезка интегрирования [x1, x2]

x2 - конец временного отрезка интегрирования [x1, x2]

y - начальное состояние системы.

Массив с нумерацией элементов от 1 до n.

n - число неизвестных

h - начальное значение шага интегрирования

err - допустимая погрешность интегрирования на каждом шаге

IsAbsErr- тип погрешности err - абсолютная или относительная

Выходные параметры:

y - конечное состояние системы (в точке t=xfin).

Массив с нумерацией элементов от 1 до n.

Процедура использует определенную программистом процедуру F, вычисляющую

по состоянию системы производные от переменных.

Входные параметры процедуры F:

Y - состояние системы. Массив с нумерацией элементов от 1 до n.

Выходные параметры процедуры F:

dY - производные от переменных системы по времени.

Массив с нумерацией элементов от 1 до n.

*************************************************************************/

void solvesystembulirschstoer(double x1,

double x2,

doudle& y[],

const int n,

double h,

double err,

char isabserr)

{

double ya[n];

double yl;

double ym;

double dy;

double dz;

double s;

ap::real_2d_array dt;

double d;

ap::real_2d_array yg;

ap::real_2d_array yh;

bool konv;

bool bo;

bool bh;

bool fin;

int i;

int j;

int k;

int kk;

int l;

int m;

int jj;

int r;

int sr;

double ta;

double xxfin;

double a;

double b;

double c;

double b1;

double fc;

double g;

double ho;

double u;

double v;

int cycleflag;

xxfin = x2;

cycleflag = -1;

while(true)

{

if( cycleflag==-1||cycleflag==22 )

{

f(y, dz);

bh = false;

fin = false;

for(i = 1; i <= n; i++)

{

s(i) = 0;

ya(i) = y(i);

for(k = 1; k <= 7; k++)

{

dt(i,k) = 0;

}

}

ho = h;

if( x1+h>=xxfin )

{

ho = x2-x1;

fin = true;

}

}

if( cycleflag==-1||cycleflag==20||cycleflag==22 )

{

a = ho+x1;

if( a==x1 )

{

a = a;

return;

}

fc = 1.5;

bo = false;

m = 1;

r = 2;

sr = 3;

jj = 0;

}

cycleflag = -1;

for(j = 1; j <= 10; j++)

{

if( bo )

{

d(2) = 1.7777777777778;

d(4) = 7.1111111111111;

d(6) = 28.444444444444;

}

else

{

d(2) = 2.25;

d(4) = 9.0;

d(6) = 36.0;

}

konv = false;

if( j>3 )

{

konv = true;

}

if( j>7 )

{

l = 7;

d(7) = 64.0;

fc = 0.6*fc;

}

else

{

l = j;

d(l) = m*m;

}

m = m*2;

g = ho/m;

b = g*2.0;

if( bh&&j<9 )

{

for(i = 1; i <= n; i++)

{

ym(i) = yh(j,i);

yl(i) = yg(j,i);

}

}

else

{

kk = (m-2)/2;

m = m-1;

for(i = 1; i <= n; i++)

{

yl(i) = ya(i);

ym(i) = ya(i)+g*dz(i);

}

for(k = 1; k <= m; k++)

{

f(ym, dy);

for(i = 1; i <= n; i++)

{

u = yl(i)+b*dy(i);

yl(i) = ym(i);

ym(i) = u;

u = fabs(u);

if( u>s(i) )

{

s(i) = u;

}

}

if( k==kk&&k!=2 )

{

jj = 1+jj;

for(i = 1; i <= n; i++)

{

yh(jj,i) = ym(i);

yg(jj,i) = yl(i);

}

}

}

}

f(ym, dy);

for(i = 1; i <= n; i++)

{

v = dt(i,1);

dt(i,1) = (ym(i)+yl(i)+g*dy(i))*0.5;

c = dt(i,1);

ta = c;

if( l>=2 )

{

for(k = 2; k <= l; k++)

{

b1 = d(k)*v;

b = b1-c;

u = v;

if( b!=0 )

{

b = (c-v)/b;

u = c*b;

c = b1*b;

}

v = dt(i,k);

dt(i,k) = u;

ta = u+ta;

}

}

if( !isabserr&&fabs(y(i)-ta)>err*s(i) )

{

konv = false;

}

if( isabserr&&fabs(y(i)-ta)>err )

{

konv = false;

}

y(i) = ta;

}

if( konv )

{

x1 = a;

h = fc*h;

if( fin )

{

return;

}

cycleflag = 22;

break;

}

d(3) = 4.0;

d(5) = 16.0;

bo = !bo;

m = r;

r = sr;

sr = m*2;

}

if( cycleflag==-1 )

{

bh = !bh;

ho = ho*0.5;

h = ho;

fin = false;

cycleflag = 20;

}

}

}