====== Методы решения дифференциальных уравнений ======

Данная статья посвящена численным методам решения дифференциальных уравнений.
Мы будем рассматривать методы решения одного обыкновенного дифференциального
уравнения первого порядка с одним начальным условием:
<code>	y' = f(x,y)							(1)
	y(x[0]) = y[0]							(2)</code>

Методы, которые здесь рассматриваются, легко обобщаются для системы уравнений
первого порядка. А уравнения высших порядков можно свести к системе уравнений
первого порядка.

=== В основном существуют два широких класса методов ===

  * Одноступенчатые методы, в которых используется только информация о самой кривой в одной точке и не производятся итерации. Один из методов - решение с помощью рядов Тейлора, но на практике он не слишком удобен для использования.\\ Практически удобные методы этого класса - методы //Рунге-Кутта//. Эти методы прямые (без итераций), но требуют многократных повторных вычислений функции.\\ При использовании данных методов трудно оценивать допускаемую ошибку.
  * Многоступенчатые методы, в которых следующую точку кривой можно найти, не производя так много повторных вычислений, но для достижения достаточной точности требуются итерации. Большинство методов этого класса называются методами прогноза и коррекции (метод //Адамса-Бошфора//).\\ Некоторые трудности, связанные с использованием итерационной процедуры и с использованием нескольких начальных точек уравновешиваются тем фактом, что оценку ошибки при использовании этого метода легко получить в качестве побочного продукта вычислений.

Как и во многих других случаях, эти два класса методов придется сочетать разумным образом, учитывая их достоинства и недостатки.


==== 1. Методы Рунге-Кутта ====

Методы Рунге-Кутта обладают следующими отличительными свойствами:
  * Эти методы одноступенчатые: чтобы найти //y[m+1]//, нужна информация только о
   предыдущей точке //x[m],y[m]//.
  * Они согласуются с рядом Тейлора вплоть до членов порядка //h^p//,
   где //p// - различна для разных методов и называется порядком метода.
  * Они не требуют вычисления производных от //f(x,y)//, а требуют только вычисления
   самой функции.

Именно благодаря //3)// эти методы удобны для практических вычислений, однако для
вычисления одной последующей точки решения нам придется вычислять //f(x,y)//
несколько раз при различных x и y.

=== Метод Эйлера ===

Этот метод, один из самых старых и широко известных, описывается формулой:
<code>	y[m+1] = y[m] + h*f(x[m],y[m]).					(3)</code>

Найденное по формуле (3) решение согласуется с разложением в ряд Тейлора вплоть до членов порядка h, т.е. данный метод является методом Рунге-Кутта первого порядка.\\ Этот метод имеет довольно большую ошибку приближения, кроме того, он очень часто оказывается неустойчивым - изначально малая ошибка (происходящая от приближения, округления или исходных данных) увеличивается с ростом x.\\ Для вычисления //y[m+1]// метод Эйлера использует наклон касательной только в точке //x[m],y[m]//. Этот метод можно усовершенствовать множеством различных способов.\\ Рассмотрим два из них.


=== Исправленный метод Эйлера ===

В исправленном методе Эйлера мы находим средний tg угла наклона касательной
для двух точек: ////x[m],y[m]// и //x[m+1],y[m]+h*y'[m]////. Соотношения, описывающие
данный метод, имеют вид:

<code>	y[m+1] = y[m] + h*F(x[m],y[m],h)				(4)
	F(x[m],y[m],h) = ( y'[m] + f(x[m]+h, y[m]+h*y'[m]) )/2		(5)
	y'[m] = f(x[m],y[m])						(6)</code>

Исправленный метод Эйлера согласуется с разложением в ряд Тейлора вплоть до
членов степени //h^2//, являясь, т.о. методом Рунге-Кутта второго порядка.


=== Модифицированный метод Эйлера ===

В данном методе мы находим tg угла наклона касательной в точке:
<code>	x = x[m] + h/2;     y = y[m] + (h/2)*y'[m]</code>

