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

1.

ЧИСЛЕННЫЕ МЕТОДЫ
РЕШЕНИЯ
ДИФФЕРЕНЦИАЛЬНЫХ
УРАВНЕНИЙ
-Постановка задачи
-Метод Эйлера
-Модифицированный метод Эйлера

2.

Численные методы решения дифференциальных уравнений
Общий вид обыкновенного дифференциального уравнения, устанавливающего связь между
независимой переменной x неизвестной функцией y и ее производными y’,y”,…,y(n), может
мыть представлен следующим образом:
F( x , y, y , y ,.....y ( n ) ) 0
Порядок наивысшей производной, входящей в уравнение, называется порядком этого
уравнения. Решение дифференциального уравнения (интегрированием) является
некоторая функциональная зависимость y=y(x), которая при подстановке в уравнение
обращает его в тождество.
Общее решение дифференциального уравнения записывается в виде:
y=y(x,c1,c2,…,cn),
где c1,c2,…,cn произвольные постоянные.
Решение, полученное из общего решения при фиксированных значениях, называется
частным решением уравнения. Постоянные c1,c2,…,cn можно определить, задав n
условий. Если эти условия заданы как совокупность значений искомой функции и всех ее
производных до (n-1)ого порядка включительно в некоторой течке x0, то задача решения
уравнения называется задачей Коши, а заданные условия:
y(x0)=y0, y’(x0)=y’0, y”(x0)=y”0,…, yn-1(x0)=yn-10 называются начальными условиями.
Если условия заданы при разных значениях x, то задача решения дифференциального
уравнения называется граничной или краевой задачей.
2

3.

Рассмотрим дифференциальное уравнение первого порядка:
F( x, y, y ) 0
соотношение часто удается записать в виде:
y f ( x, y)
Последнее уравнение называется дифференциальным уравнением, разрешенным
относительно производной. Значение производной равно тангенсу угла наклона
касательной к графику функции в точке (x, y). Функцию f (x, y) будем называть правой
частью дифференциального уравнения.
Общим решением уравнения будет являться семейство функций y=y(x,c1)
различающихся значение постоянной c1. Задав, одно начальное условие y(x0)=y0
определяющее значение c1, получаем конкретное частное решение – задача Коши.
ПРИМЕР. Для дифференциального уравнения y’=3x2, общее решение имеет вид y=x3+c.
Подставим в общее решение начальное условие при x0=1, y0=2, вычислим с=1 и
определим частное решение как: y=x3+1
3

4.

-3
-2
10
9
8
7
6
5
4
3
y0 2
1
0
-1
-1
-2 0
-3
-4
-5
-6
-7
-8
-9
x0
1
2
3
Метод Эйлера
Дано дифференциальное уравнение y’=f (x, y), удовлетворяющее начальному условию
y(x0)=y0 . Требуется найти решение на отрезке [a,b]. Разобьем отрезок интегрирования на n
равных частей: x0=a, x1= a + h, x2=x1+h,…,xi=xi–1+h,…, xn=b, тогда величина шага
b a
интегрирования будет равна:
h
n
4

5.

Метод Эйлера
Постановка задачи.
Дано дифференциальное уравнение y’=f (x, y),
удовлетворяющее начальному условию y(x0)=y0
Требуется найти решение на отрезке [a,b]. Разобьем
отрезок интегрирования на n равных частей: x0=a,
x1= a + h, x2=x1+h,…,xi=xi–1+h,…, xn=b, тогда
величина шага интегрирования будет равна:
h=(b-a) / n

6.

y
y2
Значение функции y1 в
y1
точке x1 можно определить
как точку пересечения
касательной проведенной к
функции y=y(x) в точке
(x0,y0) с вертикальной
прямой проходящей через
y0
точку x1.
y=y(x)
β
x
x0
x1
x2
xi
xn-1
xn
Тангенс угла наклона касательной есть значение производной в точке (x0,y0) и задается
правой частью дифференциального уравнения, т.е. tg( β)=f(x0,y0). С другой стороны из
геометрического представления метода можно записать:
y y
y y
y y
tg( ) 1 0 1 0 Следовательно 1 0 f ( x 0 , y 0 )
h
x1 x 0
h
y1 y 0 h f ( x 0 , y 0 ) y 2 y1 h f ( x1, y1 ) и т.д.
Откуда
x 2 x1 h
x1 x 0 h
Решение будет заключаться в последовательном применении формул:
yi yi 1 h f ( x i 1, yi 1 )
x i x i 1 h
где i = 1, 2, 3, …, n
Результат будет представлен функцией заданной таблицей.
6

7.

Пример
y
C
y ; a 1; b 3; n 4; h 0.5 x 0 1; y 0 2; y
C 2
x
x
y1 = -2+0.5*(-(-2/1)) = -1
x1 = 1+0.5 =1.5
y2 = -1+0.5*(-(-1/1.5)) = -0.667 x2 = 1.5+0.5 =2
и т.д.
X
Y(Эйлер)
Y(теор)
1
-2
-2
1.5
-1
-1.333
2
-0.667
-1
2.5
-0.5
-0.8
3
-0.4
-0.667
7

8.

