Численное решение ОДУ первого порядка: ряды Тейлора и метод Рунге-Кутты / Хабр

Постановка задачи

Если не всем, то многим, кто обучался на направлениях из УГС 01/02/03, доводилось в курсах дисциплин “высшая математика” или “дифференциальные уравнения” видеть обыкновенное дифференциальное уравнение первого порядка вида

$y' = f(x,y),\ \ y=y(x),\ \ y(x_0)=y_0,$

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

Метод Эйлера;
Метод Рунге—Кутты;
Метод степенных рядов;
Метод аппроксимирующих многочленов;
Метод Галёркина;

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

Ряды Тейлора и остатки

Хитрые трюки

Если рассмотреть задачу

$y' = f(x,y),\ \ y=y(x),\ \ y(x_0)=y_0,$

то можно заметить, что искомую функцию можно искать в виде ряда Тейлора

$y(x) = y_0 + \sum\limits_{n=1}^{+\infty} \frac{(x-x_0)^n}{n!} y^{(n)}_0 ,\ \ y^{(n)}_0 = y^{(n)}(x_0) ,$

где, как это является очевидным,

$y' = f(x,y);\\ y'' = \frac{d}{dx} y(x);\\ \ldots\\y^{(q)} = \frac{d}{dx} y^{(q-1)}.$

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

$y_{i+1}=y_i+h_x f(x_i , y_i) ,$

и метод Рунге-Кутты, при котором из соответствующего ряда Тейлора берутся только первые несколько членов этого ряда с последующей заменой производных делёнными разностями и достаточно хитрой перестановкой слагаемых и строятся давно прописанные в разных книгах/руководствах/методичках формулы(к примеру, четвёртого порядка):

$\textbf{y}_{n+1} = \textbf{y}_n + {h \over 6}(\textbf{k}_1 + 2\textbf{k}_2 + 2\textbf{k}_3 + \textbf{k}_4), \\\textbf{k}_1 = \textbf{f} \left( x_n, \textbf{y}_n \right) ,\\\textbf{k}_2 = \textbf{f} \left( x_n + {h \over 2}, \textbf{y}_n + {h \over 2} \textbf{k}_1 \right), \\\textbf{k}_3 = \textbf{f} \left( x_n + {h \over 2}, \textbf{y}_n + {h \over 2} \textbf{k}_2 \right) , \\ \textbf{k}_4 = \textbf{f} \left( x_n + h, \textbf{y}_n + h\ \textbf{k}_3 \right) .$

От приближённого-приближённого к приближённому

Немного философии

В методе Рунге—Кутты, как ранее была указано, используются не точные значения производных, а приближённые(получаемые путём применения делённых разностей). Таким образом погрешность, накладываемая тем, что берётся лишь некоторая часть от ряда Тейлора, совмещается с тем, что берутся приближённые значения производных.

Производные

Из уравнения

становится очевидно, что

$y'' = \frac{d}{dx} y' = \frac{d}{dx} f(x,y(x)) = \frac{}{dx}\left[ \frac{\partial f}{\partial x} dx + \frac{\partial f}{\partial y} dy \right] = \frac{}{dx}\left[ \frac{\partial f}{\partial x} dx + \frac{\partial f}{\partial y} y' dx \right] = f'_x + f'_y f ,$

а тогда ещё и выходит так:

$y''' = \frac{d}{dx}y'' = \frac{d}{dx}\left[ f'_x + f'_y f \right]=\\= \frac{}{dx}\left[ \underbrace{\left\{ \frac{\partial f'_x}{dx}dx + \frac{\partial f'_x}{dy}fdx \right\}}_{d f'_x} + \underbrace{\left \{ \overbrace{\left( \frac{ \partial f'_y}{\partial x} + \frac{\partial f'_y}{\partial y} f \right) dx}^{d(f'_y)} f + f'_y \overbrace{\left( \frac{\partial f}{\partial x} + \frac{\partial f}{\partial y} f \right)dx}^{d(f)} \right\} }_{d(f'_y f)} \right]=\\=f''_{xx}+f''_{xy}f+f''_{yx}f+f''_{yy}f^2 + f'_y f'_x +(f'_y)^2 f$

Таким образом, можно искать приближённое решение через приближение

$y\approx y_0 + (x-x_0) y'_0+\frac{(x-x_0)^2}{2}y''_0+\frac{(x-x_0)^3}{6}y'''_0$$display.$$$

