Метод Рунге-Кутта решения диф. уравнений и их систем.
Метод позволяет решать системы обыкновенных дифференциальных уравнений (ОДУ) первого порядка следующего вида:
которые имеют решение:
где 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;
Метод Рунге — Кутты
Этот онлайн калькулятор реализует классический метод Рунге — Кутты (встречается также название метод Рунге — Кутта) четвертого порядка точности. Метод используется для решения дифференциальных уравнений первой степени с заданным начальным значением
Калькулятор ниже находит численное решение дифференциального уравнения первой степени методом Рунге-Кутты (иногда встречается название метод Рунге-Кутта, а в поисковиках бывает ищут «метод рунге кута», «метод рунги кутта» и даже «метод рунги кута»), который также известен как классический метод Рунге — Кутты (потому что есть на самом деле семейство методов Рунге-Кутты) или метод Рунге — Кутты четвертого порядка.
Для того, чтобы использовать калькулятор, вам надо привести дифференциальное уравнение к форме
и ввести правую часть уравнения f(x,y) в поле y’ калькулятора.
Также вам понадобится ввести начальное значение
и указать точку в которой вы хотите получить численное решение уравнения .
Последнее параметр калькулятора — размер шага с которым вычисляется следующее приближение по графику функции.
Описание метода можно найти под калькулятором.
Решение ОДУ в Matlab
Доброго времени суток! Сегодня мы поговорим о решении ОДУ (обыкновенных дифференциальных уравнений) в Matlab. Перед тем как мы начнём обсуждать данную тему, советую вам ознакомиться с темой: Численное дифференцирование в Matlab, чтобы лучше понимать теоретическую составляющую решения ОДУ.
Обыкновенные дифференциальные уравнения
С помощью дифференциальных уравнений можно описать разные задачи: движения системы, взаимодействующих материальных точек, химической кинетики и т.д. Различают три типа задач для систем диф. уравнений:
- Задача Коши
- Краевая задача
- Задача на собственные значения
Кратко расскажу о их сути:
Задача Коши предполагает дополнительные условия в виде значения функции в определённой точке.
Краевая задача подразумевает поиск решения на заданном отрезке с краевыми (граничными) условиями в концах интервала или на границе области.
Задача на собственные значения — помимо искомых функций и их производных, в уравнение входят дополнительное несколько неизвестных параметров, которые являются собственными значениями.
Методы решения дифференциальных уравнений
Решение ОДУ в Matlab и не только, в первую очередь, сводится к выбору порядка численного метода решения. Порядок численного метода не связан с порядком дифференциального уравнения. Высокий порядок у численного метода означает его скорость сходимости.
В случае большого интервала, с помощью алгоритмов с низким порядком сжимают интервал с решениями и находят приблизительные корни, а затем уже уточняют корни с помощью методов с высоким порядком.
Решение обыкновенных дифференциальных уравнений в Matlab можно реализовать «своими ручками», прописав алгоритм по разным схемам. Но также в Matlab есть встроенные функции, выполняющие все стандартные задачи.
Метод Рунге-Кутта первого порядка
Методы Рунге-Кутта представляют собой разложения в ряд Тейлора и от количества использованных элементов ряда зависит порядок этого метода. Следовательно, помимо Рунге-Кутта первого порядка, вы сможете увидеть методы других порядков. Иногда их называют другими именами.
Например, Метод Рунге-Кутта первого порядка, также известен как Метод Эйлера или Метод ломаных. Информацию о его математическом и графическом представлении советую поискать в гугл. Мы же поговорим о том, как Метод Рунге-Кутта первого порядка реализуется в Matlab для решения ОДУ. Например:
Решить и привести график ошибки уравнения y’ = y*x методом Рунге-Кутта первого порядка (Методом Эйлера, Методом ломаных).
Погрешность Метода Рунге-Кутта 1 порядка
» data-medium-file=»https://i2.wp.com/codetown.ru/wp-content/uploads/2017/02/Рунге-1-погрешность.png?fit=300%2C236&ssl=1″ data-large-file=»https://i2.wp.com/codetown.ru/wp-content/uploads/2017/02/Рунге-1-погрешность.png?fit=622%2C489&ssl=1″ loading=»lazy» src=»https://i2.wp.com/codetown.ru/wp-content/uploads/2017/02/%D0%A0%D1%83%D0%BD%D0%B3%D0%B5-1-%D0%BF%D0%BE%D0%B3%D1%80%D0%B5%D1%88%D0%BD%D0%BE%D1%81%D1%82%D1%8C.png?resize=622%2C489″ alt=»Погрешность метода 1 порядка» width=»622″ height=»489″ srcset=»https://i2.wp.com/codetown.ru/wp-content/uploads/2017/02/Рунге-1-погрешность.png?w=629&ssl=1 629w, https://i2.wp.com/codetown.ru/wp-content/uploads/2017/02/Рунге-1-погрешность.png?resize=300%2C236&ssl=1 300w» sizes=»(max-width: 622px) 100vw, 622px» data-recalc-dims=»1″ />
На данном графике показана зависимость величины ошибки от шага.
Метод Рунге-Кутта второго порядка
Также известен как Метод Эйлера-Коши. Как видите, во второй части уравнения происходит обращения к следующему шагу. Но как тогда быть, если нам ещё не известен следующий шаг? Всё просто. Метод Рунге-Кутта второго порядка — это всё тот же метод первого порядка, однако, на половине шага происходит нахождение «первичного» решения, а затем происходит его уточнение. Это позволяет поднять порядок скорости сходимости до двух.
Решить и привести график ошибки уравнения u’ = u*x методом Рунге-Кутта второго порядка.
По сравнению с Рунге-Куттом первого порядка изначальная ошибка уже гораздо меньше.
Мы не будем говорить о третьем порядке, потому что задачи на третий порядок встречаются редко, но если будет необходимо, пишите в комментариях, выложу.
Метод Рунге-Кутта четвёртого порядка
Метод Рунге-Кутта четвёртого порядка считается самым распространённым. Тем не менее, работает он аналогично второму и третьему порядку.
Решить и привести график ошибки уравнения u’ = u*x методом Рунге-Кутта четвёртого порядка.
Как видите, на последней картинке размерность ошибки на столько мала, что пришлось воспользоваться loglog() для лучшей видимости.
Решение ОДУ в Matlab стандартными средствами
Стоит отметить, что мы с вами разобрали только один самый известный метод решения ОДУ с разными порядками. Однако, методов очень много.
Для решения дифференциальных уравнений и систем в MATLAB предусмотрены следующие функции:
ode45 (f, interval, X0, [options])
ode23 (f, interval, X0, [options])
ode113 (f, interval, X0, [options])
ode15s (f, interval, X0, [options])
ode23s (f, interval, X0, [options])
ode23t (f, interval, X0, [options])
ode23tb (f, interval, X0, [options])
Входными параметрами этих функций являются:
- f — вектор-функция для вычисления правой части уравнения системы уравнений;
- interval — массив из двух чисел, определяющий интервал интегрирования дифференциального уравнения или системы;
- Х0 — вектор начальных условий системы дифференциальных уравнений;
- options — параметры управления ходом решения дифференциального уравнения или системы.
Все функции возвращают:
- массив Т — координаты узлов сетки, в которых ищется решение;
- матрицу X, i-й столбец которой является значением вектор-функции решения в узле Тi.
В функции ode45 реализован метод Рунге-Кутта 4-5 порядка точности, в функции ode23 также реализован метод Рунге-Кутта, но 2-3 порядка, а функция ode113 реализует метод Адамса.
Для решения жёстких систем предназначены функция ode15s, в которой реализован метод Гира, и функция ode23s, реализующая метод Розенброка. Для получения более точного решения жёсткой системы лучше использовать функцию ode15s. Для решения системы с небольшим числом жёсткости можно использовать функцию ode23t, а для грубой оценки подобных систем служит функция ode23tb.
Символьное решение обыкновенных дифференциальных уравнений произвольного порядка осуществляет функция dsolve r = dsolve(‘eq1,eq2,…’, ‘cond1,cond2,…‘, ‘v’)
Пример использования:
На этом мы закончим. Если остались вопросы, задавайте их в комментариях. Также вы можете скачать исходники чтобы лучше понять тему: «Решение ОДУ в Matlab».
http://planetcalc.ru/8400/
http://codetown.ru/matlab/reshenie-odu/