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

Vladimir_S · 19.05.2014, 11:55

Поскольку на наш форум вообще и на меня в частности в апреле-мае 2014 года обрушилась лавина студенческих стенаний по поводу этой задачки, разберем ее подробно.
Иногда уравнение y'=f(x,y) решается аналитически, чаще же - нет, и тут на помощь приходят численные методы.
Начну с небольшого "экскурса в историю". Когда где-то 40 лет назад я был студентом, компьютеров не было (были только ЭВМ, занимавшие целые залы и, в общем, практически недоступные), а потому основным методом численного решения задачи был т.н. "метод изоклин". Суть его в том, что (тут и далее мы будем держать в голове геометрическую интерпретацию), зная производную y' в любой точке плоскости XY, мы можем изрисовать эту самую плоскость коротенькими отрезками прямой, такими, что тангенс угла наклона каждого из них есть f(x,y), а потом, вооружившись карандашом и опираясь на интуицию, изобразить семейство кривых-решений задачи так, чтобы эти отрезки были касательными в каждой точке и чтобы не было изломов.
Сегодня, ясное дело, такое привидится лишь в страшном сне, у всех есть компьютеры, все студенты блестяще владеют языками программирования и... ой, что это я? Ну ладно. В общем, вспомнили о численных методах, разработанных столетия назад великими математиками Эйлером, Рунге, Куттой и др. Основное достоинство этих методов - простота алгоритмизации. Вот их и рассмотрим.
Но прежде - очень важное замечание.
Для того, чтобы можно было применять указанные методы, одного уравнения мало. Обязательно требуется "опорная точка х0", в которой должно быть известно значение искомой функции y0(x0). Кроме того, следует задать шаг h, с которым мы будем двигаться вдоль оси Х. В общем случае этот шаг не обязательно должен быть константой, может быть и переменным, но мы ограничимся случаем постоянного шага.
Ну вот, теперь можно и приступать.

Метод Эйлера

Самый простой, грубый, но, в то же время, самый наглядный. Итак, представьте себе, что мы стартуем из точки х0, значение функции у0 в которой мы знаем. Знаем также и значение производной в точке (х0,у0), равное f(x0,y0). Но производная - это ни что иное, как тангенс угла наклона касательной в точке (х0,у0), а потому, предполагая, что на участке от x0 до x1=(x0+h) искомая функция не слишком сильно отклоняется от этой самой касательной, легко, просто из школьной геометрии, находим приращение функции при переходе от x0 к x1, или, окончательно,
y1 = y0 + h*f(x0,y0).
Дальше "опорной" становится точка (х1,у1), по которой мы находим у2, ну и т.д.
Разумеется, "ежу понятно", что здесь более, чем вероятно, накопление ошибки из-за слишком грубых допущений, а потому был изобретен т.н. "модифицированный" метод Эйлера с пересчетом (он же метод Рунге-Кутты второго порядка), несколько более сложный, но зато существенно повышающий точность. Как выяснилось (для меня недавно), есть несколько вариантов этого метода, различающихся способом внесения коррекций в процессе вычислений. Итак:

Модифицированный метод Эйлера с пересчетом (Рунге-Кутты второго порядка)

Вариант 1
Начинаем, как и в методе Эйлера, с определения у1, соответствующего искомой функции в точке x1=x0+h, но здесь рассматриваем y1, как предварительное значение. Обозначим его у1a.
Находим производную в точке (x1,y1a), которая, естественно, равна f(x1,y1a). Далее берем среднее арифметическое производных в точках (х0,у0) и (х1,у1а) и его-то и считаем "истинным" тангенсом угла наклона хорды, проведенной через точки (х0,у0) и искомую (х1,у1). Отсюда
y1a = x0 + h*f(x0,y0)
y1 = x0 + h*(f(x0,y0) + f(x1, y1a))/2
Дальше, естественно, по цепочке находим у2, у3 и т.д.