Модифицированный метод Эйлера
Графическая интерпретация.
h
Определяем точку x 1 x 0
2
2
и вычисляем значение
функции в этой точке
h
y 1 y0 f ( x 0 , y0 )
2
2
Значение функции y1 в точке x1 определяем, как точку пересечения касательной,
вычисленной в точке (x1/2,y1/2) и проведенной к функции y=y(x) в точке (x0,y0) , с
вертикальной прямой проходящей через точку x1.
y
y=y(x)
y1/2
y1
β
y0
x
x0 x
1/2
x1
x2
xi
xn-1
xn
h
h
y1 y0 h f ( x 1 , y 1 ) или y1 y 0 h f ( x 0 , y 0 f ( x 0 , y0 ))
2
2
2
2
x1 x 0 h
8

9.

произвольную точку определим
yi yi 1 h f ( x i 1/ 2 , yi 1/ 2 ) или
x i x i 1 h
h
h
yi 1 h f ( x i 1 , yi 1 f ( x i 1, yi 1 ))
2
2
где i = 1, 2, 3, …, n
Пример
y
C
y ; a 1; b 3; n 4; h 0.5 x 0 1; y 0 2; y
C 2
x
x
y1/2 = -2+0.25*(-(-2/1)) = -1.5; x1/2 = 1+0.25 = 1.25
y1 = -2+0.5*(-(-1.5/1.25)) = -1.4; x1 = 1+0.5 = 1.5
y3/2 = -1.4 +0.25*(-(-1.4/1.5)) = -1.1667;
x3/2 = 1.5+0.25 = 1.75
y2 = -1.4 +0.5*(-(-1.1667/1.75)) = -1.0667; y2 = 1.5+0.5 = 2
и т.д.
X
Y(мод)
Y(теор)
1
-2
-2
1.5
-1.4
-1.3333
2
-1.0667
-1
2.5
-0.8593
-0.8
3
-0.7187
-0.6667
9

10.

0
0
1
2
3
4
-0.5
-1
-1.5
Y(Эйлер)
Y(Мод)
Y(теор)
-2
-2.5
10

11.

Begin
euler (x, y, h)
x0,y0,xn,n
h:=(xn-x0)/n
x(1):=x0; y(1):=y0
i=0
function y=euler(x,y,h)
function y=m_euler(x,y,h)
dfdx(x,y)
dfdx(x,y)
y=y+h*dfdx(x,y)
y=y+h*dfdx(x+h/2,y+h/2*dfdx(x,y))
End
End
function y=us_euler(x,y, h)
function y=rungekutt(x,y,h)
dfdx(x,y)
dfdx(x,y)
y=y+h/2*(dfdx(x,y)+dfdx(x+h,y+h*dfdx(x,y))
i=2 step 1 to n+1
End
y(i)=euler(x(i-1),y(i-1),h)
x(i):=x(i-1)+h
plot(x,y)
End
k0:=dfdx(x,y)
k1:=dfdx(x+h/2,y+h*k0/2)
k2:=dfdx(x+h/2,y+h*k1/2)
k3:=dfdx(x+h,y+h*k2)
y=y+h/6*(k0+2*k1+2*k2+k3)
End
function f=dfdx(x,y)
dfdx = 2 * (x ^ 2 + y)
End
11

12.

Аналитический вывод формул
dy
f ( x , y)
dx
Необходимо найти значения функции y(x) в заданных точках x1, x2,… xn, если
известны начальные значения (x0, y0 ), где y0=y(x0) . Преобразуем уравнение
dy f ( x, y)dx
Проинтегрируем левую и правую части уравнения между xi и xi+1 точкой
yi 1
x i 1
x i 1
dy f ( x, y)dx
yi
yi 1 yi f ( x, y)dx
xi
xi
Интегрируем методом прямоугольники вперед
x i 1
yi 1 yi f ( x i , yi )dx yi f ( x i , yi ) ( x i 1 x i ) yi f ( x i , yi ) h
xi
Интегрируем методом в среднем
x i 1
yi 1 yi f ( x i 1/ 2 , yi 1/ 2 )dx yi f ( x i 1/ 2 , yi 1/ 2 ) h
xi
x i 1/ 2 x i h / 2
yi 1/ 2 yi h / 2 f ( x i , yi )
Интегрируем методом трапеций
yi 1 yi h / 2 (f ( x i , yi ) f ( x i 1, yi 1 ))
x i 1 x i h
yi 1 yi h f ( x i , yi )
12

13.

Система дифференциальных уравнений
y1 ' f ( x, y1, y 2 , y3 ,..., y n )
y 2 ' f ( x, y1, y 2 , y3 ,..., y n )
y n ' f ( x, y1, y 2 , y3 ,..., y n )
где x — независимая переменная, а y1(x), y2(x), ..., yn(x) — неизвестные функции, n —
порядок системы.
y1 ' ( x )
y2 ' (x)
Y'
y ' ( x )
n
Обозначив
f1 ( x, y1, y 2 , y n
f ( x , y , y , y
1 2
n
F ( x, Y ) 2
f ( x , y , y , y
n
1 2
n
или
Y' F ( x, Y )
Решением системы называется вектор-функция Y ( x ), которая определена и непрерывно
дифференцируема на интервале (a, b) и удовлетворяет системе, т.е. для всех x o (a, b)
справедливо
Y' F ( x, Y )
дифференциальные уравнения n-ого порядка
y n f ( x , y, y' , y' ' ,..., y n 1 )
Приводим к системе дифференциальных уравнений
y y1
y1 ' y 2
y 2 ' y3
............
y n ' f ( x, y1, y 2 , y3 ,..., y n )
13