Реализация метода Рунге-Кутты

program CallExponent;
{программа вычисляет значение константы e путём
 решения дифференциального уравнения y'=y, y(0)=1}
uses
  RungeKutta;
 
function Equation(x,y:extended):extended;
begin
  Result:=y
end;
 
begin
  WriteLn(GetSolutionAtPoint(1,0.001,0,1,Equation)) 
  // 2.71828182845902
  // погрешность в 14-ом знаке после запятой
end.

unit RungeKutta;
 
interface
 
{.$DEFINE Cache}
 
uses
  System.SysUtils, System.Math{$IFDEF Cache}, CacheUnit{$ENDIF};
 
type
  TRKEquation=reference to function(x,y:extended):extended;
 
  /// <summary>
  /// Решение ОДУ методом Рунге-Кутты
  /// </summary>
  IRungeKutta=interface
    /// <summary>
    /// Собственно решение
    /// </summary>
    function y(x:extended):extended;
  end;
 
/// <summary>
/// Создание интерфейса <see cref="ACod.Math.RungeKutta.IRungeKutta"/>
/// </summary>
function InitialRK(h,x0,y0:extended; RKEquation:TRKEquation
  {$IFDEF Cache}; Cache:TCacheLimit=0{$ENDIF}):IRungeKutta;
/// <summary>
/// Получение решения уравнения в точке X с шагом H
/// </summary>
function GetSolutionAtPoint(x,h,x0,y0:extended; RKEquation:TRKEquation)
  :extended;
 
implementation
 
type
  TRungeKutta=class(TInterfacedObject,IRungeKutta)
  private
    FH:extended;
    FX0:extended;
    FY0:extended;
    FRKEquation:TRKEquation;
    FYN,FXN:extended;
    FStepCount:UInt64;
    {$IFDEF Cache}
    FCache:boolean;
    FCacheMass:TCache<extended,extended>;
    {$ENDIF}
    procedure SetH(const Value:extended);
    procedure SetRKEquation(const Value:TRKEquation);
    procedure SetX0(const Value:extended);
    procedure SetY0(const Value:extended);
    procedure Step(dx:extended);
  public
    constructor Create{$IFDEF Cache}(Cache:TCacheLimit){$ENDIF};
    destructor Free;
    property RKEquation:TRKEquation read FRKEquation write SetRKEquation;
    property x0:extended read FX0 write SetX0;
    property y0:extended read FY0 write SetY0;
    property h:extended read FH write SetH;
    function y(x:extended):extended;
  end;
 
{ TRungeKutta }
 
constructor TRungeKutta.Create{$IFDEF Cache}(Cache:TCacheLimit){$ENDIF};
begin
  {$IFDEF Cache}
  FCache:=Cache>0;
  if FCache then
    FCacheMass:=TCache<extended,extended>.Create(Cache);
  {$ENDIF}
end;
 
destructor TRungeKutta.Free;
begin
  {$IFDEF Cache}
  if FCache then
    FCacheMass.Free;
  {$ENDIF}
end;
 
procedure TRungeKutta.SetH(const Value:extended);
begin
  FH:=Value;
end;
 
procedure TRungeKutta.SetRKEquation(const Value:TRKEquation);
begin
  FRKEquation:=Value;
end;
 
procedure TRungeKutta.SetX0(const Value:extended);
begin
  FX0:=Value;
end;
 
procedure TRungeKutta.SetY0(const Value:extended);
begin
  FY0:=Value;
end;
 
procedure TRungeKutta.Step(dx:extended);
var
  k1,k2,k3,k4:extended;
begin
  k1:=dx*RKEquation(FXN,FYN);
  k2:=dx*RKEquation(FXN+dx/2,FYN+k1/2);
  k3:=dx*RKEquation(FXN+dx/2,FYN+k2/2);
  k4:=dx*RKEquation(FXN+dx,FYN+k3);
  FYN:=FYN+(k1+2*k2+2*k3+k4)/6;
  FXN:=FXN+dx;
end;
 
function TRungeKutta.y(x:extended):extended;
var
  tmp:UInt64;
  pls:boolean;
begin
  Result:=NaN;
  if IsNaN(x)or IsInfinite(x) then
    raise Exception.CreateFmt('Неподходящее для вычислений значение x=%f',[x]);
  {$IFDEF Cache}
  if FCache then
    if FCacheMass.HasItem(x) then
      Exit(FCacheMass[x]);
  {$ENDIF}
  FXN:=FX0;
  FYN:=FY0;
  FStepCount:=round(abs(x-FX0)/FH);
  pls:=x<FX0;
  if pls then
    FH:=-FH;
  tmp:=0;
  while tmp<(FStepCount-1) do
    begin
      Step(FH);
      inc(tmp);
    end;
  {new:23.08.2012 10:37}
  Step(x-(FH*(FStepCount-1)+FX0));
  {/new}
  Result:=FYN;
  if pls then
    FH:=-FH;
  {$IFDEF Cache}
  if FCache then
    FCacheMass[x]:=FYN;
  {$ENDIF}
end;
 
{ Funtions }
 
function InitialRK(h,x0,y0:extended; RKEquation:TRKEquation{$IFDEF Cache}; Cache:TCacheLimit{$ENDIF})
  :IRungeKutta;
var
  tmp:TRungeKutta;
begin
  if CompareValue(h,0)<=0 then
    raise Exception.CreateFmt('Величина шага (%f) должна быть больше 0!',[h]);
  tmp:=TRungeKutta.Create{$IFDEF Cache}(Cache){$ENDIF};
  tmp.h:=h;
  tmp.x0:=x0;
  tmp.y0:=y0;
  tmp.RKEquation:=RKEquation;
  Result:=tmp;
end;
 
function GetSolutionAtPoint(x,h,x0,y0:extended; RKEquation:TRKEquation)
  :extended;
begin
  Result:=InitialRK(h,x0,y0,RKEquation{$IFDEF Cache},0{$ENDIF}).y(x);
end;