1. Задачи с начальными условиями

Дифференциальные уравнения описывают соотношения между неизвестной функцией и ее производными. Решить дифференциальное уравнение означает: найти функцию, удовлетворяющую этим соотношениям и некоторым дополнительным условиям. В этом курсе мы, в основном, будем иметь дело с одним классом задач (и, соответственно, с одним видом дополнительных условий) — задачами с начальными условиями1.

В самом общем виде поведение системы в задаче с начальными условиями описывается обыкновенным дифференциальным уравнением (ОДУ) в форме

$$ \dot{\textbf{x}} = f(\textbf{x},t), $$

где \(f\) — известная функция (т.е. нечто, что мы можем вычислить, зная \(\textbf{x}\) и \(t\)), \(\textbf{x}\)состояние системы, \(\dot{\textbf{x}}\) — производная \(\textbf{x}\) по времени. Как правило, \(\textbf{x}\) и \(\dot{\textbf{x}}\) являются векторами. Как подсказывает само название, в задачах с начальными условиями мы располагаем значением \(\textbf{x}\) в начальный момент времени \(t_0\): \(\textbf{x}(t_0) = \textbf{x}_0\). Пользуясь им нужно определить состояние системы \(\textbf{x}\) во все последующие моменты времени.

fig1.png

Рис. 1. Производная функция \(f(\textbf{x},t)\) определяет векторное поле

Задачу с начальными условиями легко представить наглядно. В двумерном случае \(\textbf{x}(t)\) представляет собой траекторию движения точки \(\textbf{p}\) на плоскости. В любой точке \(\textbf{x}\) можно вычислить функцию \(f\) и получить в результате двумерный вектор, так что \(f\) определяет векторное поле на плоскости (рис. 1). Этот вектор представляет собой скорость, с которой должна двигаться точка \(\textbf{p}\), если она будет проходить через \(\textbf{x}\) (она может проходить через эту точку, а может и не проходить). Представим, что \(f\) переносит \(\textbf{p}\) из точки в точку, подобно океанскому течению. Как только мы поместим \(\textbf{p}\) на плоскость, "течение" подхватит ее. Место, куда в конечном итоге попадет \(\textbf{p}\), зависит от того, куда мы поместили ее в начале, но как только это сделано, все последующее движение точки \(\textbf{p}\) определяется функцией \(f\). Траектория движения \(\textbf{p}\) под действием \(f\) называется интегральной кривой векторного поля (см. рис. 2).

fig2.png

Рис. 2. Задача с начальным условием. Стартуя из точки \(\textbf{x}_0\), частица \(\textbf{p}\) движется со скоростью, определяемой векторным полем

Мы записали функцию \(f\) как зависящую от \(\textbf{x}\) и от \(t\), но производная функция может как зависеть, так и не зависеть от времени. Если она зависит от времени, то в нашем наглядном представлении движется не только \(\textbf{p}\), движется все векторное поле, так что скорость \(\textbf{p}\) зависит не только от того, где она находится, но и от того, каким путем она попала в это место. Зависимость \(\dot{\textbf{x}}\) в этом случае проявляется двойственно: во-первых, вектор производной изменяется сам, во-вторых, точка \(\textbf{p}\), двигаясь по траектории \(\textbf{x}(t)\), "видит" разные векторы производных в разные моменты времени. Эта двойственная зависимость не сбивать с толку, если представлять себе частицу, движущуюся по "волнистому" векторному полю.

2. Численные решения

Обычно во вводных курсах по теории дифференциальных уравнений рассматривают уравнения, имеющие аналитические решения, т.е. предполагается, что искомую функцию можно записать в виде формулы. Например, для дифференциального уравнения \(\dot{x}=-kx\), где \(\dot{x}\) означает производную от \(x\) времени, решением является \(x=e^{-kt}\).

