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

Данная статья посвящена численным методам решения дифференциальных уравнений. Мы будем рассматривать методы решения одного обыкновенного дифференциального уравнения первого порядка с одним начальным условием:

	y' = f(x,y)							(1)
	y(x[0]) = y[0]							(2)

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

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

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

Методы Рунге-Кутта обладают следующими отличительными свойствами:

• Эти методы одноступенчатые: чтобы найти y[m+1], нужна информация только о

предыдущей точке x[m],y[m].

• Они согласуются с рядом Тейлора вплоть до членов порядка h^p,

где p - различна для разных методов и называется порядком метода.

• Они не требуют вычисления производных от f(x,y), а требуют только вычисления

самой функции.

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

Этот метод, один из самых старых и широко известных, описывается формулой:

	y[m+1] = y[m] + h*f(x[m],y[m]).					(3)

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

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

	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)

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

В данном методе мы находим tg угла наклона касательной в точке:

	x = x[m] + h/2;     y = y[m] + (h/2)*y'[m]

Соотношения, описывающие модифицированный метод Эйлера имеют вид:

	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)

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

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

	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)

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

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

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

Для прогноза используем формулу второго порядка

	y^(0)_[m+1] = y[m-1] + 2*h*f(x[m],y[m])				(15)

, где (0) - означает исходное приближение y[m+1], т.е. предсказанное значение. Непосредственно из написанной формулы следует, что с ее помощью нельзя вычислить y[1], т.к. для его вычисления потребовалась бы точка, расположенная перед начальной точкой y[0]. Чтобы начать решение с помощью метода прогноза и коррекции, для нахождения y[1] необходимо использовать метод типа Рунге-Кутта. Для коррекции возьмем формулу, похожую на исправленный метод Эйлера (4)-(6):

	y^(i)_[m+1] = y[m] + h(f(x[m],y[m])+f(x[m+1],y^(i-1)_[m+1]))/2	(16)

для i=1,2,3, …

Итерационный процесс прекращается, когда

	| y^(i+1)_[m+1] - y^(i)_[m+1] | < Eps				(17)

для некоторого Eps>0.

{$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.

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