4.3 Постановка задачи о численном решении обыкновенного дифференциального уравнения.

Обыкновенным дифференциальным уравнением называется равенство

f(x, y, y', y'', ..., y⁽ⁿ⁾) = 0 (4.8),

в котором x — независимая переменная, изменяющаяся в некотором отрезке [a, b], а y — неизвестная функция от x, которую и надо найти. Различают два типа обыкновенных дифференциальных уравнений — уравнения без начальных условий и уравнения с начальными условиями. Уравнения без начальных условий — это как раз то, что было только что определено. А уравнение с начальными условиями — это записанное выше уравнение относительно функции y, но в котором требуется найти лишь такую функцию y, которая удовлетворяет при некоторомследующим условиям:

y(c) = y₀, y'(c) = y'₀, y''(c) = y''₀, ..., y^(n-1)(c) = y^(n-1)₀, т.е. в точке c функция y и ее первые n-1 производных принимают наперед заданные значения. В этой ситуации число n называется порядком уравнения.

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

Предположим, что порядок уравнения (4.8) равен 1 и, более того, его

можно представить виде:

y' = f(x,y) (4.9)

В случае первого порядка начальные условия превращаются в единственное равенство: y(c) = y₀.

Искомая функция при численном подходе к проблеме решения всегда воспринимается через какую-то таблицу своих значений; иными словами, по заданным значениям аргумента надо найти такие y = y₀, y₁, ..., y_n, чтобы таблица

x	x1	x2	x3	...	xn
f(x)	y1	y2	y3	...	yn

оказалась таблицей значений искомой функции. Для этого, в свою очередь, достаточно по заданному значению уметь находить значение в любой точке. Вот метод Эйлера для этого последнего поиска:

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 + h·f(x_i-1, y_i-1), i = 1, 2, ..., n, в которой, напомним, y₀ = y(c). Обозначим y_n через U.

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

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

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

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

В той же ситуации, что и выше, т.е. при поиске числа 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"

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

double f(int i, double x, double& y)

{

/*Функция которую будем решать*/

return;

}

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

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

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

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

Алгоритм решает систему дифференциальных уравнений 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;

}

}