Помощь - Поиск - Пользователи - Календарь
Полная версия: Численные методы решения СОДУ
Форум «Всё о Паскале» > Pascal, Object Pascal > Задачи > FAQ
hiv
Модуль для решения системы обыкновенных дифференциальных уравнений методом Рунге-Кутты 2-го порядка точности.
Сам метод реализован в функции:
function Tsodu.Calc(t0, h0: TRealType; X0: TVector; var XX: TVector):TRealType;
const {константы метода}
c1:TRealType=1;
c2:TRealType=0.5;
a10:TRealType=1;
a20:TRealType=0.25;
a21:TRealType=0.25;
b0:TRealType=0.5;
b1:TRealType=0.5;
bp0:TRealType=1/6;
bp1:TRealType=1/6;
bp2:TRealType=2/3;
var
k0,k1,k2,Xh: TVector;
error,err,xe,hh,th: TRealType;
i :integer;
begin
error:=-1;
if (N>0) then
begin
Xh.count:=N;
hh:=h0;
F(t0,X0,k0);
for i:=1 to N do Xh.x[i]:=X0.x[i]+hh*a10*k0.x[i];
th:=t0+hh*c1;
F(th,Xh,k1);
for i:=1 to N do Xh.x[i]:=X0.x[i]+hh*(a20*k0.x[i]+a21*k1.x[i]);
th:=t0+hh*c2;
F(th,Xh,k2);
for i:=1 to N do
begin
{вычисляем приближенное решение}
xe:=X0.x[i]+hh*(b0*k0.x[i]+b1*k1.x[i]);
{вычисляем более точное решение}
Xh.x[i]:=X0.x[i]+hh*(bp0*k0.x[i]+bp1*k1.x[i]+bp2*k2.x[i]);
if abs(xe)<1 {вычисляем ошибку}
then err:=ABS(Xh.x[i]-xe) {абсолютная ошибка}
else err:=ABS((Xh.x[i]-xe)/Xh.x[i]); {относительная ошибка}
if i=1 then error:=err
else if error<err then error:=err;
end;
XX:=Xh;
end;
Calc:=error;
end;

Пример использования модуля:
program test_sodu;
uses sodu;

var D :Tsodu;
i:integer;

begin
with D do
begin
Init;
X.x[1]:=1; {задаем начальные условия Коши}
X.x[2]:=0;
t:=0;
h:=0.15;
eps:=0.000001;
writeln('Решение:');
writeln(t:12:8,' ',X.x[1]:12:8,' ',X.x[2]:12:8,' ',cos(t):12:8,' ',sin(t):12:8,' ',(sin(t)-X.x[2]):12:8);
for i:=1 to 21 do
begin
FullCalc(i*0.15);
writeln(t:12:8,' ',X.x[1]:12:8,' ',X.x[2]:12:8,' ',cos(t):12:8,' ',sin(t):12:8,' ',(sin(t)-X.x[2]):12:8);
end;
Done;
end;
end.

Сам модуль: Нажмите для просмотра прикрепленного файла
hiv
Используем модуль Нажмите для просмотра прикрепленного файла для решения системы дифференциальных уравнений типа "хищник-жертва" методом Рунге-Кутты 5-го порядка точности:
program test_RKF5;
uses sodu;

type
TRKF5 = object (Tsodu)
{определяем новый метод вычислений - Рунге-Кутта 5-го порядка точности}
function Calc(t0, h0: TRealType; X0: TVector; var XX: TVector):TRealType; virtual;
{задаем новую систему ОДУ}
function F(tf:TRealType; Xf:TVector; var Yf:TVector):integer; virtual;
end;

function TRKF5.Calc(t0, h0: TRealType; X0: TVector; var XX: TVector):TRealType;
const {константы метода}
c1:Trealtype=1/5;
c2:Trealtype=0.3;
c3:Trealtype=4/5;
c4:Trealtype=8/9;
c5:Trealtype=1;
c6:Trealtype=1;
a10:Trealtype=1/5;
a20:Trealtype=3/40;
a21:Trealtype=9/40;
a30:Trealtype=44/45;
a31:Trealtype=-56/15;
a32:Trealtype=32/9;
a40:Trealtype=19372/6561;
a41:Trealtype=-25360/2187;
a42:Trealtype=64448/6561;
a43:Trealtype=-212/729;
a50:Trealtype=9017/3168;
a51:Trealtype=-355/33;
a52:Trealtype=46732/5247;
a53:Trealtype=49/176;
a54:Trealtype=-5103/18656;
a60:Trealtype=35/384;
a61:Trealtype=0;
a62:Trealtype=500/1113;
a63:Trealtype=125/192;
a64:Trealtype=-2187/6784;
a65:Trealtype=11/84;
b0:Trealtype=35/384;
b1:Trealtype=0;
b2:Trealtype=500/1113;
b3:Trealtype=125/192;
b4:Trealtype=-2187/6784;
b5:Trealtype=11/84;
bp0:Trealtype=5179/57600;
bp1:Trealtype=0;
bp2:Trealtype=7571/16695;
bp3:Trealtype=393/640;
bp4:Trealtype=-92097/339200;
bp5:Trealtype=187/2100;
bp6:Trealtype=1/40;