Мы, напротив сосредоточимся исключительно на численных решениях, которые заключаются в выполнении дискретных шагов по времени, начиная с начального момента \(\textbf{x}(t_0)\). Шаг по времени выполняется так: вычисляется функция \(f\), чтобы найти приблизительную величину приращения \(\textbf{x}\) за время \(\Delta t\)\(\Delta\textbf{x}\), затем \(\Delta\textbf{x}\) прибавляется к \(\textbf{x}\), чтобы получить новое состояние системы. Теперь нам известно состояние системы \(\textbf{x}(t_1)\) в момент \(t_1 = t+\Delta t\). Принимая его за начальное, выполним следующий шаг по времени и найдем состояние системы в момент \(t_2 = t+2\Delta t\) и т.д. Функцию \(f\) можно рассматривать как черный ящик, на вход которого подаются значения \(\textbf{x}\) и \(t\), а на выходе получается значение \(\dot{\textbf{x}}\). Различные методы численного решения используют одно или несколько вычислений производной функции на каждом шаге по времени.

2.1 Метод Эйлера

Простейший метод численного решения ОДУ называется методом Эйлера. Пусть начальное условие в момент времени \(t_0\) обозначается как \(\textbf{x}_0=\textbf{x}(t_0)\), а оценка \(\textbf{x}\) в момент времени \(t_0+h\) — как \(\textbf{x}(t_0+h)\), где \(h\)длина шага. Сделаем шаг в направлении, указанном вектором производной, и вычислим \(\textbf{x}(t_0+h)\)

$$ \textbf{x}(t_0+h) = \textbf{x}_0 + h \dot{\textbf{x}}(t_0) $$

— в этом и состоит метод Эйлера.

Чтобы понять как работает метод Эйлера, представьте картинку двумерного векторного поля (рис. 3). Вместо реальной интегральной кривой, точка \(\textbf{p}\) будет двигаться по ломаной, каждое звено которой соответствует шагу численного решения. Направление звена определяется оценкой вектора \(f\), сделанной в начале шага, длина звена равна длине этого вектора, умноженной на \(h\).

fig3.png

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

Метод Эйлера прост, однако недостаточно точен. Рассмотрим двумерную функцию \(f\), интегральные кривые которой представляют собой концентрические окружности (рис. 4). Тогда точка \(\textbf{p}\), движимая \(f\), должна все время оставаться на одной и той же окружности. Вместо этого, с каждым шагом расчета по методу Эйлера, \(\textbf{p}\) будет двигаться по прямой в сторону окружности большего радиуса, так что ее траектория превратится в разворачивающуюся спираль. Уменьшив длину шага, можно замедлить этот процесс, но не исключить его полностью.

fig4.png

Рис. 4. Вверху: реальные интегральные кривые являются концентрическими окружностями, а метод Эйлера дает вместо них разворачивающиеся спирали, потому что каждый шаг, сделанный по касательной к окружности, переводит точку \(\textbf{p}\) на окружность большего радиуса. Уменьшение шага не решает проблему, а только уменьшает скорость накопления ошибки.

Внизу: слишком большая длина шага заставляет метод Эйлера расходиться.

Кроме того, метод Эйлера может быть неустойчивым. Рассмотрим функцию \(f=-kx\), которая, как мы видели, дает интегральные кривые экспоненциально убывающие к нулю. Для достаточно малых длин шага мы получим приемлемые приближения, однако, если \(h > 1/k\), то мы получим \(|\Delta x| > |x|\), и приближенное решение будет колебаться в окрестности нуля. Если же длина шага \(h = 2/k\), то эти колебания будут возрастающими. Как говорят в таких случаях, система "взрывается".

Наконец, метод Эйлера не является даже эффективным. В большинстве методов численного решения ОДУ основное время тратится на вычисление производной, так что вычислительная "стоимость" одного шага определяется тем, сколько раз на нем вычислялась производная. В методе Эйлера производная на шаге вычисляется всего один раз, однако реальная эффективность метода определяется не только эффективностью одного шага, но и числом сделанных шагов — при этом надо сохранить устойчивость и обеспечить требуемую точность расчетов. Поэтому более продвинутые методы численного решения выигрывают у метода Эйлера за счет того, что хотя и требуют четырех или пяти вычислений производной за один шаг, но зато позволяют делать значительно более длинные шаги.

