Система дифференциальных уравнений методом рунге кутта 4 порядка

Метод Рунге-Кутта решения диф. уравнений и их систем.

Метод позволяет решать системы обыкновенных дифференциальных уравнений (ОДУ) первого порядка следующего вида:

которые имеют решение:

где t — независимая переменная (например, время); X, Y и т.д. — искомые функции (зависимые от t переменные). Функции f, g и т.д. — заданы. Также предполагаются заданными и начальные условия, т.е. значения искомых функций в начальный момент.

Одно диф. уравнение — частный случай системы с одним элементом. Поэтому, далее речь пойдет для определенности о системе уравнений.

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

Метод Рунге-Кутта заключается в рекурентном применении следующих формул:

Реализация Метода Рунге-Кутта на Delphi может выглядеть так (привожу полностью модуль):

Модуль полностью работоспособен. Возвращаемое функцией Runge_Kutt значение — код ошибки. Вы можете дополнить список ошибок по своему усмотрению. Рассчитанные функции системы помещаются в массив Res. Чтобы не загромождать код, в модуле опущены проверки (типа блоков try). Рекомендую их добавить по своему усмотрению.

Ниже приводится описание функции Runge_Kutt и типов, использующихся в модуле.

• FunArray — вектор функций (правых частей уравнений системы);
• First, Last — начальная и конечная точки расчетного интервала;
• Steps — число шагов по расчетному интервалу;
• InitArray — вектор начальных значений
• var Res — матрица результатов включая независимую переменную.

В модуле описаны типы:

Функция возвращает коды ошибок:

• 0 — нет ошибок;
• 100 — число уравнений не равно числу начальных условий.

Решение содержится в переменной-матрице Res. Первый индекс матрицы относится к переменной (0 — независимая переменная, 1 — первая зависимая и т.д.), второй — к номеру расчетной точки (0 — начальная точка).

Рассмотрим один пример использования модуля. Создадим новое приложение и подключим к нему модуль. На форме приложения разместим кнопку Button1 и область текста Memo1. Поместим в приложение две функции и обработчик нажатия кнопки:

Нажатие кнопки приведет к расчету точек системы, которые будут выведены в текстовую область.

Модуль с примером и справкой можно скачать бесплатно по адресу RK.zip (ZIP, 15,3Kb) (русский вариант). Английский вариант (условно-бесплатный) можно скачать по адресу RK_Eng.zip (ZIP, 23.4Kb)

Ссылки

• http://sadovoya.narod.ru/RK.zip (русский вариант).
• http://sintreseng.narod.ru/RK_Eng.zip (английский, условно-бесплатный вариант)

Оставить комментарий

Комментарии

Скачала по Вашей ссылке русский вариант, изменила для своей системы диф. уравнений, но при запуске выдаёт ошибку :
Project Ex.exe raised exception class EOverflow with message ‘ Floating point overflow ‘
Помогите, пожалуйста .

Вот изменённый мною модуль:

unit Unit1;
interface
uses
SysUtils, Forms, StdCtrls, Controls, Classes, Dialogs, Math;
type
TForm1 = class(TForm)
Memo1: TMemo;
rk_But: TButton;
procedure rk_ButClick(Sender: TObject);
private
< Private declarations >
public
< Public declarations >
end;
var
Form1: TForm1;
pn,k,ro,Pzv: Extended;

implementation
uses rk_method, Windows;

<$R *.dfm>
procedure Syst (var t: TFloat; var X: TFloatVector;
var RP: TFloatVector);
const
fdr1=0.503;
fdr2=0.503;
fdr3=0.196;
W1=179.8928;
W2=3773.8568;
W3=2504.1203;
b1=55.9203;
b2=98.6;
b3=98.6;
Ls1=3.78;
Ls2=9;
Ls3=15.3;
Svidj2=1352.438;
Svidj3=1352.438;
my=0.62;
vk=30;
m=1.2;
L1=30.969;
L2=42.131;
delta1=0;

begin
pn:=2.5*Power(10,4);
k:=6*Power(10,-7);
ro:=8.5*Power(10,-7);
Pzv:=3.919*Power(10,7);