var
k0,k1,k2,k3,k4,k5,k6,Xh: TVector;
error,err,xe,hh,th: TRealType;
i :integer;
begin
error:=-1;
if (N>0) then
begin
Xh.count:=N;
hh:=h0;
F(t0,X0,k0);
for i:=1 to N do Xh.x[i]:=X0.x[i]+hh*a10*k0.x[i];
th:=t0+hh*c1;
F(th,Xh,k1);
for i:=1 to N do Xh.x[i]:=X0.x[i]+hh*(a20*k0.x[i]+a21*k1.x[i]);
th:=t0+hh*c2;
F(th,Xh,k2);
for i:=1 to N do Xh.x[i]:=X0.x[i]+hh*(a30*k0.x[i]+a31*k1.x[i]+a32*k2.x[i]);
th:=t0+hh*c3;
F(th,Xh,k3);
for i:=1 to N do Xh.x[i]:=X0.x[i]+hh*(a40*k0.x[i]+a41*k1.x[i]+a42*k2.x[i]+a43*k3.x[i]);
th:=t0+hh*c4;
F(th,Xh,k4);
for i:=1 to N do Xh.x[i]:=X0.x[i]+hh*(a50*k0.x[i]+a51*k1.x[i]+a52*k2.x[i]+a53*k3.x[i]+a54*k4.x[i]);
th:=t0+hh*c5;
F(th,Xh,k5);
for i:=1 to N do Xh.x[i]:=X0.x[i]+hh*(a60*k0.x[i]+a61*k1.x[i]+a62*k2.x[i]+a63*k3.x[i]+a64*k4.x[i]+a65*k5.x[i]);
th:=t0+hh*c6;
F(th,Xh,k6);
for i:=1 to N do
begin
{вычисляем приближенное решение}
xe:=X0.x[i]+hh*(bp0*k0.x[i]+bp1*k1.x[i]+bp2*k2.x[i]+bp3*k3.x[i]+bp4*k4.x[i]+bp5*k5.x[i]+bp6*k6.x[i]);
{вычисляем более точное решение}
Xh.x[i]:=X0.x[i]+hh*(b0*k0.x[i]+b1*k1.x[i]+b2*k2.x[i]+b3*k3.x[i]+b4*k4.x[i]+b5*k5.x[i]);
if abs(xe)<1 {вычисляем ошибку}
then err:=ABS(Xh.x[i]-xe) {абсолютная ошибка}
else err:=ABS((Xh.x[i]-xe)/Xh.x[i]); {относительная ошибка}
if i=1 then error:=err
else if error<err then error:=err;
end;
XX:=Xh;
end;
Calc:=error;
end;

function TRKF5.F(tf:TRealType; Xf:TVector; var Yf:TVector):integer;
const
FN=2; {количество уравнений в системе}
var
i :integer;
begin
F:=FN;
if X.count=FN then
begin
Yf.count:=FN; {сисема ОДУ "хищник-жертва"}
Yf.x[1]:=4*Xf.x[1]-2.5*Xf.x[1]*Xf.x[2]; {первое уравнение: X1'(t)=4*X1(t)-2.5*X1(t)*X2(t) }
Yf.x[2]:=Xf.x[1]*Xf.x[2]-2*Xf.x[2]; {второе уравнение: X2'(t)=X1(t)*X2(t)-2*X2(t) }
end;
end;


var D :TRKF5;
i:integer;

begin
with D do
begin
Init;
X.x[1]:=2; {задаем начальные условия Коши}
X.x[2]:=1;
t:=0;
h:=0.1;
eps:=0.000001;
writeln('Решение:');
writeln(t:12:8,' ',X.x[1]:12:8,' ',X.x[2]:12:8);
for i:=1 to 21 do
begin
FullCalc(i*0.1);
writeln(t:12:8,' ',X.x[1]:12:8,' ',X.x[2]:12:8);
end;
Done;
end;
end.
Это текстовая версия — только основной контент. Для просмотра полной версии этой страницы, пожалуйста, нажмите сюда.