Решение ОДУ с помощью численных методов
Anastasia Makarevich
ПОСТАНОВКА ЗАДАЧИ
Задано обыкновенное дифференциальное уравнение первого порядка: \[ \frac{dy}{dx} = cos(x e ^x) + y x \] И начальное значение: \(y(0)=1\) (то есть нам предстоит решить задачу Коши). Решить уравенение с помощью методов Рунге-Кутта, в класс которых входят:
- Исправленные метод Эйлера (метод Рунге-Кутта 2-го порядка)
- Модифицированный метод Эйлера
- Метод Рунге-Кутты
Сравнить методы, затем сравнить результаты для разного количества точек и разного размера шага:
1) 10 точек с шагом h = 0.2
2) 20 точек с шагом h = 0.1
ХОД РАБОТЫ
Общая идея метода Эйлера
Суть метода Эйлера в самой примитивной его форме сводится к тому, что мы берем начальную точку, вычисляем в ней значение производной - находим таким образом тангенс угла наклона, затем находим уравнение прямой с таким тангенсом угла наклона и проходящей через нашу начальную точку. делаем небольшой шаг h, получая таким образом новый х. Используя уравнение прямой: вычисляем новый y. Затем подставляем уже эти х и у, чтобы найти новый тангенс угла наклона и так далее, пока не достигнем конца интересующего нас интервала. В итоге, очевидно, получается ломаная линия, причем чем меньше шаг, тем точнее будет приближение и тем больше она будет походить на реальную функцию. Однако данный метод в чистом виде используется редко. Эйлер сам же его доработал, представив улучшенную версию своего метода.
Исправленный (усовершенствованный) метод Эйлера
Чтобы повысить точность решения, Эйлером было предложено брать тангенс угла наклона не в одной точке, а в двух и затем считать среднее значение. ТО есть сначала мы считаем производную в точке \((x_i, y_i)\),и получаем тангенс угла наклона первой касательной, а затем в \((x_{i+1}, y^{(*)}_{i+1}) = (x_i + h, y_i + h y_i')\), и получаем тангенс угла наклона второй касательной. А потом считаем среднее арифметическое. Этот метод имеет степень точности \(O(h^2)\).
В качестве отправной точки в нашем случае используется начальное значение \((x_0, y_0) = (0,1)\), а в качестве второй точки \((x_{1}, y^{(*)}_{1}) = (x_0 + h, y_0 + h y_0') = (0.2, 1.2)\) \[ y'(x_0, y_0) = cos(x_0 e ^{x_0}) + y_0 x_0 = cos(0 e^{0}) + 1*0 = 1\] \[ y'(x_1, y^{(*)}_1) = cos(x_1 e ^{x_1}) + y_1 x_1 = cos(0.2 e^{0.2}) + 1.2*0.2 = 1.2399909\] Затем мы находим среднее между этими значениями по формуле \(\frac{1 + 1.2399909}{2} = 1.1199955\) и получаем тангенс угла наклона новой прямой - \(\overline{L}\), проходящей через точку \((x_1, y_1)\). В общем виде эта формула записывается так: \[\Phi(x_i,y_i,h) = \frac{1}{2} [f(x_i, y_i) + f(x_i + h, y_i + h y_i')]\] После этого мы проводим через точку \((x_i, y_i)\) прямую L, параллельную \(\overline{L}\). Уравнение этой прямой имеет вид: \[y = y_i + (x-x_i)\Phi(x_i,y_i,h)\] В нашем конкретном случае уравнение этой прямой будет: \[y = 1 + (x-0)*1.1199955\] Получив это уравнение, легко можно найти следующую точку- \((x_1, y_1)\): \[x_1 = x_0+ h = 0 + 0.2 = 0.2\] \[y_1 = 1 + (x_1-0)*1.1199955 = 1 + h 1.1199955 = 1 + 0.2*1.1199955 = 1.2239991\] Затем операция повторяется сначала, но в качестве отправной точки используется только что найденная \((x_1, y_1)\). Напишем универсальную функцию, которая будет считать следующее значение у.
# Original derivative function
dydx <- function(x,y) {
cos((x*exp(x))*pi/180) + y*x
}
#Function that calculates the next y for (x+h)
#Arguments: xi, yi - current point, h - step, f - derivative
euler1 <- function(xi,yi,h,f) {
#Calcuate derivative at (xi, yi):
d1 = f(xi,yi)
#Calculate the "next"" point
xiplus1 = xi + h
yiplus1 = yi + h*f(xi,yi)
#Calculate derivative at the second point:
d2 = f(xiplus1, yiplus1)
tgavg = (d1+d2)/2
#Calculate the next point:
newy = yi + h*tgavg
newy
}Теперь напишем универсальную функцию, которая будет считать значения у на заданном интервале заданным методом:
library(knitr)
generateValues <- function(x0, y0, h, n, f, method) {
#Generate the array of x-values
ubound <- x0+h*n
xvals <- seq(x0,ubound,h)
# Initialize the array of y-vals with the initial value already known
yvals <- vector()
y = y0
for (x in xvals) {
yvals = c(yvals,y)
y = round(method(x,y,h,f),6)
}
cbind(xvals,yvals)
}
xy1 = as.data.frame(generateValues(0,1,0.2,10, dydx,euler1))
xy2 = as.data.frame(generateValues(0,1,0.1,20, dydx, euler1))
kable(generateValues(0,1,0.2,10, dydx,euler1))| xvals | yvals |
|---|---|
| 0.0 | 1.000000 |
| 0.2 | 1.223999 |
| 0.4 | 1.507391 |
| 0.6 | 1.877341 |
| 0.8 | 2.374122 |
| 1.0 | 3.059279 |
| 1.2 | 4.029362 |
| 1.4 | 5.439583 |
| 1.6 | 7.545549 |
| 1.8 | 10.778560 |
| 2.0 | 15.884640 |
kable(generateValues(0,1,0.1,20, dydx, euler1))| xvals | yvals |
|---|---|
| 0.0 | 1.000000 |
| 0.1 | 1.105500 |
| 0.2 | 1.223193 |
| 0.3 | 1.355638 |
| 0.4 | 1.505895 |
| 0.5 | 1.677658 |
| 0.6 | 1.875431 |
| 0.7 | 2.104748 |
| 0.8 | 2.372457 |
| 0.9 | 2.687093 |
| 1.0 | 3.059361 |
| 1.1 | 3.502775 |
| 1.2 | 4.034498 |
| 1.3 | 4.676470 |
| 1.4 | 5.456907 |
| 1.5 | 6.412330 |
| 1.6 | 7.590313 |
| 1.7 | 9.053226 |
| 1.8 | 10.883378 |
| 1.9 | 13.190113 |
| 2.0 | 16.119675 |
Построим получившиеся ломаные линии:
library(ggplot2)
ggplot(xy1, aes(x = xy1[,1], y = xy1[,2])) +
geom_line() +
geom_line(data = xy2, aes(x = xy2[,1], y = xy2[,2]),color="chartreuse4", linetype=2)+
labs(title = "Euler2", y = "y", x="x")
Как видим, они практически совпадают, но при меньшем шаге функция получается более “гладкой”.
Модифицированный метод Эйлера
Следующий метод тоже является одним из методов Рунге-Кутты, и его идея похожа на идею исправленного метода, только в этот раз мы усредняем не тангенс угла наклона, а точки. Сначала мы как обычно находим производную в известной начальной точке: \[ y'(x_0, y_0) = cos(x_0 e ^{x_0}) + y_0 x_0 = cos(0 e^{0}) + 1*0 = 1\] Затем делаем шажок на \(\frac{h}{2}\) и находим тангенс угла наклона в новой точке:
\((x_0 + \frac{h}{2}, y_0 + \frac{h}{2} y'(x_0, y_0)))\)
x0 = 0
y0 = 1
h = 0.2
x2 = x0 + h/2
y2 = y0+(h/2)*d0
d2 = dydx(x2,y2)
Phi <- function(x0,y0,h) {
x2 = x0 + h/2
y2 = y0+(h/2)*d0
dydx(x2,y2)
}
newy <- function(x,y,h) {
y + h*Phi(x,y,h)
}euler2 <- function(xi,yi,h,f) {
#Calcuate derivative at (xi, yi):
d1 = f(xi,yi)
#Calculate the "next"" point
xiplus1 = xi + h/2
yiplus1 = yi +(h/2)*d1
#Calculate derivative at the second point:
d2 = f(xiplus1, yiplus1)
#Calculate the next point:
newy = yi + h*d2
newy
}\[ y'(x^{(*)}_1, y^{(*)}_1) = cos((x_0 + \frac{h}{2})(e ^{x_1+\frac{h}{2}}) + (y_0 + \frac{h}{2} y'(x_0, y_0)) (x + \frac{h}{2}) = cos(0.1 e^{0.1}) + 1.1*0.1 = 1.1099981\]
“Проведем”" через эту новую точку \((0.1, 1.1)\) прямую \(L^{(*)}\). Затем через предыдущую точку проведем прямую, параллельную \(L^{(*)}\) - \(L_0\). Уравнение этой прямой будет выглядеть следующим образом: \[y = y_i + (x - x_i) \Phi(x_i, y_i, h)\], где \(\Phi(x_i,y_i,h) = f(x_i + \frac{h}{2}, y_i+\frac{h}{2} y_m')\) Или, в данном конкретном случае для первого шага: \[y = 1 + (x - 0) 1.1099981\] Таким образом, новое значение у: \[y_1 = 1 + (x_1 - 0) 1.1099981 = 1 + (0.2 - 0) 1.1099981 = 1.2219996\]
xy3 = as.data.frame(generateValues(0,1,0.2,10, dydx,euler2))
xy4 = as.data.frame(generateValues(0,1,0.1,20, dydx,euler2))
kable(generateValues(0,1,0.2,10, dydx,euler2))| xvals | yvals |
|---|---|
| 0.0 | 1.000000 |
| 0.2 | 1.222000 |
| 0.4 | 1.502781 |
| 0.6 | 1.869049 |
| 0.8 | 2.360353 |
| 1.0 | 3.037048 |
| 1.2 | 3.993656 |
| 1.4 | 5.381853 |
| 1.6 | 7.450925 |
| 1.8 | 10.620613 |
| 2.0 | 15.615316 |
kable(generateValues(0,1,0.1,20, dydx, euler2))| xvals | yvals |
|---|---|
| 0.0 | 1.000000 |
| 0.1 | 1.105250 |
| 0.2 | 1.222661 |
| 0.3 | 1.354782 |
| 0.4 | 1.504657 |
| 0.5 | 1.675963 |
| 0.6 | 1.873181 |
| 0.7 | 2.101816 |
| 0.8 | 2.368680 |
| 0.9 | 2.682259 |
| 1.0 | 3.053195 |
| 1.1 | 3.494917 |
| 1.2 | 4.024477 |
| 1.3 | 4.663665 |
| 1.4 | 5.440497 |
| 1.5 | 6.391223 |
| 1.6 | 7.563046 |
| 1.7 | 9.017830 |
| 1.8 | 10.837185 |
| 1.9 | 13.129486 |
| 2.0 | 16.039623 |
Сравним графики для этого метода для разного шага:
ggplot(xy3, aes(x = xy3[,1], y = xy3[,2])) +
geom_line() +
geom_line(data = xy4, aes(x = xy4[,1], y = xy4[,2]),color="chartreuse4", linetype=2)+
labs(title = "Euler 2", y = "y", x="x")
Метод Рунге-Кутта четвертого порядка
В основе этого метода - усовершенствованные методы Эйлера. Он также является одним из наиболее употребимых при решении ОДУ. Метод Рунге-Кутты описывается следующей системой соотношений: \[y_{i+1} = y_i + \frac{h}{6} \Phi(k_1+2 k_2+ 2 k_3 + k_4)\] \[k_1 = f(x_i,y_i)\] \[k_2 = f(x_i + \frac{h}{2},y_i+\frac{h}{2} k_1)\] \[k_3 = f(x_i + \frac{h}{2},y_i+\frac{h}{2} k_2)\] \[k_4 = f(x_i + h,y_i + h k_3)\]
rungeKutta <- function(xi,yi,h,f) {
k1 = dydx(xi,yi)
k2 = dydx(xi + h/2, yi + (h*k1)/2)
k3 = dydx(xi + h/2, yi + (h*k2)/2)
k4 = dydx(xi + h, yi + h*k3)
yi + (h/6)*(k1+2*k2+2*k3+k4)
}
xy5 = as.data.frame(generateValues(0,1,0.2,10, dydx,rungeKutta))
xy6 = as.data.frame(generateValues(0,1,0.1,20, dydx,rungeKutta))
kable(generateValues(0,1,0.2,10, dydx,rungeKutta))| xvals | yvals |
|---|---|
| 0.0 | 1.000000 |
| 0.2 | 1.222888 |
| 0.4 | 1.505311 |
| 0.6 | 1.874646 |
| 0.8 | 2.371681 |
| 1.0 | 3.059107 |
| 1.2 | 4.035940 |
| 1.4 | 5.462654 |
| 1.6 | 7.606084 |
| 1.8 | 10.921673 |
| 2.0 | 16.207941 |
kable(generateValues(0,1,0.1,20, dydx, rungeKutta))| xvals | yvals |
|---|---|
| 0.0 | 1.000000 |
| 0.1 | 1.105346 |
| 0.2 | 1.222888 |
| 0.3 | 1.355189 |
| 0.4 | 1.505312 |
| 0.5 | 1.676960 |
| 0.6 | 1.874649 |
| 0.7 | 2.103932 |
| 0.8 | 2.371687 |
| 0.9 | 2.686493 |
| 1.0 | 3.059120 |
| 1.1 | 3.503178 |
| 1.2 | 4.035973 |
| 1.3 | 4.679654 |
| 1.4 | 5.462745 |
| 1.5 | 6.422223 |
| 1.6 | 7.606333 |
| 1.7 | 9.078441 |
| 1.8 | 10.922336 |
| 1.9 | 13.249567 |
| 2.0 | 16.209676 |
Сравним графики для этого метода для разного шага:
ggplot(xy5, aes(x = xy5[,1], y = xy5[,2])) +
geom_line() +
geom_line(data = xy6, aes(x = xy6[,1], y = xy6[,2]),color="chartreuse4", linetype=2)+
labs(title = "Runge-Kutta", y = "y", x="x")
Ошибка при вычислении методом Рунге-Кутта составляет \(e_T = О(h^5)\), то есть имеем четвертый порядок точности. При этом существенным недостатком метода является то, что функция вычисляется 4 раза для каждой точки. Другой недостаток заключается в том, что нет простых способов оценить ошибку. Впрочем, это верно для всех методов класса Рунге-Кутта. Важный момент при применении метода Рунге-Кутта - выбор шага интегрирования. Один из методов, позволяющих уменьшить ошибку выглядит следующим образом. На каждом шаге мы оцениваем соотношение:
\[ \theta = \left| \frac{k_2 - k_3}{k_1-k_2} \right| \] И если \[ \theta > 0.01 \], тогда шаг следует уменьшить.
Теперь, когда все три метода реализованы, сравним графики всех методов по группам: для каждого размера шага.
Для шага h = 0.2 графики ломаных для всех трех методов будут выглядеть так: 
Для шага h = 0.1 графики ломаных для всех трех методов будут выглядеть так:
И, наконец, сравним все полученные решения и убедимся, что в целом мы получили одно решение с разной степенью точности.
Выводы
Все перечисленные методы относятся к классу одноступенчатых методов, то есть для вычисления следующего значения у достаточно информации лишь об одной точке. Большим плюсом этих методов является простота для понимания, однако с повышением точности повышается количество операций, то есть повышается вычислительная сложность. Существенным недостатком методов данного класса является, как уже было сказано, то, что не существует простых способо оценить ошибку. Среди других недостатков можно выделить также то, что, например, метод Рунге-Кутта иногда может дать искаженный результат при неблагоприятных условиях (например, плохо подобран шаг или же для данной области нельзя получить корректное решение).