Хочу отметить, что именно этот алгоритм и рассматривается в учебниках и руководствах, как метод Рунге-Кутты с коррекцией, однако, как выяснилось, некоторые преподы под этим подразумевают несколько иной способ, который мы и рассмотрим. Итак,

Вариант 2
Здесь в качестве y1a выступает значение искомой функции для точки, находящейся на середине интервала между х0 и х1, т.е. для х0+h/2, каковое определяется, естественно, по производной f(x0,y0), т.е.
y1a = x0 + h/2*f(x0,y0)
Далее значение f(x0+h/2,y1a) мы принимаем в качестве "истинного" тангенса, и отсюда
y1 = x0 + h*f(x0+h/2,y1а)
Дальше - всё та же цепочка.

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

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

Здесь для получения у1 следует предварительно сосчитать 4 вспомогательных параметра K1, K2, K3 и K4, причем каждый следующий получается из предыдущего:
K1 = f(x0, y0)
K2 = f(x0+h/2, y0+h/2*K1)
K3 = f(x0+h/2, y0+h/2*K2)
K4 = f(x0+h, y0+h*K3)
И тогда у1 вычисляется, как
y1 = y0 + h/6*(K1 + 2*K2 + 2*K3 + K4)
Дальше, как и в предыдущих методах, идем по цепочке.

В следующем сообщении будет показана работа всех рассмотренных алгоритмов на конкретном примере с приложением программы на Паскале.

Vladimir_S · 19.05.2014, 13:05

(продолжение)

В качестве примера рассмотрим уравнение
y' = x + y +1
y(1) = 1
Оно имеет точное решение:
y = 4*Exp(x-1) - x - 2
каковое мы используем для сравнения с численными.
Ниже приведен листинг программы и результат. Столбцы соответствуют результатам решения по рассмотренным алгоритмам, в последнем даны точные значения.

Код:

Const
 h=0.1;
 N=10;

Var
 X,Y1,Y2,Y3,Y4:Array[0..N] of Real;
 i:Integer;

Function F(Xf,Yf:Real):Real;
begin
 F:=1.0+Xf+Yf;
end;

Function Y_acc(X_acc:Real):Real;
begin
 Y_acc:=4.0*Exp(X_acc-1)-(X_acc+2);
end;

Procedure Euler;
begin
 X[0]:=1;
 Y1[0]:=1;
 for i:=1 to N do
  begin
   X[i]:=X[i-1]+h;
   Y1[i]:=Y1[i-1]+h*F(X[i-1],Y1[i-1]);
  end;
end;

Procedure Runge_Kutt_2ord_Var1;
var Z:real;
begin
 X[0]:=1;
 Y2[0]:=1;
 for i:=1 to N do
  begin
   X[i]:=X[i-1]+h;
   Z:=Y2[i-1]+h*F(X[i-1],Y2[i-1]);
   Y2[i]:=Y2[i-1]+h*(F(X[i-1],Y2[i-1])+F(X[i],Z))/2;
  end;
end;

Procedure Runge_Kutt_2ord_Var2;
var Z:real;
begin
 X[0]:=1;
 Y3[0]:=1;
 for i:=1 to N do
  begin
   X[i]:=X[i]+h;
   Z:=Y3[i-1]+h/2*F(X[i-1],Y3[i-1]);
   Y3[i]:=Y3[i-1]+h*F(X[i-1]+h/2,Z);
  end;
end;

Procedure Runge_Kutt_4ord;
var K1,K2,K3,K4:Real;

function fK1(Xk,Yk:Real):Real;
begin
 fK1:=F(Xk,Yk);
end;

function fK2(Xk,Yk,Hk,Q:Real):Real;
begin
 fK2:=F(Xk+Hk/2,Yk+Hk/2*Q);
end;

function fK3(Xk,Yk,Hk,Q:Real):Real;
begin
 fK3:=F(Xk+Hk/2,Yk+Hk/2*Q);
end;

function fK4(Xk,Yk,Hk,Q:Real):Real;
begin
 fK4:=F(Xk+Hk,Yk+Hk*Q);
end;