Если разбить исходный отрезок, на котором ищется решение на равные отрезки(каждый из которых равен ), то выходит такая формула:

${ y_i = y_{i-1} + h y'_{i-1}+\frac{h^2}{2}y''_{i-1}+\frac{h^3}{6}y'''_{i-1} , }$

где ${ y_i }$ это приближённое значение искомой функции в точке с номером а ${ y'_{i-1}, y''_{i-1}, y'''_{i-1} }$ значения первой, второй и третьей производной искомой функции в точке ${(x_i , y_i) ,}$ где ${x_i = h\cdot i + x_0 . }$

Код программы с реализацией

Ниже представлен код с реализацией описанного метода:

from matplotlib import pyplot as plt
import math

y_true = lambda x: 2* math.exp(math.sin(x))

f = lambda x,y: math.cos(x)*y

f_x = lambda x,y: -math.sin(x)*y

f_y = lambda x,y: math.cos(x)

f_xx = lambda x,y: -f(x,y)

f_yy = lambda x,y:0

f_xy = f_yx = lambda x,y:-math.sin(x)

x0 = -1
x1 = 10

y0 = y_true(x0)

x=[x0]
y = [y0]

n = 100

dx = (x1-x0)/(n-1.0)

while len(x) < n:
    y_i, x_i = y[-1], x[-1]
    S = f_xx(x_i,y_i)+f_xy(x_i,y_i)*f(x_i,y_i)+f_yx(x_i,y_i)*f(x_i,y_i)
    + f_yy(x_i,y_i)*(f(x_i,y_i))**2 + f_x(x_i,y_i)*f_y(x_i,y_i)+(f_y(x_i,y_i))**2 * f(x_i,y_i)
    S = S*dx/6 + (f_x(x_i,y_i)+f_y(x_i,y_i)*f(x_i,y_i))/2
    S = S*dx + f(x_i,y_i)
    S = S*dx + y_i
    x.append(x_i+dx)
    y.append(S)
plt.plot(x,[y_true(xi) for xi in x], color="blue")
plt.plot(x,y,color="red",linestyle="dashed")

    
plt.show()

А какие могут быть проблемы?

Проблемы с данным методом могут быть на том основании, что ряд Тейлора в общем случае не сходится поточечно к функции. Такие примеры существуют( к примеру, пример Коши, о котором рассказывают в курсе матанализа на профильных направлениях ВУЗов). Таким образом встаёт вопрос: а как быть? Автор данной статьи предлагает сделать иной шаг(если таковой вообще возможен): брать не ряд Тейлора, где идёт разложение в степенной ряд с коэффициентами, состоящими из дроби, где в знаменателе факториал, а в числителе значение производной соответствующего порядка в определённой точке, а в ряд с многочленами Чебышева, раскладывая функцию по многочленам Чебышева.

Одно интересное применение метода Рунге—Кутты

Если рассмотреть задачу

${y'=f(x),\ \ y(0)=0 ,}$

то можно заметить, что решение данного ОДУ даёт такую функцию, что

${y(x)=\int_{0}^{x} f(x)dx ,}$

а тогда выходит, что применение метода Рунге—Кутты находит своё место в численном интегрировании.

Заключение и дальнейшие планы

Как уже удалось понять, уважаемый читатель данной статьи, целью данной статьи был разбор основания метода Рунге—Кутты. Данный метод позволяет решать приближённо отдельные задачи.

В дальнейшем мной планируется применение данного метода к “задаче n тел”.

Автор данной статьи выражает надежду, что данная статья будет полезна(хотя бы студентам/обучающимся первых-вторых курсов разных университетов/институтов/академий).

Только зарегистрированные пользователи могут участвовать в опросе. Войдите, пожалуйста.

40%Да2

60%Нет3

Проголосовали 5 пользователей. Воздержавшихся нет.