Соотношения, описывающие модифицированный метод Эйлера имеют вид:
<code>	y[m+1] = y[m] + h*F(x[m],y[m],h)				(7)
	F(x[m],y[m],h) = f( x[m]+h/2, y[m]+(h/2)*y'[m] )		(8)
	y'[m] = f(x[m],y[m])						(9)</code>

Модифицированный метод Эйлера также согласуется с разложением в ряд Тейлора
вплоть до членов степени //h^2//, и также является методом Рунге-Кутта второго
порядка.

=== Метод Рунге-Кутта четвертого порядка ===

Данный метод является одним из самых употребительных методов интегрирования
дифференциальных уравнений. Этот метод применяется настолько широко, что в
литературе его просто называют "методом Рунге-Кутта" без всяких указаний на его
тип или порядок. Этот классический метод описывается системой следующих пяти
соотношений:

<code>	y[m+1] = y[m] + h*(k[1] + 2*k[2] + 2*k[3] + k[4])/6		(10)
	k[1] = f(x[m], y[m])						(11)
	k[2] = f(x[m]+h/2, y[m]+h*k[1]/2)				(12)
	k[3] = f(x[m]+h/2, y[m]+h*k[2]/2)				(13)
	k[4] = f(x[m]+h, y[m]+h*k[3])					(14)</code>

Ошибка приближения для этого метода равна //e[t]=k*h^5//. Заметим, что при
использовании этого метода функцию необходимо вычислять четыре раза.

Один из серьезных недостатков методов Рунге-Кутта состоит в отсутствии
простых способов оценки их ошибки.

====	 2. Методы прогноза и коррекции ====

Отличительной чертой методов Рунге-Кутта является то, что при вычислении
следующей точки //(x[m+1],y[m+1])// используется информация только об одной
предыдущей точке //(x[m],y[m])//, но не нескольких. Кроме того, для методов
Рунге-Кутта отсутствуют достаточно простые способы оценки ошибки, что приводит
к необходимости рассмотрения некоторых дополнительных методов решения ДУ.\\ Отличительное свойство этих методов состоит в том, что с их помощью нельзя
начать решение уравнения т.к. в них необходимо использовать информацию о
предыдущих точках решения. Чтобы начать решение, имея только одну точку,
определяемую начальными условиями, или для того, чтобы изменить шаг //(h)//, необходим метод типа Рунге-Кутта. Поэтому приходится использовать разумное
сочетание этих двух методов.\\ Методы, которые мы рассмотрим, известны под общим названием методов прогноза и корректировки. Как ясно из названия вначале "предсказывается" значение
//y[m+1]//, а затем используется тот или иной метод его "корректировки".\\  Среди множества возможных формул прогноза и коррекции выберем по одному
примеру, применимому ко многим практическим задачам.

=== Метод Адамса-Бошфора ===

Для прогноза используем формулу второго порядка
<code>	y^(0)_[m+1] = y[m-1] + 2*h*f(x[m],y[m])				(15)</code>
, где (0) - означает исходное приближение //y[m+1]//, т.е. предсказанное значение.
Непосредственно из написанной формулы следует, что с ее помощью нельзя вычислить //y[1]//, т.к. для его вычисления потребовалась бы точка, расположенная перед
начальной точкой //y[0]//. Чтобы начать решение с помощью метода прогноза и коррекции, для нахождения //y[1]// необходимо использовать метод типа Рунге-Кутта.
Для коррекции возьмем формулу, похожую на исправленный метод Эйлера (4)-(6):
<code>	y^(i)_[m+1] = y[m] + h(f(x[m],y[m])+f(x[m+1],y^(i-1)_[m+1]))/2	(16)</code>
для i=1,2,3, ...

Итерационный процесс прекращается, когда
<code>	| y^(i+1)_[m+1] - y^(i)_[m+1] | < Eps				(17)</code>
для некоторого //Eps>0//.

<code pascal>
{$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.</code>

Автор статьи: Labinskiy Nikolay aka e-moe (c) 2005