begin
 X[0]:=1;
 Y4[0]:=1;
 for i:=1 to N do
  begin
   X[i]:=X[i-1]+h;
   K1:=fK1(X[i-1],Y4[i-1]);
   K2:=fK2(X[i-1],Y4[i-1],h,K1);
   K3:=fK3(X[i-1],Y4[i-1],h,K2);
   K4:=fK4(X[i-1],Y4[i-1],h,K3);
   Y4[i]:=Y4[i-1]+h/6*(K1+2.0*K2+2.0*K3+K4);
  end;
end;

Begin
 Euler;
 Runge_Kutt_2ord_Var1;
 Runge_Kutt_2ord_Var2;
 Runge_Kutt_4ord;
 Writeln(' X      Euler       R-K_2_V1    R-K_2_V2     R_K_4      Accur.');
 for i:=0 to N do
  begin
   write(X[i]:3:1);
   write(Y1[i]:12:5);
   write(Y2[i]:12:5);
   write(Y3[i]:12:5);
   write(Y4[i]:12:5);
   writeln(Y_acc(X[i]):12:5);
  end;
 ReadLn;
End.

Видно, что процедура Рунге-Кутты четвертого порядка дает практически точный результат.

AЛХИМИК · 07.06.2015, 17:33

Замечание:

Тема несет информационный характер. Не для обсуждения.
При желании задать вопрос - создавайте отдельную тему.

