Эксперименты лаба10,11(2курс)
.pdf1. Обыкновенные дифференциальные уравнения
ОДУ называются такие уравнения, которые содержат одну или несколько производных от искомой функции y = y(x):
F (x, y, y0, ..., y(n)) = 0. |
(1) |
Наивысший порядок n входящей в уравнение производной называется порядком дифференциального уравнения. Уравнение первого порядка
F (x, y, y0) = 0,
второго порядка
F (x, y, y0, y00) = 0.
Если из уравнения удается выразить старшую производную и привести
его к виду
y(n) = F (x, y, y0, ..., y(n−1)),
то такая форма записи называется уравнением, разрешенным относительно старшей производной.
Линейным дифференциальным уравнением называется уравнение, линейное относительно искомой функции и ее производных.
Решением дифференциального уравнения называется всякая n раз дифференцируемая функция y = ϕ(x), которая после подстановки в уравнение превращает его в тождество.
Общее решение уравнения n-го порядка содержит n произвольных постоянных C1, C2, ... , Cn:
y = ϕ(x, C1, C2, ..., Cn). |
(2) |
Причем любое решение нашего уравнения можно представить в виде (2) при определенных значениях констант.
Для уравнения первого порядка общее решение зависит от одной произвольной постоянной:
y = ϕ(x, C).
Если постоянной придать определенное значение C = C0, то получаем частное решение
y = ϕ(x, C0).
Геометрическая интерпретация общего решения уравнения первого порядка состоит в том, что оно описывает бесконечное семейство интегральных кривых y = y(x, C) с параметром C, а частному решению
соответствует одна кривая этого семейства, причем через каждую точку (x0, y0) проходит одна и только одна интегральная кривая.
Для уравнений высших порядков – через каждую точку (x0, y0) проходит не одна интегральная кривая, поэтому для выделения некоторого частного решения надо задать дополнительные условия по числу произвольных постоянных.
1.1. Приближенные методы решения ОДУ. Метод Эйлера
Наиболее универсальным методом решения ОДУ является метод конечных разностей. Он заключается в том, что область непрерывного изменения аргумента заменяется дискретным множеством точек – узлами, образующими разностную сетку. Искомая функция непрерывного аргумента заменяется функцией дискретного аргумента на данной сетке – сеточной функцией. Решение ОДУ сводится к отысканию значений сеточной функции в узлах сетки.
Например, дифференциальное уравнение (задача Коши)
dy |
= f(x, y), y(x0) = A |
dx |
вводя равномерную сетку с шагом h и приняв в качестве узлов сетки x0, x1 = x0 + h,..., можно свести к разностному, приближая производную в i-ом узле с помощью правых разностей
yi+1 − yi |
= fi, i = 0, 1, ..., y0 = A. |
|
h |
||
|
Здесь fi = f(xi, yi). Таким образом мы определили разностную схему. Разностной схемой называется замкнутая система разностных уравнений вместе с дополнительными условиями – начальными или краевыми.
Выражая yi+1, получаем
yi+1 = yi + hfi = yi + hf(xi, yi), i = 0, 1, ... |
(3) |
В результате приходим к методу Эйлера:
y1 = y0 + hf(x0, y0), y2 = y1 + hf(x1, y1),
...........
yn = yn−1 + hf(xn−1, yn−1).
2
Очевидно, что в данном случае искомая интегральная кривая приближается ломаной, наклон которой на элементарном участке [xi, xi+1] определяется наклоном интегральной кривой исходного уравнения, выходящей из точки (xk, yk).
Неявный метод Эйлера состоит в приближении производной в окрестности i-го узла с помощью левой разности
yi − yi−1 |
= fi, i = 0, 1, ..., y0 = A. |
h |
|
Очевидно, что в этом методе искомая величина yi входит в уравнение в общем случае нелинейным образом. Поэтому необходимо применять известные методы решения нелинейных уравнений.
Можно показать, что явный и неявный метод – методы первого порядка точности. То есть погрешность в точке xi
δi = y(xi) − yi,
равная разности между точным значением искомой функции y(xi) и значением сеточной функции в узле
δi . h.
Можно построить другие методы решения задачи Коши. Например, приблизим производную в окрестности i-го узла с помощью центральных разностей:
yi+1 − yi−1 |
= fi, i = 0, 1, ..., y0 = A. |
2h |
|
Полученная система уравнений не замкнута – необходимо доопределить y1, что можно сделать с помощью метода Эйлера
y1 = y0 + hf(x0, y0).
Можно показать, что построенная схема имеет второй порядок точности. Рассмотрим модифицированный метод Эйлера. Запишем систему урав-
нений |
− yi |
|
|
|
yi+1 |
= fi+1/2, i = 0, 1, ..., y0 = A, |
|
|
|
||
|
h |
||
|
|
||
где
h
fi+1/2 = f(xi+1/2, yi+1/2) = f(xi + 2 , yi+1/2),
а значение yi+1/2 вычисляем по методу Эйлера
yi+1/2 |
= yi + |
h |
fi = yi + |
h |
f(xi, yi), i = 0, 1, ... . |
|
|
||||
|
2 |
2 |
|
||
3
В результате получаем следующую разностную схему:
yi+1 − yi |
= f(xi + |
h |
, yi + |
h |
f(xi, yi)), i = 0, 1, ..., y0 = A, |
h |
|
|
|||
2 |
2 |
|
|||
которая, как можно показать, имеет второй порядок точности.
1.2. Метод Эйлера с пересчетом (предиктор – корректор)
Предиктор – предсказание результата, корректор – уточнение результата. Заменим правую часть нашего уравнения на среднее значений в соседних узлах:
yi+1 − yi |
= |
1 |
[f(xi, yi) + f(xi+1, yi+1)] . |
(4) |
|||
|
|
|
|||||
h |
2 |
||||||
|
|
|
|||||
Как видно, наклон интегральной кривой посередине отрезка [xi, xi+1] приближенно заменяется средним арифметическим наклонов на границах этого отрезка. Отсюда
yi+1 |
= yi + |
h |
[f(xi, yi) + f(xi+1, yi+1)] . |
(5) |
|
||||
|
2 |
|
|
|
Искомое значение yi+1 входит и в правую сторону тоже, и его нельзя в общем случае выразить явно. Но его можно найти по формуле метода Эйлера:
|
|
|
yi+1 = yi + hf(xi, yi) – предиктор |
(6) |
и, подставляя в правую часть, получаем |
|
|||
yi+1 |
= yi + |
h |
[f(xi, yi) + f(xi+1, (yi + hf(xi, yi)))] – |
корректор. (7) |
|
||||
|
2 |
|
|
|
1.3.Методы Рунге – Кутта
Рассмотренные ранее метод Эйлера и его модификация относятся к классу методов Рунге – Кутта. Суть методов – для вычисления yi+1 используется yi, а также значения функции f(x, y) в некоторых специальных точках. Широко распространен метод Рунге – Кутта четвертого
порядка. Алгоритм этого метода имеет вид |
|
|
|||
|
h |
|
|
||
yi+1 = yi + |
|
(k0 + 2k1 + 2k2 + k3), |
i = 0, 1, ... , |
||
6 |
|||||
k0 = f(xi, yi), k1 = f xi + 2 , yi + 20 |
, |
||||
|
|
|
h |
hk |
|
4
k2 |
= f |
xi + 2 , yi + |
21 |
, |
|
|
|
|
h |
hk |
|
k3 = f(xi + h, yi + hk2).
1.4.Многоточечные методы
Многоточечные, или многошаговые методы решения задачи Коши отличаются тем, что решение в текущем узле зависит от данных не только в одном предыдущем узле, но и в ряде предшествующих.
Такие методы можно строить, используя метод неопределенных коэффициентов. Для этого запишем разностное соотношение для i-го узла в виде
a0yi + a1yi−1 + ... + amyi−m |
= b |
f |
+ b |
f |
i−1 |
+ ... + b f |
i−l |
. |
(8) |
h |
0 |
i |
1 |
|
l |
|
|
Коэффициенты ai и bi подбираются так, чтобы получить аппроксимацию с нужным порядком точности. Если
a0 = −a1 = 1, a2 = ... = am = 0,
то получаем методы Адамса. Уравнение y0 = f(x, y) эквивалентно интегральному соотношению
xi+1 |
|
|
yi+1 − yi = Z |
f dx . |
(9) |
xi
Если решение в узлах вплоть до i-го уже вычислено, то по известным значениям fi = f(xi, yi) можно интерполировать подинтегральную функцию полиномами различной степени. Далее, вычисляя интеграл от выбранного интерполяционного полинома, получаем формулы Адамса.
Приближая функцию f на отрезке [xi, xi+1] полиномом в форме Ньютона
f = fi + fi − fi−1 (x − xi)+ h
fi − 2fi−1 + fi−2
+ 2h2 (x − xi)(x − xi−1) + ... (10)
и учитывая первые два слагаемых при вычислении интеграла, получаем метод Адамса второго порядка точности
yi+1 − yi |
= |
3 |
f |
i − |
|
1 |
f |
i−1 |
. |
(11) |
h |
|
|
||||||||
|
2 |
2 |
|
|||||||
5
Учитывая в (10) три слагаемых, приходим к методу третьего порядка точности
yi+1 − yi |
= |
|
1 |
(23f |
i − |
16f |
+ 5f |
) |
(12) |
|
h |
12 |
|||||||||
|
|
i−1 |
i−2 |
|
||||||
и случае четырех слагаемых в (10) метод Адамса четвертого порядка
yi+1 − yi |
= |
1 |
(55f |
i − |
59f |
i−1 |
+ 37f |
i−2 − |
9f |
i−3 |
). |
(13) |
|
h |
24 |
||||||||||||
|
|
|
|
|
|
Метод Адамса по сравнению с методом Рунге – Кутты требует меньших затрат при нахождении очередного значения yi+1, т.к. требуется найти лишь fi, а fi−1, ... ,fi−3 уже известны к этому моменту, по формулам Рунге – Кутты на каждом шаге надо находить четыре значения f.
Формулы Рунге – Кутты позволяют проводить вычисления с переменным шагом. При использовании многошаговых методов может неконтролируемо возрастать погрешность вычислений, т.е. возникать неустойчивость.
1.5.Задания для самостоятельной работы
Решить задачу Коши для обыкновенных дифференциальных уравнений
а) методом Эйлера; б) методом Эйлера с пересчетом;
в) методом Рунге – Кутты; г) методом Адамса.
1. |
y0 |
= 2x − 3 + y, y(0) = 0, a = 0, b = 1, h = 0, 1 |
|||||||||||
2. |
y0 |
= |
1 |
|
, y(−1) = 0, a = −1, b = 1, h = 0, 1 |
||||||||
|
|
||||||||||||
x + 2y |
|||||||||||||
3. |
y0 |
= |
2 − y |
, |
|
y(0) = 3, a = 0, b = 1, |
h = 0, 1 |
||||||
ctg x |
|||||||||||||
4. |
y0 |
|
2x4 + 2y |
|
|
|
|||||||
= |
|
|
|
, y(−1) = 2, a = −1, b = 1, h = 0, 2 |
|||||||||
x |
|
||||||||||||
5. |
y0 |
= |
4x + 2y |
|
, y(0) = 2, |
a = 0, |
b = 1, h = 0, 1 |
||||||
2x + 1 |
|||||||||||||
|
|
|
|
|
|
|
|
|
|||||
6. |
y0 |
= |
x + 2y |
, |
|
y(2) = 2, |
a = 2, |
b = 3, |
h = 0, 1 |
||||
|
|
||||||||||||
|
|
|
x |
|
|
|
|
|
|
||||
7. |
y0 |
= |
|
|
y |
|
|
, y(0) = 1, |
a = 0, |
b = 1, h = 0, 1 |
|||
|
|||||||||||||
2x + y − 4 ln x |
|||||||||||||
|
|
|
|
|
|
|
|
|
|
6 |
|
|
|
8. |
y0 |
= 2x(x2 + y), |
y(0) = −1, |
a = 0, b = 1, h = 0, 1 |
||||
9. |
y0 |
= |
2xy − y2 |
, |
y(2) = 4, a = 2, b = 3, |
h = 0, 1 |
||
|
|
|
x2 |
|
|
|
|
|
10. |
y0 |
= |
(x+y) ln (x+y)− ln x+y |
, y(1) = e |
− |
1, a = 1, b = 2, h = 0, 1 |
||
|
|
|
|
x |
|
|
||
1.6.Примеры процедур в среде Maple
1.6.1.Метод Эйлера
>restart;
>Euler:=proc(f,y_a,a,b,n)
# f функция
# y_a начальное условие задачи Коши
# a, b границы
# n количество точек разбиения промежутка [a,b] local i,h,x,y;
# i переменная циклов
# h шаг
# x узлы
# y приближения в решению задачи y[0]:=y_a;
h:=(b-a)/n;
for i from 0 to n-1 do
# задаем конечное множество точек на отрезке [a,b] x[i]:=a+h*i;
# вычисление приближенных значений y(x[i]) y[i+1]:=y[i]+h*f(x[i],y[i]);
end do;
for i from 0 to n-1 do print(x=x[i],y=y[i]);
end do; end proc:
>f:=(x,y)->y+(1+x)*yˆ 2:
>dy/dx=f(x,y);
Проверим работу процедуры > Euler(f,-1,1,1.5,5);
x = 1., y = -1
7
x = 1.100000000, y = -0.9000000000 x = 1.200000000, y = -0.8199000000 x = 1.300000000, y = -0.7539980778 x = 1.400000000, y = -0.6986398723
Сравним полученные результаты с решением данной задачи встроенными командами Maple
>proverka:=proc(n,a,b) local h,ODU,i,t; h:=abs(b-a)/n;
ODU:=di (y(x),x)-f(x,y(x)); for t from 0 to n do
i:=a+h*t; print(evalf(eval(dsolve(ODU),_C1=0,x=i)));
end do; end proc:
>proverka(5,1,1.5);
y(1.) = -1.
y(1.100000000) = -0.9090909091
y(1.200000000) = -0.8333333333
y(1.300000000) = -0.7692307692
y(1.400000000) = -0.7142857143
y(1.500000000) = -0.6666666667
1.6.2.Метод Эйлера с пересчетом
>restart;
>EulerPereschet:=proc(f,y_a,a,b,n)
# f функция
# y_a начальное условие задачи Коши
# a, b границы
# n количество точек разбиения промежутка [a,b] local i,h,x,y,y1;
# i переменная циклов
# h шаг
# x узлы
# y приближенное решение задачи
# y1 скорректированное решение задачи y1[0]:=y_a;
h:=(b-a)/n;
8
for i from 0 to n do
#задаем конечное множество точек на отрезке [a,b] x[i]:=a+h*i;
#вычисление первого приближения (предиктора) y1[i+1]:=y1[i]+h*f(x[i],y1[i]);
end do; y[0]:=y1[0];
for i from 0 to n do
# уточнение результата (корректор) y[i+1]:=y[i]+(h/2)*(f(x[i],y[i])+f(x[i+1],(y1[i+1]))); print(x=x[i], y=y[i]);
end do; end proc:
>f:=(x,y)->y+(1+x)*yˆ 2:
>dy/dx=f(x,y);
Проверим работу процедуры
>EulerPereschet(f,-1,1,1.5,5);
x = 1., y = -1
x = 1.100000000, y = -0.9099500000 x = 1.200000000, y = -0.8355555936 x = 1.300000000, y = -0.7728574240 x = 1.400000000, y = -0.7191700795 x = 1.500000000, y = -0.6725981326
1.6.3.Метод Рунге – Кутты
>restart;
>RungeKutt:=proc(f,y_a,a,b,n)
# f функция
# y_a начальное условие задачи Коши
# a, b границы
# n количество точек разбиения промежутка [a,b] local i,h,x,y,dy,K1,K2,K3,K4;
y[0]:=y_a; h:=(b-a)/n;
for i from 0 to n do
# задаем конечное множество точек на отрезке [a,b] x[i]:=a+h*i;
end do;
9
for i from 0 to n do K1[i]:=h*f(x[i],y[i]); K2[i]:=h*f(x[i]+h/2,y[i]+K1[i]/2); K3[i]:=h*f(x[i]+h/2,y[i]+K2[i]/2); K4[i]:=h*f(x[i]+h,y[i]+K3[i]);
dy[i]:=(1/6)*(K1[i]+2*K2[i]+2*K3[i]+K4[i]);
y[i+1]:=y[i]+dy[i]; end do;
for i from 0 to n do print(x=x[i], y=y[i]);
end do; end proc:
>f:=(x,y)->y+(1+x)*yˆ 2:
>dy/dx=f(x,y);
Проверим работу процедуры > RungeKutt(f,-1,1,1.5,5);
x = 1., y = -1
x = 1.100000000, y = -0.9090933148 x = 1.200000000, y = -0.8333367499 x = 1.300000000, y = -0.7692344925 x = 1.400000000, y = -0.7142893912 x = 1.500000000, y = -0.6666701276
1.6.4.Метод Адамса
>restart;
>Adams:=proc(f,y_a,a,b,n)
# f функция
# y_a начальное условие задачи Коши
# a, b границы
# n количество точек разбиения промежутка [a,b] local i,h,x,y,dy,K1,K2,K3,K4;
y[0]:=y_a; h:=(b-a)/n;
for i from 0 to n do
# задаем конечное множество точек на отрезке [a,b] x[i]:=a+h*i;
end do; y[1]:=y[0]+h*f(x[0],y[0]);
10