RP[0] := (1/(k*W1))*(my*fdr1*sqrt(2/ro)*sqrt(Abs(pn-X[0]))-my*fdr2*sqrt(2/ro)*sqrt(Abs(X[0]-X[1]))-(delta1*delta1*delta1*b1)/(12*ro*vk*Ls1)*X[0]); // dp1/dt
RP[1] := (1/(k*W2))*(my*fdr2*sqrt(2/ro)*sqrt(Abs(X[0]-X[1]))-my*fdr3*sqrt(2/ro)*sqrt(Abs(X[1]-X[2]))-(X[4]*X[4]*X[4]*b2)/(12*ro*vk*Ls2)*X[1]); // dp2/dt
RP[2] := (1/(k*W3))*(my*fdr3*sqrt(2/ro)*sqrt(Abs(X[1]-X[2]))-(X[6]*X[6]*X[6]*b3)/(12*ro*vk*Ls3)*X[2]); // dp3/dt;
RP[3] := (((Svidj2*X[1]*(L1+L2))/L1)-Pzv)*(2/m); // dv2/dt
RP[4] := X[3]; // d delta2/dt
RP[5] := (((Svidj3*X[2]*(L1+L2))/L2)-Pzv)*(2/m); // dv3/dt
RP[6] := X[5]; // d delta3/dt
end;

procedure TForm1.rk_ButClick(Sender: TObject);
var
I, t1, t2: Cardinal;
tOut, InitConds: TFloatVector;
XOuts: TFloatMatrix;
Points: Cardinal;
First, Last: TFloat;
StepsFact: Cardinal;
Count: Word;
begin
Memo1.Clear;
First := 0.0;
Last := 10.0;
Count:= 7;
Points:=10+1; //11 points for output
StepsFact:=1000000; //all steps inside function = 10*StepsFact

try
SetLength(InitConds, Count);
InitConds[0]:=0.0; //x0(0)=0
InitConds[1]:=0.0; //x1(0)=0
InitConds[2]:=0.0; //x2(0)=0
InitConds[3]:=0.0; //x3(0)=0
InitConds[4]:=0.0; //x4(0)=0
InitConds[5]:=0.0; //x5(0)=0
InitConds[6]:=0.0; //x6(0)=0

SetLength(tOut, Points);
SetLength(XOuts, Count, Points);
except
ShowMessage(‘Out of memory. ‘);
exit;
end;

Метод Рунге-Кутты 4-го порядка для решения дифференциального уравнения

Учитывая следующие входные данные,

• Обычное дифференциальное уравнение, которое определяет значение dy / dx в виде x и y.
• Начальное значение y, т. Е. Y (0)

Таким образом, мы приведены ниже.

Задача состоит в том, чтобы найти значение неизвестной функции y в заданной точке x.

Метод Рунге-Кутты находит приблизительное значение y для данного x. Только обыкновенные дифференциальные уравнения первого порядка могут быть решены с помощью метода 4-го порядка Рунге Кутты.

Ниже приведена формула, используемая для вычисления следующего значения y _{n + 1} из предыдущего значения y _n . Значения n равны 0, 1, 2, 3,… (x — x0) / h. Здесь h — высота шага, а x _{n + 1} = x ₀ + h

, Меньший размер шага означает большую точность.

Формула в основном вычисляет следующее значение y _{n + 1,} используя текущее значение y _n плюс средневзвешенное значение четырех приращений.

• k ₁ — приращение, основанное на наклоне в начале интервала, используя y
• K ₂ представляет собой приращение на основе наклона в средней точке интервала, с использованием Y + ₁ HK / 2.
• к ₃ снова приращение на основе наклона в средней точке, используя при помощи у + кк _2/2.
• k ₄ — это приращение, основанное на наклоне в конце интервала, с использованием y + hk ₃ .

Этот метод является методом четвертого порядка, это означает, что локальная ошибка усечения имеет порядок O (h 5 ), в то время как общая накопленная ошибка составляет порядок O (h 4 ).

