unit SODU;
{ модуль для решения системы ОДУ с первой производной  }
{ метод вычислений - Рунге-Кутта 2-го порядка точности }

{$N+}   {включаем сопроцессор}

interface

const
  Max_Vector = 64;   {максимальное количество элементов в векторе}
                     {или максимальное количество уравнений в системе ОДУ}

type
  TRealType = extended;  {задается тип данных для всего модуля}
                         {от этого зависит точность вычислений}
  TVector = record     {вектор значений}
    count :integer;
    X :array[1..Max_Vector] of TRealType;
  end;

  {решение системы ОДУ}
  Tsodu = object
    public
      N :integer;      {количество уравнений}
      X :TVector;      {текущее значение решения}
      t :TRealType;    {текущее значение аргумента}
      h :TRealType;    {текущее значение шага вычислений}
      eps :TRealType;  {точность вычислений}

      constructor Init;  {инициализация объекта}
      destructor Done;   {уничтожение объекта}
       {функция в которой задана система ОДУ}
      function F(tf:TRealType; Xf:TVector; var Yf:TVector):integer; virtual;
       {сам метод вычислений - Рунге-Кутта 2-го порядка точности}
      function Calc(t0, h0: TRealType; X0: TVector; var XX: TVector):TRealType; virtual;
       {пошаговое вычисление (без отслеживания точности eps)}
      procedure StepCalc;
       {вычисление с корректировкой шага h для достижения точности eps}
      procedure StepCalcAdapt; virtual;
       {непрерывное вычисление до аргумента t_end с заданной точностью eps}
      procedure FullCalc(t_end :TRealType);
  end;


implementation


function min(a,b:Trealtype):Trealtype;
begin
  if a>b then min:=b else min:=a;
end;

function max(a,b:Trealtype):Trealtype;
begin
  if a>b then max:=a else max:=b;
end;


{ Tsodu }

function Tsodu.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]:=-Xf.x[2]; {первое уравнение: X1'(t)=X2(t)  }
    Yf.x[2]:=Xf.x[1];  {второе уравнение: X2'(t)=-X1(t) }
  end;
end;

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;

destructor Tsodu.Done;
begin
  N:=0;
end;

constructor Tsodu.Init;
var
  Y: TVector;
  i: integer;
begin
  t:=0;
  X.count:=0;
  h:=0.001;
  eps:=0.001;
  N:=F(t,X,Y);
  X.count:=N;
  Y.count:=N;
  for i:=1 to N do X.x[i]:=0;
end;

procedure Tsodu.StepCalc;
begin
  eps:=Calc(t,h,X,X);  {просто считаем}
  t:=t+h;
end;

procedure Tsodu.StepCalcAdapt;
var  newX :TVector;
     error :TRealType;
begin
  if eps=0 then eps:=0.001;

  error:=Calc(t,h,X,newX);
  while error>eps do  {пока не достигнем требуемой точности}
  begin
    h:=h*min(1.6,max(0.2,0.9*exp((ln(eps/error))/4))); {уменьшаем шаг}
    error:=Calc(t,h,X,newX);
  end;

  X:=newX;
  t:=t+h;

  if (10*error)<eps then  {увеличим шаг если точность вычислений слишком высока для требуемой}
    h:=h*min(1.6,max(0.2,0.9*exp((ln(eps/error))/4)));
end;

procedure Tsodu.FullCalc(t_end :TRealType);
begin
  h:=t_end-t;
  if eps=0 then eps:=0.001;

  repeat
    if (t+h)>t_end then h:=t_end-t;   {поправляем - если перелетели за t_end}
    StepCalcAdapt;
  until t>(t_end-h/2);
end;

end.