Решение дифференциальных уравнений четвертого порядка методом Эйлера и методом Рунге-Кутта

{$E+,N+,X+}
 
const
  eps = 1e-5;
  x0  = 3;
  y0  = 3;
 
type
  TFunc = function(x,y: extended): extended;
 
function func(x,y: extended): extended; far;
begin
  func := y*cos(x) - 2*sin(2*x);
end;
 
function Euler(f: TFunc; x0_,x_end,y0_: extended; n: word): extended;
{ где x0_ и y0_  - начальное условие
      x_end - точка, в которой необходимо вычислить результат
      n - количество шагов для вычисления результата }
var
  i   : word;      { счетчик цикла }
  x,h : extended;  { текущая точка и длина шага }
  res : extended;  { переменная для накопления конечного результата функции}
begin
  h:= (x_end - x0_)/n; { Находим длину шага }
  res:= y0_; { устанавливаем начальные значения}
  x:=x0_;
  for i:=1 to n do
    begin { вычисляем результат по методу Эйлера }
      res:=res+h*f(x,res);
      x:=x+h; { переходим к следующей точке }
    end;
  Euler:=res;  { присваиваем конечный результат функции }
end;
 
function Euler2(f: TFunc; x0_,x_end,y0_: extended; n: word): extended;
{ где x0_ и y0_  - начальное условие
      x_end - точка, в которой необходимо вычислить результат
      n - количество шагов для вычисления результата }
var
  i   : word;      { счетчик цикла }
  x,h : extended;  { текущая точка и длина шага }
  res : extended;  { переменная для накопления конечного результата функции}
begin
  h:= (x_end - x0_)/n; { Находим длину шага }
  res:= y0_; { устанавливаем начальные значения}
  x:=x0_;
  for i:=1 to n do
    begin { вычисляем результат по исправленному методу Эйлера }
      res:=res+h*(f(x,res)+f(x+h,res+h*f(x,res)))/2;
      x:=x+h; { переходим к следующей точке }
    end;
  Euler2:=res; { присваиваем конечный результат функции }
end;
 
function Euler3(f: TFunc; x0_,x_end,y0_: extended; n: word): extended;
{ где x0_ и y0_  - начальное условие
      x_end - точка, в которой необходимо вычислить результат
      n - количество шагов для вычисления результата }
var
  i   : word;      { счетчик цикла }
  x,h : extended;  { текущая точка и длина шага }
  res : extended;  { переменная для накопления конечного результата функции}
begin
  h:= (x_end - x0_)/n; { Находим длину шага }
  res:= y0_; { устанавливаем начальные значения}
  x:=x0_;
  for i:=1 to n do
    begin { вычисляем результат по модифицированному методу Эйлера }
      res:=res+h*f(x+h/2,res+(h/2)*f(x,res));
      x:=x+h; { переходим к следующей точке }
    end;
  Euler3:=res; { присваиваем конечный результат функции }
end;
 
function RungeKutt(f: TFunc; x0_,x_end,y0_: extended; n: word): extended;
{ где x0_ и y0_  - начальное условие
      x_end - точка, в которой необходимо вычислить результат
      n - количество шагов для вычисления результата }
var
  i   : word;      { счетчик цикла }
  x,h : extended;  { текущая точка и длина шага }
  res : extended;  { переменная для накопления конечного результата функции }
  k1,k2,k3,k4: extended; { вспомогательные переменные вычисления результата }
begin
  h:= (x_end - x0_)/n; { Находим длину шага }
  res:= y0_; { устанавливаем начальные значения}
  x:=x0_;
  for i:=1 to n do
    begin { вычисляем результат по методу Рунге-Кутта 4го порядка }
      k1:=f(x,res);
      k2:=f(x+h/2,res+h*k1/2);
      k3:=f(x+h/2,res+h*k2/2);
      k4:=f(x+h,res+h*k3);
      res:=res+h*(k1+2*k2+2*k3+k4)/6;
      x:=x+h; { переходим к следующей точке }
    end;
  RungeKutt:=res; { присваиваем конечный результат функции }
end;
 
function AdmsBoshf(f: TFunc; x0_,x_end,y0_: extended; n: word): extended;
{ где x0_ и y0_  - начальное условие
      x_end - точка, в которой необходимо вычислить результат
      n - количество шагов для вычисления результата }