Ниже приведена реализация приведенной выше формулы.

// C программа для реализации метода Рунге Кутты
#include

// Пример дифференциального уравнения «dy / dx = (x — y) / 2»

float dydx( float x, float y)

// Находит значение y для заданного x, используя размер шага h
// и начальное значение y0 в x0.

float rungeKutta( float x0, float y0, float x, float h)

// Подсчитать количество итераций, используя размер шага или

int n = ( int )((x — x0) / h);

float k1, k2, k3, k4, k5;

// Итерация по количеству итераций

// Применить формулы Рунге Кутты, чтобы найти

// следующее значение у

k2 = h*dydx(x0 + 0.5*h, y + 0.5*k1);

k3 = h*dydx(x0 + 0.5*h, y + 0.5*k2);

k4 = h*dydx(x0 + h, y + k3);

// Обновить следующее значение y

y = y + (1.0/6.0)*(k1 + 2*k2 + 2*k3 + k4);;

// Обновляем следующее значение x

float x0 = 0, y = 1, x = 2, h = 0.2;

printf ( «\nThe value of y at x is : %f» ,

rungeKutta(x0, y, x, h));

// Java-программа для реализации метода Рунге Кутты

double dydx( double x, double y)

// Находит значение y для заданного x, используя размер шага h

// и начальное значение y0 в x0.

double rungeKutta( double x0, double y0, double x, double h)

differential d1 = new differential();

// Подсчитать количество итераций, используя размер шага или

int n = ( int )((x — x0) / h);

double k1, k2, k3, k4, k5;

// Итерация по количеству итераций

