Численные методы решения СОДУ, система обыкновенных дифф-ных уравнений Опции

[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.

Сам модуль:  SODU.PAS ( 4.41 килобайт ) Кол-во скачиваний: 1473


Сообщение отредактировано: volvo - 2.11.2006 22:21

[hiv]

Используем модуль  SODU.PAS ( 4.41 килобайт ) Кол-во скачиваний: 1473
для решения системы дифференциальных уравнений типа "хищник-жертва" методом Рунге-Кутты 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.