Чтобы понять, как можно улучшить метод Эйлера, рассмотрим откуда берется ошибка, производимая этим методом. Ключ к пониманию дает разложение в ряд Тейлора. Предположив, что \(\textbf{x}(t)\) — гладкая функция, выразим ее значение в конце шага через бесконечную сумму, состоящую из значения этой функции и ее производных в начале шага

$$ \textbf{x}(t_0+h) = \textbf{x}_0 + h \dot{\textbf{x}}(t_0) + \frac{h^2}{2!} \ddot{\textbf{x}}(t_0) + \frac{h^3}{3!} \frac{\partial^3\textbf{x}}{\partial t^3} + \ldots + \frac{h^n}{n!} \frac{\partial^n\textbf{x}}{\partial t^n} + \ldots $$

Как видно, формула метода Эйлера получается, если обрезать этот ряд, сохранив только первые два слагаемых в правой части равенства. Это означает, что метод Эйлера будет точен только в случае, если все производные порядков выше первого будут равны нулю, то есть когда функция \(\textbf{x}(t)\) линейна. Погрешность метода Эйлера равна "обрезанной" части разложения в ряд Тейлора и определяется, в первую очередь, слагаемым \(({h^2}/{2!}) \ddot{\textbf{x}}(t_0)\). Следовательно, мы можем записать погрешность шага метода Эйлера как \(O(h^2)\) (читается: "порядка \(h\) в квадрате").

Предположим, что мы уменьшили шаг вдвое и теперь он равен \(h/2\). Хотя при этом мы получаем ошибку, равную лишь четверти той, что была при длине шага \(h\), но теперь нам нужно сделать вдвое больше шагов, чтобы покрыть такой же по времени интервал. Это означает, что погрешность, накапливаемая на интервале от \(t_0\) к \(t_1\), зависит от длины шага \(h\) линейно — \(O(h)\) — меньше делается шаг и пропорционально уменьшается погрешность. На практике это приводит к тому, что для достижения хорошей точности нужно сделать слишком много шагов.

2.2 Метод средней точки2

Если бы мы, наряду с \(\dot{\textbf{x}}\), могли вычислить и \(\ddot{\textbf{x}}\), то можно было бы достичь точности \(O(h^3)\) вместо \(O(h^2)\), сохранив в разложении в ряд Тейлора еще одно слагаемое

\begin{equation} \textbf{x}(t_0+h) = \textbf{x}_0 + h \dot{\textbf{x}}(t_0) + \frac{h^2}{2} \ddot{\textbf{x}}(t_0) + O(h^3). \label{eq:midpoint1} \end{equation}

Вспомним, что производная по времени \(\dot{\textbf{x}}\) задается функцией \(f(\textbf{x}(t),t)\). Для простоты будем считать, что \(f\) не зависит явно от времени, так что \(\dot{\textbf{x}} = f(\textbf{x}(t))\). Правило дифференцирования сложной функции дает в этом случае

$$ \ddot{\textbf{x}} = \frac{\partial f}{\partial \textbf{x}}\dot{\textbf{x}} = f^\prime f. $$

Чтобы избежать вычисления \(f^\prime\), которое может быть сложным и требовать большого расхода времени, найдем приближенное выражение для \(\ddot{\textbf{x}}\) через \(f\) (не содержащее \(f^\prime\)) и подставим его в \eqref{eq:midpoint1}. При этом порядок точности \(O(h^3)\) сохраниться.

Для этого разложим \(f\) в ряд Тейлора

\begin{equation} f(\textbf{x}_0 + \Delta\textbf{x}) = f(\textbf{x}_0) + \Delta\textbf{x} f^\prime (\textbf{x}_0) + O(\Delta\textbf{x}^2). \label{eq:midpoint2} \end{equation}