for ( int i = 1 ; i

// Применить формулы Рунге Кутты, чтобы найти

// следующее значение у

k1 = h * (d1.dydx(x0, y));

k2 = h * (d1.dydx(x0 + 0.5 * h, y + 0.5 * k1));

k3 = h * (d1.dydx(x0 + 0.5 * h, y + 0.5 * k2));

k4 = h * (d1.dydx(x0 + h, y + k3));

// Обновить следующее значение y

y = y + ( 1.0 / 6.0 ) * (k1 + 2 * k2 + 2 * k3 + k4);

// Обновляем следующее значение x

public static void main(String args[])

differential d2 = new differential();

double x0 = 0 , y = 1 , x = 2 , h = 0.2 ;

System.out.println( «\nThe value of y at x is : «

+ d2.rungeKutta(x0, y, x, h));

// Этот код предоставлен Prateek Bhindwar

# Программа Python для реализации метода Рунге Кутты
# Пример дифференциального уравнения «dy / dx = (x — y) / 2»

# Находит значение y для заданного x, используя размер шага h
# и начальное значение y0 при x0.

def rungeKutta(x0, y0, x, h):

# Подсчитать количество итераций, используя размер шага или

n = ( int )((x — x0) / h)

# Итерировать по количеству итераций

for i in range ( 1 , n + 1 ):

«Apply Runge Kutta Formulas to find next value of y»

k1 = h * dydx(x0, y)

k2 = h * dydx(x0 + 0.5 * h, y + 0.5 * k1)

k3 = h * dydx(x0 + 0.5 * h, y + 0.5 * k2)

k4 = h * dydx(x0 + h, y + k3)

# Обновить следующее значение y

y = y + ( 1.0 / 6.0 ) * (k1 + 2 * k2 + 2 * k3 + k4)

# Обновить следующее значение x

print ‘The value of y at x is:’ , rungeKutta(x0, y, x, h)

# Этот код предоставлен Prateek Bhindwar

// C # программа для реализации Runge
// метод Кутты

static double dydx( double x, double y)

// Находит значение y для данного x

// используя размер шага h и начальный

// значение y0 в x0.

static double rungeKutta( double x0,

double y0, double x, double h)

// Подсчитать количество итераций используя

// размер шага или высота шага h

int n = ( int )((x — x0) / h);

double k1, k2, k3, k4;

// Итерация по количеству итераций

for ( int i = 1; i

// Применяем формулы Рунге Кутты

// найти следующее значение у

k1 = h * (dydx(x0, y));

k2 = h * (dydx(x0 + 0.5 * h,

k3 = h * (dydx(x0 + 0.5 * h,

k4 = h * (dydx(x0 + h, y + k3));

// Обновить следующее значение y

y = y + (1.0 / 6.0) * (k1 + 2

// Обновляем следующее значение x

public static void Main()

double x0 = 0, y = 1, x = 2, h = 0.2;

Console.WriteLine( «\nThe value of y»

+ rungeKutta(x0, y, x, h));

// Этот код предоставлен Sam007.

// PHP-программа для реализации
// метод Рунге Кутта

// Пример дифференциального уравнения
// «dy / dx = (x — y) / 2»

function dydx( $x , $y )

return (( $x — $y ) / 2);

// Находит значение y для
// дано х, используя размер шага h
// и начальное значение y0 в x0.

function rungeKutta( $x0 , $y0 , $x , $h )

// Подсчитать количество итераций

// используя размер шага или шаг

$k1 ; $k2 ; $k3 ; $k4 ; $k5 ;

// Итерация по номеру

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

// Применить Рунге Кутта

// формулы для поиска

// следующее значение у

$k1 = $h * dydx( $x0 , $y );

$k2 = $h * dydx( $x0 + 0.5 * $h ,

$k3 = $h * dydx( $x0 + 0.5 * $h ,

$k4 = $h * dydx( $x0 + $h , $y + $k3 );

// Обновить следующее значение y

$y = $y + (1.0 / 6.0) * ( $k1 + 2 *

// Обновляем следующее значение x

echo «The value of y at x is : » ,

rungeKutta( $x0 , $y , $x , $h );

// Этот код предоставлен anuj_67.
?>

Сложность по времени вышеупомянутого решения составляет O (n), где n — (x-x0) / ч.

Некоторые полезные ресурсы для подробных примеров и большего количества объяснения.
http://w3.gazi.edu.tr/

Эта статья предоставлена Арпит Агарвал . Если вам нравится GeeksforGeeks и вы хотите внести свой вклад, вы также можете написать статью и отправить ее по почте на contrib@geeksforgeeks.org. Смотрите свою статью, появляющуюся на главной странице GeeksforGeeks, и помогите другим вундеркиндам.

Пожалуйста, напишите комментарии, если вы обнаружите что-то неправильное, или вы хотите поделиться дополнительной информацией по обсуждаемой теме

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

На практике наибольшее распространение получил метод Рунге-Кутта 4-го порядка, в котором усреднение проводится по трём точкам, формула Эйлера на каждом отрезке используется 4 раза: в начале отрезка, дважды в его середине и в конце отрезка.

Расчетные формулы метода для дифференциального уравнения (7.3) имеют вид:

, (7.8)

Где i = 0, 1, …., n-1 — номер узла;

Xi = a + i×h — координата узла;

У0 = у(х0) — начальное условие.

Погрешность метода dМ = О(h5).

Схема алгоритма решения ОДУ методом Рунге-Кутта 4-го порядка отличается алгоритмом расчёта новой точки (Рис. 7.5).

Пример 7.4. Решение ранее рассмотренного уравнения (пример 7.1) методом Рунге-Кутта 4 порядка.

Y’ — 2×y + x2 = 1, x Î [0;1], y(0) = 1.

Пусть n = 10 , h = (1 — 0)/10 = 0,1.

Начальная точка x0 = 0, y0 = 1.

Рассчет первой точки.

Сначала вычислим значения C0, C1, C2, C3:

Вычислим значение y1:

Аналогично можно вычислить значения функции во 2, 3, . , 10 точках.

Рис. 7.7. Схема алгоритма расчета новой точки методом Рунге-Кутта 4-го порядка.

Общая характеристика методов:

1. Все методы являются Одношаговыми, то есть для вычисления значения функции в новой точке используется ее значение в предыдущей точке. Это свойство называется Самостартованием.

2. Все методы легко обобщаются на системы дифференциальных уравнений 1-го порядка.