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

Дифференциальными называются функциональные уравнения (т.е. искомым является функция), содержащие кроме искомой функции также ее производные. Рассмотрим явное дифференциальное уравнение первого порядка, имеющее общий вид: $$y'(t)=f(t,y).$$ Например, таковым явлется уравнение $y'(t)=y(t)$, явное решение которого известно: $y(t)=c\cdot e^t$, где $c$ - константа интегрирования. Уже из примера видно, что чтобы получить единственное решение дифференциального уравнения первого порядка, нужно зафиксировать начальное условие, т.е. зафиксировать значение $y(0)$. Для нетривиальных диф. уравнений аналитическое решение в элементарных функциях обычно не существует. Однако решение можно найти приближенно, воспользовавшись методами численного интегрирования. Рассмотрим простейший метод численного интегрирования - метод Эйлера.

Решения дифференциального уравнения являются дифференциируемыми функциями, поэтому для них справедливо следующее представление в малой окрестности каждой точки $t$: $$y(t+\delta)=y(t)+\delta y'(t)+o(\delta),\quad \delta\to0.$$ Отбрасывая бесконечно малое, можно получить приближенную формулу, по которой можно находить значение в "следующий" момент времени, зная значение в "предыдущий" момент времени. При численном решении мы не можем найти решение во все моменты времени, вместо этого мы вычисляем значения решения в дискретные моменты времени $t_k$, которые мы считаем упорядоченными $t_k<t_{k+1}$. В простейшем случае можно считать приращении времени постоянным, тогда $t_k=t_0+\Delta\cdot k$, где $\Delta$ - шаг времени, а $t_0$ - начальный момент времени.

Зная значение функции и производной в момент времени $t_k$, можно найти приближенное значение решения в момент времени $t_{k+1}$: $$y_a(t_{k+1})=y_a(t_k)+(t_{k+1}-t_k)y_a'(t_k).$$ Подставляя значение производной из дифференциального уравнения, получаем формулу нахождения решения по методу Эйлера: $$y_a(t_{k+1})=y_a(t_k)+(t_{k+1}-t_k)f(t_k,y_a(t_k)).$$ Найденное таким образом решение $y_a(t)$ не совпадает с точным решением $y(t)$ уравнения. Если выбрать значение $y_a(t)$ совпадающим с $y(t)$, то величина $h(\delta)=y_a(t+\delta)-y(t+\delta)$ характеризует точность одного шага данного метода интегрирования. Точность методы выражаеют порядком метода, который выражает скорость уменьшения погрешности с уменьшением шага интегрирования. Так на одном шаге интегрирования метод имеет порядок $\alpha$, если $h(\delta)\sim \delta^\alpha$, т.е. $h(\delta)$ пропорционально степени $\alpha$ шага интегрирования $\alpha$, Метод Эйлера имеет второй порядок точности одного шага интегрирования, если решение достаточно гладкое.

Реализуем метод Эйлера. Ограничимся только автономными системами, т.е. уравнениями вида $$y'(t)=f(y(t)),$$ где функция $f$ явным образом не зависит от времени, а зависит только от $y$.

На практике величина шага интегрирования $\Delta$ является вспомогательной величиной, изменение которой позволяет контролировать точность нахождения результата. Интересуют же нас обычно значения функции на неком интервале $t\in[t_0,t_0+T]$. Если мы используем равномерные приращения по времени $\Delta$, т.е. $t_k=t_0+k\cdot\Delta$, то для достижения правого конца $t_0+T$ интервала потребуется $N$ шагов длины $T/N$. Каждый шаг интегрирования добавляет погрешность, а количество шагов растет с уменьшением длины шага, поэтому можно ожидать, что точность нахождения численным методом решения уравнения на интервале отличается от точности одного шага интегрирования. Пусть $y_a$ приближенное решение уравнения, полученное интегрированием с шагом $\delta$, а $y$ - точное решение уравнения, тогда величину $H(\delta)=y_a(t_0+T)-y(t_0+T)$ будем называть погрешностью интегрирования на интервале. Вычислим эту погрешность для метода Эйлера.

Как мы видим, на интервале метод Эйлера имеет первый порядок точности. Существует однако ситуация, когда метод Эйлера дает точный ответ: это происходит, когда производная постоянна.