Введем в него \(\ddot{\textbf{x}}\), подобрав величину шага равной

$$ \Delta\textbf{x} = \frac{h}{2}f(\textbf{x}_0), $$

так что \eqref{eq:midpoint2} будет иметь вид

$$ f(\textbf{x}_0 + \frac{h}{2}f(\textbf{x}_0)) = f(\textbf{x}_0) + \frac{h}{2}f(\textbf{x}_0) f^\prime (\textbf{x}_0) + O(h^2) = f(\textbf{x}_0) + \frac{h}{2} \ddot{\textbf{x}}(t_0) + O(h^2). $$

Отсюда, умножая обе части равенства на \(h\) (что превращает O(h^2) в O(h^3)) и перегруппировывая слагаемые, получим

$$ \frac{h^2}{2} \ddot{\textbf{x}}(t_0) + O(h^3) = h(f(\textbf{x}_0 + \frac{h}{2}f(\textbf{x}_0))) - f(\textbf{x}_0). $$

Подстановка правой части этого равенства в \eqref{eq:midpoint1} дает формулу для численного решения

$$ \textbf{x}(t_0+h) = \textbf{x}_0 + h(f(\textbf{x}_0 + \frac{h}{2}f(\textbf{x}_0))). $$

В соответствии с нею сначала выполняется шаг по методу Эйлера, а затем вычисляется оценка второй производной в средней точке шага для получения нового значения \(\textbf{x}\). Поэтому этот метод называют методом средней точки. Погрешность шага этого метода оценивается как \(O(h^3)\), но для этого требуется два вычисления производной. Иллюстрация работы метода приведена на рис. 5.

fig5.png

Рис. 5. Метод средней точки относится к методам 2-го порядка, т.е. его погрешность на интервале расчета — \(O(h^2)\). Шаг метода состоит из a) выполнения шага по методу Эйлера; b) оценки производной в средней точке шага; c) вычисления новой точки интегральной кривой.

Достигнутая погрешность \(O(h^3)\) — далеко не предел. Вычислив \(f\) еще несколько раз во время шага, мы сможем исключить производные более высоких порядков. Самый популярный метод, в котором используется эта идея, называется методом Рунге-Кутта 4-го порядка и дает погрешность на шаге \(O(h^5)\) (метод средней точки можно назвать методом Рунге-Кутта 2-го порядка). Мы не будем выводить, а просто приведем здесь формулы метода Рунге-Кутта 4-го порядка

$$ \begin{array} k_1 &= hf(\textbf{x}_0, t_0), \\ k_2 &= hf(\textbf{x}_0+\frac{k_1}{2}, t_0+\frac{h}{2}), \\ k_3 &= hf(\textbf{x}_0+\frac{k_2}{2}, t_0+\frac{h}{2}), \\ k_4 &= hf(\textbf{x}_0+k3, t_0+h), \\ \textbf{x}(t_0+h) &= \textbf{x}_0+\frac{k_1}{6}+\frac{k_2}{3}+\frac{k_3}{3}+\frac{k_4}{6}. \end{array} $$

3. Адаптивный выбор шага

Какой бы метод численного решения мы не использовали, остается вопрос выбора длины шага. В идеале, хотелось бы сделать \(h\) настолько большим, насколько это возможно, чтобы обеспечить заданную точность расчетов и избежать неустойчивости метода. С фиксированной длиной шага мы можем двигаться вперед не быстрее, чем нам позволит это сделать наиболее трудная для вычислений часть \(\textbf{x}(t)\). Хотелось бы, чтобы длина шага на интервале расчетов уменьшалась, когда это нужно для обеспечения точности, а затем, когда трудное место пройдено, вновь увеличивалась. Т.е. длина шага должна выбираться адаптивно.

