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

[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 килобайт ) Кол-во скачиваний: 2117


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