06. Метод Рунге – Кутта
Метод Рунге – Кутта является более точным по сравнению с методом Эйлера.
Суть уточнения состоит в том, что искомое решение представляется в виде разложения в ряд Тейлора. (См. Формула Тейлора.)
Если в этой формуле ограничиться двумя первыми слагаемыми, то получим формулу метода Эйлера. Метод Рунге – Кутта учитывает четыре первых члена разложения.
.
В методе Рунге – Кутта приращения DYi предлагается вычислять по формуле:
Где коэффициенты Ki вычисляются по формулам:
Пример. Решить методом Рунге – Кутта дифференциальное уравнение 
при начальном условии у(0) = 1 на отрезке [0; 0,5] с шагом 0,1.
Для i = 0 вычислим коэффициенты Ki.
Последующие вычисления приводить не будем, а результаты представим в виде таблицы.
Решим этот же пример методом Эйлера.
Применяем формулу 
Производя аналогичные вычисления далее, получаем таблицу значений:
Применим теперь уточненный метод Эйлера.
Для сравнения точности приведенных методов численного решение данного уравнения решим его аналитически и найдем точные значения функции У на заданном отрезке.
Уравнение 
является линейным неоднородным дифференциальным уравнением первого порядка. Решим соответствующее ему однородное уравнение.
Решение неоднородного уравнения имеет вид 
Общее решение: 
C учетом начального условия: 
Частное решение: 
Для сравнения полученных результатов составим таблицу.
Уточненный метод Эйлера
Метод Рунге — Кутта
Как видно из полученных результатов метод Рунге – Кутта дает наиболее точный ответ. Точность достигает 0,0001. Кроме того, следует обратить внимание на то, ошибка (расхождение между точным и приближенным значениями) увеличивается с каждым шагом вычислений. Это обусловлено тем, что во – первых полученное приближенное значение округляется на каждом шаге, а во – вторых – тем, что в качестве основы вычисления принимается значение, полученное на предыдущем шаге, т. е. приближенное значение. Таким образом происходит накопление ошибки.
Это хорошо видно из таблицы. С каждым новым шагом приближенное значение все более отличается от точного.
При использовании кмпьютерной версии “Курса высшей математики” возможно запустить программу, которая решает любое дифференциальное уравнение первого порядка рассмотренным выше методом Рунге — Кутта. Программа подробно выводит результаты вычислений на каждом шаге.
 |
Для запуска программы дважды щелкните на значке
Примечание: Для запуска программы необходимо чтобы на компьютере была установлена программа Maple (Ó Waterloo Maple Inc.) любой версии, начиная с MapleV Release 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 0 + h
, Меньший размер шага означает большую точность.
Формула в основном вычисляет следующее значение y n + 1, используя текущее значение y n плюс средневзвешенное значение четырех приращений.
- k 1 — приращение, основанное на наклоне в начале интервала, используя y
- K 2 представляет собой приращение на основе наклона в средней точке интервала, с использованием Y + 1 HK / 2.
- к 3 снова приращение на основе наклона в средней точке, используя при помощи у + кк 2/2.
- k 4 — это приращение, основанное на наклоне в конце интервала, с использованием y + hk 3 .
Этот метод является методом четвертого порядка, это означает, что локальная ошибка усечения имеет порядок 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, и помогите другим вундеркиндам.
Пожалуйста, напишите комментарии, если вы обнаружите что-то неправильное, или вы хотите поделиться дополнительной информацией по обсуждаемой теме
http://codenet.ru/progr/alg/Runge-Kutt-Method/
http://espressocode.top/runge-kutta-4th-order-method-solve-differential-equation/