Рассмотрим как осуществляется адаптивный выбор шага в методе Эйлера. Предположим, что у нас задана длина шага \(h\) и нужно определить, насколько ее следует увеличить или уменьшить не теряя точности.

Допустим, мы вычислили две оценки величины \(\textbf{x}(t_0+h)\). Первая (\(\textbf{x}_a\)) получена вычислением шага \(h\) от \(t_0\) к \(t_0 + h\), вторая (\(\textbf{x}_b\)) — вычислением двух последовательных шагов длины \(h/2\) от \(t_0\) к \(t_0 + h\). Обе оценки \(\textbf{x}_a\) и \(\textbf{x}_b\) отличаются от истинного значения \(\textbf{x}(t_0+h)\) не более чем на \(O(h^2)\). Следовательно и друг от друга \(\textbf{x}_a\) и \(\textbf{x}_b\) отличаются не более чем на \(O(h^2)\).

Определим меру текущей погрешности \(e\) как

$$ e = |\textbf{x}_a - \textbf{x}_b|. $$

Предположим, что требуется обеспечить погрешность на шаге расчета не более \(10^{-4}\), а текущая погрешность \(e\) составляет всего \(10^{-8}\). Т.к. погрешность пропорциональна \(h^2\), мы можем увеличить длину шага до

$$ \left( \frac{10^{-4}}{10^{-8}} \right)^\frac{1}{2} h = 100h. $$

Напротив, если текущая погрешность составляет \(10^{-3}\) мы должны уменьшить шаг до

$$ \left( \frac{10^{-4}}{10^{-3}} \right)^\frac{1}{2} h \approx 0.316h. $$

Адаптивный выбор шага — очень полезная техника.

4. Реализация

ОДУ, которые мы будем решать, могут представлять собой модели самых разных физических систем: материальных точек, соединенных пружинами, твердых тел или деформируемых объектов. Мы хотим реализовать решатели3 ОДУ так, чтобы они не зависели от моделей, для которых они применяются. Это позволит легко заменять один решатель на другой и упростит повторное использование их кода4. К счастью, достичь такого уровня модульности нетрудно, т.к. каждый решатель можно представить в виде короткого и однотипного набора операций.

Обычно, система моделируемых с помощью ОДУ объектов представляет собой структуру (класс). Наш подход заключается в том, чтобы написать для этой структуры код, выполняющий необходимые решателю операции, а затем написать решатель, алгоритм которого использует эти операции.

С точки зрения решателя, система с которой он работает представляет собой функцию-"черный ящик" \(f(\textbf{x},t)\). Решатель должен иметь возможность вычислить \(f\) при любых значениях \(\textbf{x}\) и \(t\) и вернуть обновленные значения \(\textbf{x}\) и \(t\) после окончания расчета. Чтобы обеспечить выполнение этих операций, объект, представляющий ОДУ, должен уметь обрабатывать следующие запросы решателя:

  • Возвратить значение \(\dim(\textbf{x})\). Поскольку \(\textbf{x}\) и \(\dot{\textbf{x}}\) — векторы, решатель должен знать их длину чтобы выделить нужное место в памяти, выполнить векторные операции и т. п.
  • Получить/задать значения \(\textbf{x}\) и \(t\). Решатель должен возвращать окончательные значения этих величин. К тому же, при использовании многошаговых методов, решатель должен задавать промежуточные значения \(\textbf{x}\) и \(t\) для вычисления производной.
  • Вычислить \(f\) при заданных \(\textbf{x}\) и \(t\).

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

Ссылки

  1. W.H. Press, B.P. Flannery, S. A. Teukolsky, and W. T. Vetterling. Numerical Recipes in C. Cambridge University Press, Cambridge, England, 1988.


  1. Другое название задач с начальными условиями — задачи Коши. 

  2. Midpoint method. 

  3. Функции, реализующие выбор шага и тот или иной численный метод решения ОДУ. 

  4. Например, решатель с адаптивным шагом может вызывать решатель с фиксированным шагом. 



Комментарии

comments powered by Disqus