19.05.2014, 11:55	#1 (permalink)
Vladimir_S Специалист Регистрация: 27.08.2008 Адрес: Санкт-Петербург Сообщений: 27,807 Сказал(а) спасибо: 340 Поблагодарили 583 раз(а) в 208 сообщениях Репутация: 113184	Численные методы решения дифференциальных уравнений Поскольку на наш форум вообще и на меня в частности в апреле-мае 2014 года обрушилась лавина студенческих стенаний по поводу этой задачки, разберем ее подробно. Иногда уравнение y'=f(x,y) решается аналитически, чаще же - нет, и тут на помощь приходят численные методы. Начну с небольшого "экскурса в историю". Когда где-то 40 лет назад я был студентом, компьютеров не было (были только ЭВМ, занимавшие целые залы и, в общем, практически недоступные), а потому основным методом численного решения задачи был т.н. "метод изоклин". Суть его в том, что (тут и далее мы будем держать в голове геометрическую интерпретацию), зная производную y' в любой точке плоскости XY, мы можем изрисовать эту самую плоскость коротенькими отрезками прямой, такими, что тангенс угла наклона каждого из них есть f(x,y), а потом, вооружившись карандашом и опираясь на интуицию, изобразить семейство кривых-решений задачи так, чтобы эти отрезки были касательными в каждой точке и чтобы не было изломов. Сегодня, ясное дело, такое привидится лишь в страшном сне, у всех есть компьютеры, все студенты блестяще владеют языками программирования и... ой, что это я? Ну ладно. В общем, вспомнили о численных методах, разработанных столетия назад великими математиками Эйлером, Рунге, Куттой и др. Основное достоинство этих методов - простота алгоритмизации. Вот их и рассмотрим. Но прежде - очень важное замечание. Для того, чтобы можно было применять указанные методы, одного уравнения мало. Обязательно требуется "опорная точка х0", в которой должно быть известно значение искомой функции y0(x0). Кроме того, следует задать шаг h, с которым мы будем двигаться вдоль оси Х. В общем случае этот шаг не обязательно должен быть константой, может быть и переменным, но мы ограничимся случаем постоянного шага. Ну вот, теперь можно и приступать. Метод Эйлера Самый простой, грубый, но, в то же время, самый наглядный. Итак, представьте себе, что мы стартуем из точки х0, значение функции у0 в которой мы знаем. Знаем также и значение производной в точке (х0,у0), равное f(x0,y0). Но производная - это ни что иное, как тангенс угла наклона касательной в точке (х0,у0), а потому, предполагая, что на участке от x0 до x1=(x0+h) искомая функция не слишком сильно отклоняется от этой самой касательной, легко, просто из школьной геометрии, находим приращение функции при переходе от x0 к x1, или, окончательно, y1 = y0 + hf(x0,y0). Дальше "опорной" становится точка (х1,у1), по которой мы находим у2, ну и т.д. Разумеется, "ежу понятно", что здесь более, чем вероятно, накопление ошибки из-за слишком грубых допущений, а потому был изобретен т.н. "модифицированный" метод Эйлера с пересчетом (он же метод Рунге-Кутты второго порядка), несколько более сложный, но зато существенно повышающий точность. Как выяснилось (для меня недавно), есть несколько вариантов этого метода, различающихся способом внесения коррекций в процессе вычислений. Итак: Модифицированный метод Эйлера с пересчетом (Рунге-Кутты второго порядка)* Вариант 1 Начинаем, как и в методе Эйлера, с определения у1, соответствующего искомой функции в точке x1=x0+h, но здесь рассматриваем y1, как предварительное значение. Обозначим его у1a. Находим производную в точке (x1,y1a), которая, естественно, равна f(x1,y1a). Далее берем среднее арифметическое производных в точках (х0,у0) и (х1,у1а) и его-то и считаем "истинным" тангенсом угла наклона хорды, проведенной через точки (х0,у0) и искомую (х1,у1). Отсюда y1a = x0 + hf(x0,y0) y1 = x0 + h(f(x0,y0) + f(x1, y1a))/2 Дальше, естественно, по цепочке находим у2, у3 и т.д. Хочу отметить, что именно этот алгоритм и рассматривается в учебниках и руководствах, как метод Рунге-Кутты с коррекцией, однако, как выяснилось, некоторые преподы под этим подразумевают несколько иной способ, который мы и рассмотрим. Итак, Вариант 2 Здесь в качестве y1a выступает значение искомой функции для точки, находящейся на середине интервала между х0 и х1, т.е. для х0+h/2, каковое определяется, естественно, по производной f(x0,y0), т.е. y1a = x0 + h/2f(x0,y0) Далее значение f(x0+h/2,y1a) мы принимаем в качестве "истинного" тангенса, и отсюда y1 = x0 + hf(x0+h/2,y1а) Дальше - всё та же цепочка. Итак, со вторым порядком вроде как разобрались. Но далеко не всегда точности, даваемой даже применением методов "с коррекцией", достаточно. И тогда следует применить метод Рунге-Кутты четвертого порядка. Он, конечно, существенно сложнее, и предполагает использование некоей четырехступенчатой рекурсивной процедуры, но зато резко повышает точность. Не будем пытаться понять, почему значения вспомогательных функций и коэффициентов именно такие, а не другие - просто поверим выдающимся умам, которые это изобрели. Итак Метод Рунге-Кутты четвертого порядка Здесь для получения у1 следует предварительно сосчитать 4 вспомогательных параметра K1, K2, K3 и K4, причем каждый следующий получается из предыдущего: K1 = f(x0, y0) K2 = f(x0+h/2, y0+h/2K1) K3 = f(x0+h/2, y0+h/2K2) K4 = f(x0+h, y0+hK3) И тогда у1 вычисляется, как y1 = y0 + h/6(K1 + 2K2 + 2K3 + K4) Дальше, как и в предыдущих методах, идем по цепочке. В следующем сообщении будет показана работа всех рассмотренных алгоритмов на конкретном примере с приложением программы на Паскале.

19.05.2014, 11:55
Helpmaster Member Регистрация: 08.03.2016 Сообщений: 0	На форуме люди обсуждали что то схожее, посмотрите Решения дифференциальных уравнений Метод Ньютона для решения системы m нелинейных уравнений Turbo Pascal. Численные методы

Ads
Яндекс Member Регистрация: 31.10.2006 Сообщений: 40200 Записей в дневнике: 0 Сказал(а) спасибо: 0 Поблагодарили 0 раз(а) в 0 сообщениях Репутация: 55070