Численное решение задачи Коши

Лабораторная работа № 6.

Численное решение задачи Коши.


Основные понятия.

Уравнения первого порядка имеют вид , где f(x,y) – непрерывная функция в области D от x,y.
Пусть точка M0(x0,y0) принадлежит области D.

Задача Коши:

Найти такое решение y=y(x) этого уравнения, которое удовлетворяет условию y(x0)= y0, т.е. найти ту интегральную кривую, которая проходит через точку M0.


Постановка задачи Коши при численном решении дифференциальных уравнений.

Пусть дано и . Требуется найти приближённое решение поставленной задачи на отрезке [a,b], где a=x0.
y

| | | | | | x
a=x0 b
Разбиваем отрезок [a,b] равноотстоящими точками , где – шаг разбиения.
Ищем приближённое значение решения задачи Коши в точках .
Если найдено, то следующие значение получаем по формуле . в каждом методе находится по формулам, соответствующим выбранному методу.









… …



Решение получаем в виде таблицы









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


В методе Эйлера:
Рассмотрим метод Эйлера:
Пусть – точное решение задачи Коши.
В точке .
Известно, что - угловой коэффициент касательной к интегральной кривой y(x) в точке .
y



y(x)

h





x



С другой стороны:
Тогда:
Таким образом, вместо искомой интегральной кривой на отрезке , рассмотрим отрезок касательной к в точке .В случае небольшого шага график функции и график касательной не успевают существенно разойтись друг от друга и можно в качестве значения решения в принять значение касательной вместо значения неизвестного точного решения. В результате многократного повторения этого действия на каждом отрезке , получаем ломанную (см. рис.).


Метод Эйлера является простейшим методом и даёт грубое приближение. Более того, этот метод даёт систематическое отклонение искомого решения при переходе от одного шага к следующему. Наиболее приемлемым для практики методом оценки точности является в данном случае способ двойного счета – с шагом и с шагом .



ПРИМЕР 1.
Изучить работу метода Эйлера на отрезке [0,2] на примере задачи Коши вида

С этой целью найти точное решение задачи средствами Maple и приближенное решение с шагом 0.1, 0.02, 0.01.
Результаты сравнить графически. Сделать выводы.
Решение.
> restart:with(plots):
1. Зададим массивы для численного решения
> y1:=array(1..200,[]);x1:=array(1..200,[]);


Зададим правую часть дифференциального уравнения.
> f:=(x,u)->2*x*u;

Зададим дифференциальное уравнение.
> ode:=diff(y(x),x)=f(x,y(x));

Зададим начальное условие.
> iq:=y(0)=1;

Найдем точное решение дифференциального уравнения.
> dsolve({ode,iq},y(x));

Зададим выражение для построения графика.
> f1:=plot(exp(x^2),x=0..2,y=0..20):
Найдем приближенное решение с шагом 0.1 по формулам

и

> h:=.1;y1[1]:=1:x1[1]:=0:

> for i from 1 to 20 do y1[i+1]:=y1[i]+h*f(x1[i],y1[i]):x1[i+1]:=x1[i]+h:end do:
Зададим выражение для построения графика решения
> l := [[ x1[m], y1[m]] $m=1..20]:f2:=plot(l,x=0..2,y=0..20,style=point,symbol=circle,color=black):
Найдем приближенное решение с шагом 0.02 и получим выражение для построения графика.
> h:=.02;y1[1]:=1:x1[1]:=0: for i from 1 to 100 do y1[i+1]:=y1[i]+h*f(x1[i],y1[i]):x1[i+1]:=x1[i]+h:end do: l := [[ x1[m1], y1[m1]] $m1=1..100]:
f3:=plot(l,x=0..2,y=0..20,style=point,symbol=circle,color=blue):

Найдем приближенное решение с шагом 0.01 и получим выражение для построения графика.
> h:=.01;y1[1]:=1:x1[1]:=0: for i from 1 to 199 do y1[i+1]:=y1[i]+h*f(x1[i],y1[i]):x1[i+1]:=x1[i]+h:end do: l := [[ x1[m1], y1[m1]] $m1=1..200]:
f4:=plot(l,x=0..2,y=0..20,style=point,symbol=circle,color=green):

Построим графики всех решений.
> plots[display]({f1,f2,f3,f4});


Выводы:
1. Решение с большим шагом может приводить к большим погрешностям;
2. Уменьшение шага на порядок приводит к резкому уменьшению погрешности результатов.
3. Для достаточно малого шага отличие точного и приближенного решений незначительное.

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

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

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

Рассмотрим геометрический смысл



Сначала с шагом находим yi+h/2 т.е. точку (xi+h/2 , y i+h/2) (как в методе Эйлера). Определяем в этой точке новое направление . Затем проводим отрезок из точки Mi именно в этом новом направлении в точку Mi+1
Точность полученного метода примерно в 3 раза выше обычного метода Эйлера.



