Известно большое количество методов численного решения дифференциальных уравнений в общем виде [1,2]. Многие из них реальзованы программно и доступны в виде готовых модулей. К сожалению, обобщённые методы не всегда дают достаточно точный результат, т.к. существует естественное ограничение вычислительных мощностей. Также, применение готовых модулей часто выводит процесс решение за пределы разрядности процессора, что сильно ограничивает диапазон вводимых значений.

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

Для рассмотрения предлагаемого метода нужно вспомнить, что процесс численного решения состоит из итераций — многократных повторений некоторого вычисления [3], а номер шага такой итерации обозначается символом \(i\). Сразу же рассмотрим пример в виде следующего дифф. уравнения: \[\dot Y = f(x)\,Y + a \qquad (1)\] Здесь: \(Y=Y(x)\), \(\dot Y\) — некоторое приращение \(Y(x)\) к изменению координаты \(x\), т.е. \[\dot Y = {dY \over dx} \approx {\Delta Y \over \Delta x}, \] \(f(x)\) — некоторая функция от этой координаты, а \(a\) — константа. К слову, если \(f(x)\) будет довольно сложной зависимостью, то аналитическое решение такого уравнения будет невозможно, а вот решение численными методами — всегда и с любой точностью.

Для этого представим \(\Delta Y\), как разность между текущим значением и значением на предыдущем шаге итерации: \[\Delta Y = Y(x_i) - Y(x_{i-1}) \qquad (2)\] При этом координата изменится на \(\Delta x = x_i - x_{i-1}\), и такое изменение координаты будет постоянным и на любом шаге. От \(\Delta x\) зависит точность всего метода и общее число итераций: \(i \in 1, 2, 3, ..., N\), а само значение координаты находится просто: \(x_i = i\,\Delta x\).

Тогда наше уравнение запишется так: \[{Y(x_i) - Y(x_{i-1}) \over \Delta x} = f(x_{i})\,Y(x_{i}) + a \qquad (3)\] Остаётся только найти текущие значения функции, которые и будут являться выходными данными численного решения для каждого шага итерации: \[Y(x_i) = {Y(x_{i-1}) + a\,\Delta x \over 1 - f(x_i)\,\Delta x} \qquad (4)\] Также, необходимо определиться со значением \(Y(x_{i-1})\) на самом первом шаге, т.е. когда \(i=1\), а \(x_0=0\). Это будет начальное значение нашей функции, которое в мат. анализе обознается, как \(Y(0)\). Тогда общая запись численного решение будет такой: \[Y(x_i) = {Y(x_{i-1}) + a\,\Delta x \over 1 - f(x_i)\,\Delta x}, \quad Y(x_0) = Y(0) \qquad (5)\] При всей простоте найденной функции — это самый быстрый и устойчивый алгоритм для численного решения этого уравнения!

Любителям свободной энергии будет интересно узнать, что в линейном уравнении Бернулли скрыт один из способов её получения. Найдём численное решение этого уравнения: \[g(x)\,\dot Y = f(x)\,Y + h(x) \qquad (6)\] Вначале предлагаем посмотреть на аналитический подход [4], а затем, по предложенному выше методу, найти численное решение. Проделав почти те же действия тожно получить на выходе следующий алгоритм: \[Y(x_i) = {g(x_i)\,Y(x_{i-1}) + h(x_i)\,\Delta x \over g(x_i) - f(x_i)\,\Delta x}, \quad Y(x_0) = Y(0) \qquad (7)\] Для сравнения, если искать численное решение этого уравнения с помощью подхода [4], то оно получится таким: \[Y(x_i) = \left[ Y(x_{i-1}) + \frac{h(x_i)}{g(x_i)} \Delta x \right] \left[ 1 + \frac{f(x_i)}{g(x_i)} \Delta x \right], \quad Y(x_0) = Y(0) \qquad (8)\] Интересно то, что на практике, программирование с помощью алгоритма (7) даёт более точные значения на выходе, чем (8), а сам процесс вычисления оказывается более устойчивым.

Это решение должно раскрыть всю мощь этого метода. Для этого возьмём следующее уравнение: \[\dot Y = f(x) \exp(a\,Y) + h(x) \qquad (9)\] и попробуем найти его численное решение по предложенной методике. Но сначала вспомним, что экспоненциальная функция, при малых изменениях, может быть разложена в ряд Тейлора [5]: \[\exp(a\,Y_i) = \exp(a\,Y_{i-1})\exp(a\,\Delta Y) = \exp(a\,Y_{i-1})(1 + a(Y_i-Y_{i-1})) \qquad (10)\] Здесь, для лучшего восприятия введены такие упрощения: \(Y_i = Y(x_i), Y_{i-1} = Y(x_{i-1})\). Теперь, по предложенной методике, составим общее уравнение: \[{Y_i - Y_{i-1} \over \Delta x} = f(x_i) \exp(a\,Y_{i-1})[1 + a(Y_i-Y_{i-1})] + h(x_i) \qquad (11)\] и выведем из него искомую функцию: \[Y_i = Y_{i-1} + {f(x_i) \exp(a\,Y_{i-1}) + h(x_i) \over 1 + a\,f(x_i) \exp(a\,Y_{i-1}) \Delta x } \Delta x \qquad (12)\]

На следующей странице мы расскажем о численном решении некоторых дифф. уравнений второго порядка.