var
  i     : word;         { счетчик цикла }
  x,h   : extended;     { текущая точка и длина шага }
  y1,y2 : extended;     { переменные для вычисления следующей точки }
  tmp   : extended;     { временная переменная для уточнения результата }
begin
  h:= (x_end - x0_)/n; { Находим длину шага }
  x:=x0_; { устанавливаем текущую точку }
  y1:=y0_+h*f(x+h/2,y0_+(h/2)*f(x,y0_)); { находим начальное значение
                                           модифицированным методом Эйлера }
  repeat { уточняем(корректируем) начальное значение }
    tmp:=y1;
    y1:=y0_+h*(f(x,y0_)+f(x+h,y1))/2;
  until abs(y1-tmp) < eps;
 
  x:=x+h; { переходим к следующей точке }
  y2:=y0_+2*h*f(x,y1); { прогнозируенм следующее значение }
 
  repeat { уточнаяем наш прогноз }
    tmp:=y2;
    y2:=y1+h*(f(x,y1)+f(x+h,y2))/2;
  until abs(y2-tmp) < eps;
 
  for i:=3 to n do { начинаем счет с 3 т.к. две точки уже найдены }
    begin
      x0_:=x; { переходим к следующей точке }
      x:=x+h;
      y0_:=y1;
      y1:=y2;
 
      y2:=y0_+2*h*f(x,y1); { прогнозируенм следующее значение }
 
      repeat { уточнаяем наш прогноз }
        tmp:=y2;
        y2:=y1+h*(f(x,y1)+f(x+h,y2))/2;
      until abs(y2-tmp) < eps;
    end;
 
  AdmsBoshf:=y2; { присваиваем конечный результат функции }
end;
 
begin
  writeln('Численное решение дифференциальных уравнений:');
  writeln(#10,'   y'' = y*cos(x) - 2*sin(2*x);   y(0)=3;   x_end=',(5*x0+3.5):5:5);
  writeln(#10,'Метод Эйлера:');
  writeln('n=5 000, result: ',Euler(func,x0,5*x0+3.5,y0,5000):5:5);
  writeln('n=10 000, result: ',Euler(func,x0,5*x0+3.5,y0,10000):5:5);
  writeln('n=25 000, result: ',Euler(func,x0,5*x0+3.5,y0,25000):5:5);
  writeln(#10,'Исправленный метод Эйлера:');
  writeln('n=50, result: ',Euler2(func,x0,5*x0+3.5,y0,50):5:5);
  writeln('n=100, result: ',Euler2(func,x0,5*x0+3.5,y0,100):5:5);
  writeln('n=250, result: ',Euler2(func,x0,5*x0+3.5,y0,250):5:5);
  writeln(#10,'Модифицированный метод Эйлера:');
  writeln('n=50, result: ',Euler3(func,x0,5*x0+3.5,y0,50):5:5);
  writeln('n=100, result: ',Euler3(func,x0,5*x0+3.5,y0,100):5:5);
  writeln('n=250, result: ',Euler3(func,x0,5*x0+3.5,y0,250):5:5);
 
  writeln(#10,'Метод Рунге-Кутта:');
  writeln('n=50, result: ',RungeKutt(func,x0,5*x0+3.5,y0,50):5:5);
  writeln('n=100, result: ',RungeKutt(func,x0,5*x0+3.5,y0,100):5:5);
  writeln('n=250, result: ',RungeKutt(func,x0,5*x0+3.5,y0,250):5:5);
 
  write(#10,'Press Enter to continue.');
  readln;
 
  writeln(#10,'Метод Адамса-Бошфора:');
  writeln('n=50, result: ',AdmsBoshf(func,x0,5*x0+3.5,y0,50):5:5);
  writeln('n=100, result: ',AdmsBoshf(func,x0,5*x0+3.5,y0,100):5:5);
  writeln('n=250, result: ',AdmsBoshf(func,x0,5*x0+3.5,y0,250):5:5);
end.

far

Euler(func, x0, 5*x0+3.5, y0, 5000)

Euler(func, 0, 1, 0, 10)

h:= (x_end - x0_)/n;

n:= round((x_end - x0_)/h);