ПРИМЕР 2.

На примере предыдущего задания изучить работу модифицированного метода Эйлера, сравнивая точное решение задачи с приближенным, подсчитанным соответственно с шагом h=0.2, 0.1, 0.05.
Решение.
> restart: with(plots):

Зададим массивы для численного решения
> y1:=array(1..200,[]);x1:=array(1..200,[]);


Зададим правую часть дифференциального уравнения.
> f:=(x,u)->2*x*u;

Зададим дифференциальное уравнение.
> ode:=diff(y(x),x)=f(x,y(x));

Зададим начальное условие.
> iq:=y(0)=1;

Найдем точное решение дифференциального уравнения.
> dsolve({ode,iq},y(x));
f1:=plot(exp(x^2),x=0..2,y=0..20):

Найдем приближенное решение с шагом 0.2, по формулам

и

> h:=.2;y1[1]:=1:x1[1]:=0:

> for i from 1 to 20 do y1[i+1]:=y1[i]+h*f(x1[i]+h/2,y1[i]+f(x1[i],y1[i])*h/2):x1[i+1]:=x1[i]+h:end do:
Зададим выражение для построения графика решения
> l := [[ x1[m], y1[m]] $m=1..20]:
f2:=plot(l,x=0..2,y=0..20,style=point,symbol=circle,color=black):
Аналогично найдем приближенное решение с шагом 0.1и получим выражение для построения графика.
> h:=.1;y1[1]:=1:x1[1]:=0: for i from 1 to 100 do y1[i+1]:=y1[i]+h*f(x1[i]+h/2,y1[i]+f(x1[i],y1[i])*h/2):x1[i+1]:=x1[i]+h:end do: l := [[ x1[m1], y1[m1]] $m1=1..100]:
f3:=plot(l,x=0..2,y=0..20,style=point,symbol=circle,color=blue):

Найдем приближенное решение с шагом 0.05 и получим выражение для построения графика.
> h:=.05;y1[1]:=1:x1[1]:=0: for i from 1 to 199 do y1[i+1]:=y1[i]+h*f(x1[i]+h/2,y1[i]+f(x1[i],y1[i])*h/2):x1[i+1]:=x1[i]+h:end do: l := [[ x1[m1], y1[m1]] $m1=1..200]:
f4:=plot(l,x=0..2,y=0..20,style=point,symbol=circle,color=green):

Построим графики всех решений.
> plots[display]({f1,f2,f3,f4});





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

Метод Рунге – Кутта четвертого порядка является одним из самых широко используемых методов решения задачи Коши.

Пусть дано и отрезок , , – шаг разбиения.

Вычисления проводят по схеме:

Для оценки погрешности решения при переходе от точки xi к точке xi+1 выбирают значение шага h и дважды рассчитывают решение в точке xi+1 : один раз с шагом h, другой раз – с шагом 2h . Рассчитывается величина
и сравнивается с заданной погрешностью . Если , то можно продолжить вычисления с тем же шагом, в противном случае необходимо вернуться к решению в точке xi,уменьшив шаг и повторить вычисления.
В Maple 6 реализован метод Рунге – Кутта 4 порядка method=rkf45 (установлен по умолчанию). Для того, чтобы найти численное решение задачи Коши или краевой задачи в команде dsolve следует указать параметр numeric . Тогда команда решения дифференциального уравнения будет иметь вид
ode:= dsolve(eq,vars, type=numeric, method=rkf45): где eq - уравнения, vars - список неизвестных функций. График численного решения можно построить с помощью команды odeplot(ode,[x,y(x)],x1..x2);
ПРИМЕР 3.

на отрезке [0,2]
Решение.

> restart:with(plots):
Warning, the name changecoords has been redefined

 ode:=dsolve({diff(y(x),x)=2*x*y(x),y(0)=1},y(x),numeric);

В строке вывода появилось сообщение о том, что при решении использован метод Рунге – Кутта (rkf45) . Если требуется получить решение при фиксированном значении x, например, при x=0,1, то следует набрать:

> ode(.1);

Построим график численного решения

 odeplot(ode,[x,y(x)],0..2);




КОНТРОЛЬНОЕ ЗАДАНИЕ 1.
Решить задачу Коши для обыкновенного дифференциального уравнения первого порядка методом Эйлера и модифицированным методом Эйлера на отрезке[0,2;1,2] с шагом 0,1 при начальном условии y(0,2)=0,25. Построить графики решений.
КОНТРОЛЬНОЕ ЗАДАНИЕ 2.
Решить задачу Коши для обыкновенного дифференциального уравнения первого порядка методом Рунге-Кутта четвертого порядка на отрезке[0,2;1,2] при начальном условии y(0,2)=0,25
1) используя в Maple 6 method=rkf45. Построить график решения.
2) Оценить погрешность решения в точке х=1,2.
Замечание.
Условие вашего задания получить у преподавателя.