При приращении аргумента на некий шаг метод Эйлера приращает значение решения таким оразом, что производная остается постоянной и равной значению производной в начале шага. Однако это вообще говоря очень грубое допущение. Чтобы получить более точный метод, мы можем ввести поправки на изменении производной. Раскладывая решение по формуле Тейлора, получаем при малых $t$: $$y(t)=y(0)+y'(0)t+y''(0)\frac{t^2}2+o(t^2).$$ Так $y'(t)=f(t,y(t))$, то дифференцируя по $t$, получаем: $$y''(t)=\frac{\partial}{\partial t}f(t,y(t))+y'(t)\frac{\partial}{\partial y}f(t,y(t)).$$ Подставляя производные в разложение, получаем: $$y(\delta)=y(0)+f(0,y(0))\delta+\left(\frac{\partial}{\partial t}f(0,y(0)) +f(0,y(0))\frac{\partial}{\partial y}f(0,y(0))\right)\frac{\delta^2}2+o(\delta^2).$$ Эта формула обновления позволяет получить решение в точке $t=\delta$, зная решение в точке $t=0$, с погрешностью меньшей $\delta^2$. Это означает, что уменьшая шаг интегрирования в два раза, мы получаем точность выше более чем в 4 раза. Тот же результат можно было получить, используя формулу Эйлера с шагом $\delta/4$.

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

Модифицированный метод Эйлера вычисляет значение решения $y(\delta)$ по данному значению $y(0)$ с помощью следующей формулы: $$y(\delta)=y(0)+f\left(\frac\delta2,y_{\frac\delta2}\right)\delta,$$ где $y_{\frac\delta2}$ - значение решение на половинном шаге $t=\delta/2$: $$y_{\frac\delta2}=y(0)+f(0,y(0))\frac\delta2.$$ Т.е. значение производной на всем шаге считается постоянным, однако для более точной оценки производной сначала находится приближенное значение решения на половине шага интегрирования, которое и используется для вычисления производной из дифферецниального уравнения.

Чтобы получить точность выше, мы должны найти больше промежуточных значений внутри одного шага интегрирования и использовать эти значения для оценки производных высоких порядков. Полученные таким образом методы относятся к семейству методов Рунге-Кутты. Самой известным методом из этого классая является классический метод Рунге-Кутты четвертого порядка, использующий для вычислений четыре спомогательные точки. Соответствующая формула интегрирования для решения задачи $$y'(t)=f(t,y(t)),\quad y(0)=y_0,$$ имеет вид: $$y(\delta)=y(0)+\frac\delta6(k_1+2k_2+2k_3+k_4),$$ $$k_1=f(0,y(0)),$$ $$k_2=f\left(\frac\delta2,y(0)+\frac\delta2k_1\right),$$ $$k_3=f\left(\frac\delta2,y(0)+\frac\delta2k_2\right),$$ $$k_4=f\left(\delta,y(0)+\delta k_3\right).$$ Классический метод Рунге-Кутты удобен тем, что промежуточные шаги интегрирования делаются последовательно, причем в памяти достаточно хранить только значение одного промежуточного шага. Существуют однако более точные и более устойчивые варианты методы Рунге-Кутты.

Рассмотренные выше методы относились к семейству явных методов, т.е. позволяли сразу вычислить значения решения по вышеприведенным формулам, так как все величины в формулах были известны заранее. Однако формулы интегрирования могут содержать искомое решение и в левой и правой частях, тогда метод называют неявным. Например, формулу интегрирования для методы Эйлера: $$y(\delta)=y(0)+f(0,y(0))\delta,$$ можно было переписать так $$y(\delta)=y(0)+f(\delta,y(\delta))\delta,$$ т.е. для приближения производной решения на шаге интегрирования используется точное значение производной в конце, а не в начале шага. Для нахождения нового значения решения $y(\delta)$ в неявной схеме требуется решить уравнение, поэтому неявные методы вычислительно более сложны, чем явные. Однако неявные методы могут обладать привлекательными свойствами, такими как устойчивость и сохранение первых интегралов, что объясняет популяроность и важность этих методов.

Мы рализуем неявный метод Эйлера, используя для решения уравнения метод Ньютона.

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

Задания

1. Объясните, что понимается под порядком метода интегрирования? Как порядок метода можно определить на графике зависимости точности решения от шага интегрирования?
2. Вычислите погрешности интегрирования на интервале для вышеприведенных методов. Какой порядок имеют рассмотренные выше методы? Как связаны ошибки интегрирования на одном шаге и на интервале?
3. Объясните, почему неявный метод Эйлера дает более точное решение, чем явный метод, но различие заметно только на больших шагах интегрирования.
4. Решите численно на интервале $t\in[0,1]$ уравнение $$y'(t)=cos(y(t)),\quad y(0)=1.$$ Обоснуйте выбор метода интегрирования и шага интегрирования.
5. Решите уравнение $$y''(t)=-y(t),\quad y(0)=1,\quad y'(0)=0.$$ Это уравнение можно свести к системе уравнений первого порядка, с помощью новой переменной $z(t)=y'(t)$: $$\begin{cases} z'(t)=-y(t),\quad z(0)=0,\\ y'(t)=z(t),\quad y(0)=1.\\ \end{cases}$$ Тогда для решения этой системы можно использовать рассмотренные выше методы, если учесть, что искомое решение будет задаваться вектором: $$Y(t)=\begin{pmatrix}y(t)\\z(t)\end{pmatrix}.$$