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

численные методы решения дифференциальных уравнений Метод Эйлера — наиболее простой численный метод решения (систем) обыкновенных  дифференциальных уравнений. Впервые описан Леонардом Эйлером в 1768 году в  работе «Интегральное исчисление»[1]. Метод Эйлера является явным, одношаговым  методом первого порядка точности, основанном на аппроксимации интегральной  кривой кусочно­линейной функцией, т. н. ломаной Эйлера. Пусть дана задача Коши для уравнения первого порядка где функция  интервале   определена на некоторой области  . На этом интервале введем узлы . Решение ищется на  Приближенное решение в узлах  , которое обозначим через   определяется по формуле Эти формулы обобщаются на случай систем обыкновенных дифференциальных  уравнений. Оценка погрешности[править | править вики­текст] Метод Эйлера является методом первого порядка. Если функция  непрерывно дифференцируема по переменной  погрешности  в   непрерывна в   и  , то имеет место следующая оценка  где   — средний шаг, то есть существует   такая, что  . Заметим, что условия гладкости на правую часть, гарантирующие единственность  решения задачи Коши, необходимы для обоснования сходимости метода Эйлера. Значение метода Эйлера[править | править вики­текст] Метод Эйлера являлся исторически первым методом численного решения задачи  Коши. О. Коши использовал этот метод для доказательства существования решения  задачи Коши. Ввиду невысокой точности и вычислительной неустойчивости для практического нахождения решений задачи Коши метод Эйлера применяется редко.  Однако в виду своей простоты метод Эйлера находит свое применение в теоретических  исследованиях дифференциальных уравнений, задач вариационного исчисленияи ряда  других математических проблем. Модифицированный метод Эйлера с пересчетом[править | править вики­текст] Вычисления по методу Эйлера с пересчетом делаются в два этапа. Прогноз: Коррекция: . . Модифицированный метод Эйлера с пересчетом имеет второй порядок точности, однако для его реализации необходимо дважды вычислять правую часть функции. Заметим, что  метод Эйлера с пересчетом представляет собой разновидность методов Рунге­ Кутты (предиктор­корректор). См. также[править | править вики­текст] Дифференциальным уравнением первого порядка называется уравнение  вида F(x,y,у')=0 или у'=f(x,y). Функция y(x), при подстановке которой  уравнение обращается в тождество, называется решением  дифференциального уравнения. Рассмотрим несколько численных методов решения дифференциальных  уравнений первого порядка. Описание численных методов приводится  для уравнения в виде у'=f(x,y). 1. Метод Эйлера. Рассмотрим два варианта вывода расчетных формул  вариант 1 (аналитический) у=f (x,y) y1=y0+h*f(x0,y0) x1=x0+h yi+1=yi+h*f(xi,yi) xi+1=xi*h Расчетные формулы для 1­го  шага Расчетные формулы для i­го  шага  вариант 2 (графический) y1=y0+f(x0,y0)*h; x1=x0+h yi+1=yi+h*f(xi,yi) Аналогично варианту 1 k1=h*f(xi,yi) yi+1=yi+ki xi+1=xi+h  Следующие расчетные формулы приводятся без вывода. 2. Модифицированный метод Эйлера (вариант 1). уi+1=уi+hf(xi+h/2, yi+hf(xi,yi)/2), xi+1=xi+h. 3. Модифицированный метод Эйлера (вариант 2). уi+1=уi+(h/2)[f(xi,yi)+f(xi,+h,yi+hf(xi,yi))], xi+1=xi+h. 4. Метод Рунге­Кутта третьего порядка. уi+1=уi+(k1+4k2+k3)/6, k1=hf(xi, yi), k2=hf(xi+h/2, yi+k1/2), k3=hf(xi+h, yi+2k2­k1), xi+1=xi+h. 5. Метод Рунге­Кутта четвертого порядка. уi+1=уi+(k1+2k2+2k3+k4)/6, k1=hf(xi,yi), k2=hf(xi+h/2, yi+k1/2), k3=hf(xi+h/2, yi+k2/2), k4=hf(xi+h, yi+k3), xi+1=xi+h, где уi+1,уi ­ значения искомой функции в точках xi+1, xi соответственно,  индекс i показывает номер шага интегрирования, h ­ шаг  интегрирования. Начальные условия при численном интегрировании  учитываются на нулевом шаге: i=0, x=x0, y=y0. Пример. Численно и аналитически решить дифференциальное  уравнение dy/dx=x2 при y|x=0 =1. Определить значение функции при xk=1, h=1. Решение задачи приведено в таблице. Таблица N Этап  программирования Выполнение 1. Постановка задачи Решить дифференциальное уравнение dy/dx=x2 при y| x=0=1. Определить знач. функции при xk=1, h=1 2. Математическое  1. Аналитическое решение. описание dy/dx=x2 y=1+x3/3, yk=y(1)=1+1/3=4/3. 2. Метод Эйлера. 3. Модифицированный метод Эйлера 1. 4. Модифицированный метод Эйлера 2. 5. Метод Рунге­Кутта четвертого порядка. 3. 4. 5. Разработка  структограммы Написание  программы Отладка и  получение  результатов Выполнить самостоятельно Выполнить самостоятельно Выполнить самостоятельно Контрольное задание. Лабораторная работа 5. Численное решение дифференциальных уравнений Задание. 1. Решить дифференциальное уравнение аналитически и численно  указанными методами для двух значений шага интегрирования  h=0.01; 0.001. Результаты расчета вывести на экран и распечатать  в виде таблицы. 2. Построить графики функций y(x) (5 графиков). Варианты уравнений и методов их решения приведены в таблице Оформление результатов расчета Таблица Решения уравнения, у(x) Численное Аналит метод 1 Метод 2 h=0.01 h=0.001 h=0.01 h=0.001                     х       Варианты уравнений и методов их решения Таблица Вар. Вид уравнения Метод Вар. Вид уравнения Метод 1 2 3 4 5 6 7 8 у'=(xy2+x)/(y­x2y) у'=(1­2x)/y2 у'=(1­x2)/xy у'=(y2­y)/x y'=(1+y)/(tg(x) у'=exp(x)­1 y'=y ln(y)/sin(x) у'=(1+y2)/(1+x2) 1,4 2,4 3,4 1,5 2,5 3,5 1,4 2,4 14 15 16 17 18 19 20 21 у'=cos(t)­y y'=exp(bx)­ay У'=­2y/(y2­6x) у'=1/(2x­y2) у'=sec(x)­ y tg(x) y'=(exp(x)­y)/x у'=1+y/(x(x+1)) у'=(y+yx2­x2)/(x(1+x2)) 3,5 1,4 2,4 3,4 1,5 2,5 3,5 1,4 9 10 11 12 13 у'=4x­2y у'=x exp(­x2)­2xy у'=2x­y у'=exp(­x)­2y у'=exp(­x)­2x 3,4 1,5 2,5 3,5 1,4 22 23 24 25 26 у'=cos(x­y) у'=3x­2y+5 у'=sin(x)­y у'=exp(x)­y у'=exp(2x)­1 2,4 3,4 1,5 2,5 3,5 Примечание. Значение параметров a, b и начальные условия y| x=x0=y0 выбрать cамостоятельно. Содержание отчета: 1. Название, цель работы и задание. 2. Математическое описание, алгоритм (структограмма) и текст  программы. 3. Результаты расчета, пять графиков зависимости y(x) и выводы по  работе. Примеры реализации численных методов решения систем дифференциальных  уравнений Рассмотрим некоторые численные методы решения заданной системы.   Программная реализация (Turbo Pascal) program Euler; var   T : Text; a, b, h, x: real; y1, y2, y3, y4, y1_1, y1_2, y1_3, y1_4 : real; function F(x, y1, y2, y3, y4: real; n:byte):real; begin     case n of     1: f:=­y2;     2: f:=y1;     3: f:=y1­y4;     4: f:=y2+y3;     end; end; begin     assign(T, 'c:\euler.txt'); rewrite(T);     a:=0; b:=2*pi; h:=pi/40;     y1:=1; y2:=0; y3:=0; y4:=0;     x:=a;     while x<=b+h do     begin        writeln(T, x:10:5, y1:10:5, y2:10:5, y3:10:5, y4:10:5);        y1_1:=y1+h*f(x, y1, y2, y3, y4, 1);        y1_2:=y2+h*f(x, y1, y2, y3, y4, 2);        y1_3:=y3+h*f(x, y1, y2, y3, y4, 3);        y1_4:=y4+h*f(x, y1, y2, y3, y4, 4);        y1:=y1_1; y2:=y1_2; y3:=y1_3; y4:=y1_4;        x:=x+h;     end; flush(T); close(T); end. Программная реализация (Turbo Pascal) program Euler_Koshi; var T:Text;     x,h,a,b,y1, y2, y3, y4, y1_1, y1_2, y1_3,y1_4,     y11, y12, y13, y14:real; function F(x, y1, y2, y3, y4: real; n:byte):real; begin     case n of     1: f:=­y2;     2: f:=y1;     3: f:=y1­y4;     4: f:=y2+y3;     end; end; begin     assign(T, 'd:\euler_k.txt');     rewrite(T);     a:=0; b:=2*pi; h:=pi/160;     y1:=1; y2:=0; y3:=0; y4:=0; x:=a;     while x<=b+h do     begin        writeln(T, x:10:5, y1:10:5, y2:10:5, y3:10:5, y4:10:5,        cos(x):10:5, sin(x):10:5, x*cos(x):10:5, x*sin(x):10:5);        y1_1:=y1+h*f(x, y1, y2, y3, y4, 1);        y1_2:=y2+h*f(x, y1, y2, y3, y4, 2);        y1_3:=y3+h*f(x, y1, y2, y3, y4, 3);        y1_4:=y4+h*f(x, y1, y2, y3, y4, 4);        y11:=y1+h/2*(f(x, y1, y2, y3, y4, 1)+f(x+h,y1_1, y1_2, y1_3, y1_4, 1));        y12:=y2+h/2*(f(x, y1, y2, y3, y4, 2)+f(x+h,y1_1, y1_2, y1_3, y1_4, 2));        y13:=y3+h/2*(f(x, y1, y2, y3, y4, 3)+f(x+h,y1_1, y1_2, y1_3, y1_4, 3));        y14:=y4+h/2*(f(x, y1, y2, y3, y4, 4)+f(x+h,y1_1, y1_2, y1_3, y1_4, 4));        y1:=y11; y2:=y12; y3:=y13; y4:=y14;        x:=x+h;     end;     flush(T); close(T); end. Программная реализация (Turbo Pascal) program rk_III; var T:Text;     x,h,a,b,y1, y2, y3, y4,     k1_1, k1_2, k1_3, k1_4,     k2_1, k2_2, k2_3, k2_4,     k3_1, k3_2, k3_3, k3_4:real; function F(x, y1, y2, y3, y4: real; n:byte):real; begin     case n of     1: f:=­y2;     2: f:=y1;     3: f:=y1­y4;     4: f:=y2+y3;     end; end; begin     assign(T, 'd:\rk_iii.txt');     rewrite(T);     a:=0; b:=2*pi; h:=pi/160;     y1:=1; y2:=0; y3:=0; y4:=0;     x:=a;     while x<=b+h do     begin        writeln(T, x:10:5, y1:10:5, y2:10:5, y3:10:5, y4:10:5,        cos(x):10:5, sin(x):10:5, x*cos(x):10:5, x*sin(x):10:5); k1_1:=h*f(x, y1, y2, y3, y4, 1);        k1_2:=h*f(x, y1, y2, y3, y4, 2);        k1_3:=h*f(x, y1, y2, y3, y4, 3);        k1_4:=h*f(x, y1, y2, y3, y4, 4);        k2_1:=h*f(x+h/2, y1+k1_1/2, y2+k1_2/2, y3+k1_3/2, y4+k1_4/2, 1);        k2_2:=h*f(x+h/2, y1+k1_1/2, y2+k1_2/2, y3+k1_3/2, y4+k1_4/2, 2);        k2_3:=h*f(x+h/2, y1+k1_1/2, y2+k1_2/2, y3+k1_3/2, y4+k1_4/2, 3);        k2_4:=h*f(x+h/2, y1+k1_1/2, y2+k1_2/2, y3+k1_3/2, y4+k1_4/2, 4);        k3_1:=h*f(x+h,y1+k2_1, y2+k2_2, y3+k2_3, y4+k2_4, 1);        k3_2:=h*f(x+h,y1+k2_1, y2+k2_2, y3+k2_3, y4+k2_4, 2);        k3_3:=h*f(x+h,y1+k2_1, y2+k2_2, y3+k2_3, y4+k2_4, 3);        k3_4:=h*f(x+h,y1+k2_1, y2+k2_2, y3+k2_3, y4+k2_4, 4);        y1:=y1+1/4*(k1_1+2*k2_1+k3_1);        y2:=y2+1/4*(k1_2+2*k2_2+k3_2);        y3:=y3+1/4*(k1_3+2*k2_3+k3_3);        y4:=y4+1/4*(k1_4+2*k2_4+k3_4);        x:=x+h;     end;     flush(T); close(T); end. Программная реализация (Turbo Pascal) program rk_iv; var T:Text;     x,h,a,b,y1, y2, y3, y4,     k1_1, k1_2, k1_3, k1_4,     k2_1, k2_2, k2_3, k2_4,     k3_1, k3_2, k3_3, k3_4,     k4_1, k4_2, k4_3, k4_4:real; function F(x, y1, y2, y3, y4: real; n:byte):real; begin     case n of     1: f:=­y2;     2: f:=y1;     3: f:=y1­y4;     4: f:=y2+y3;     end; end; begin assign(T, 'd:\rk_iv.txt');     rewrite(T);     a:=0; b:=2*pi; h:=pi/80;     y1:=1; y2:=0; y3:=0; y4:=0;     x:=a;     while x<=b+h do     begin        writeln(T, x:10:5, y1:10:5, y2:10:5, y3:10:5, y4:10:5,        cos(x):10:5, sin(x):10:5, x*cos(x):10:5, x*sin(x):10:5);        k1_1:=h*f(x, y1, y2, y3, y4, 1);        k1_2:=h*f(x, y1, y2, y3, y4, 2);        k1_3:=h*f(x, y1, y2, y3, y4, 3);        k1_4:=h*f(x, y1, y2, y3, y4, 4);        k2_1:=h*f(x+h/2, y1+k1_1/2, y2+k1_2/2, y3+k1_3/2, y4+k1_4/2, 1);        k2_2:=h*f(x+h/2, y1+k1_1/2, y2+k1_2/2, y3+k1_3/2, y4+k1_4/2, 2);        k2_3:=h*f(x+h/2, y1+k1_1/2, y2+k1_2/2, y3+k1_3/2, y4+k1_4/2, 3);        k2_4:=h*f(x+h/2, y1+k1_1/2, y2+k1_2/2, y3+k1_3/2, y4+k1_4/2, 4);        k3_1:=h*f(x+h/2, y1+k2_1/2, y2+k2_2/2, y3+k2_3/2, y4+k2_4/2, 1);        k3_2:=h*f(x+h/2, y1+k2_1/2, y2+k2_2/2, y3+k2_3/2, y4+k2_4/2, 2);        k3_3:=h*f(x+h/2, y1+k2_1/2, y2+k2_2/2, y3+k2_3/2, y4+k2_4/2, 3);        k3_4:=h*f(x+h/2, y1+k2_1/2, y2+k2_2/2, y3+k2_3/2, y4+k2_4/2, 4);        k4_1:=h*f(x+h,y1+k3_1, y2+k3_2, y3+k3_3, y4+k3_4, 1);        k4_2:=h*f(x+h,y1+k3_1, y2+k3_2, y3+k3_3, y4+k3_4, 2);        k4_3:=h*f(x+h,y1+k3_1, y2+k3_2, y3+k3_3, y4+k3_4, 3);        k4_4:=h*f(x+h,y1+k3_1, y2+k3_2, y3+k3_3, y4+k3_4, 4); y1:=y1+1/6*(k1_1+2*k2_1+2*k3_1+k4_1);        y2:=y2+1/6*(k1_2+2*k2_2+2*k3_2+k4_2);        y3:=y3+1/6*(k1_3+2*k2_3+2*k3_3+k4_3);        y4:=y4+1/6*(k1_4+2*k2_4+2*k3_4+k4_4);        x:=x+h;     end;     flush(T); close(T); end.   ­ See more at: http://www.toehelp.ru/theory/informat/lecture13.html#sthash.vw5jkjpi.dpuf