Известно, что с помощью дифференциальных уравнений можно описать задачи движения системы взаимодействующих материальных точек, химической кинетики, электрических цепей, и др. Конкретная прикладная задача может приводить к дифференциальному уравнению любого порядка, или к системе уравнений любого порядка. Так как любое ОДУ порядка p

u (p) (x) = f(x, u, u / , u // , … u (p+1) ) заменой u (k) (x) = y_k(x) можно свести к эквивалентной системе из p уравнений первого порядка, представленных в каноническом виде:

Покажем такое преобразование на примере уравнения Бесселя: .

В окрестности точки x₀ функцию u(x) разложим в ряд Тейлора:

Идея этого и последующих методов основывается на том, что если f (x , u) имеет q непрерывных производных, то в разложении можно удержать члены вплоть до O(h q+1 ), при этом стоящие в правой части производные можно найти, дифференцируя (3) требуемое число раз. В случае метода Эйлера ограничимся только двумя членами разложения.

Пусть h – малое приращение аргумента. Тогда (4) превратится в

Так как в соответствии с (3) u / (x₀) = f(x₀, u₀), то .

Теперь приближенное решение в точке x₁ = x₀+h можно вновь рассматривать как начальное условие, т.е. организуется расчет по следующей рекуррентной формуле:

где y₀ = u₀, а все y_k – приближенные значения искомой функции (см. рисунок). В методе Эйлера происходит движение не по интегральной кривой, а по касательной к ней. На каждом шаге касательная находится уже для новой интегральной кривой (что и дало название методу – метод ломаных), таким образом ошибка будет возрастать с отдалением x от x₀.

При приближенное решение сходится к точному равномерно с первым порядком точности. То есть, метод дает весьма низкую точность вычислений: погрешность на элементарном шаге h составляет ½ h 2 y // (½(x_k+x_k+1)), а для всей интегральной кривой порядка h 1 . При h = const для оценки апостериорной погрешности может быть применена первая формула Рунге, хотя для работы метода обеспечивать равномерность шага в принципе не требуется.

Метод Эйлера легко обобщается для систем ОДУ. При этом общая схема процесса (5) может быть записана так:

где i = 1… m – число уравнений, k – номер предыдущей вычисленной точки.

Данный метод базируется на предыдущем, однако здесь апостериорная погрешность контролируется на каждом шаге вычисления. Как и в предыдущем случае, рассматриваем задачу Коши на сетке с постоянным шагом h. Грубое значение вычисляется по формуле (5) и затем итерационным циклом уточняется по формуле

где m – номер итерации. Итерационный цикл повторяется до тех пор, пока . Данная формула также легко обобщается на решение систем ОДУ. Априорная погрешность метода при m = 1 на каждом шаге порядка h 3 .

Увеличение точности решения ОДУ из предыдущей задачи при заданном шаге h может быть достигнуто учетом большего количества членов разложения функции в ряд Тейлора. Для метода Рунге-Кутты второго порядка следует взять три первых коэффициента, т.е. обеспечить:

Переходя к приближенному решению y ≈ u и заменяя производные в (8) конечными разностями, получаем в итоге следующее выражение:

где 0 ≤ α ≤ 1– свободный параметр. Можно показать, что если f(x,u) непрерывна и ограничена вместе со своими вторыми производными, то решение, полученное по данной схеме, равномерно сходится к точному решению с погрешностью порядка h 2 .

Для параметра α наиболее часто используют следующие значения:

. Графически это уточнение можно интерпретировать так: сначала по схеме ломаных делается шаг h /2, и находится значение . В найденной точке определяется наклон касательной к интегральной кривой, который и будет определять приращение функции для целого шага, т.е. отрезок [ AB ] (см. рисунок) будет параллелен касательной, проведенной в точке (x_k + h/2, y(x_k+h/2)) к интегральной кривой.

. Можно представить, что в этом случае по методу Эйлера сначала вычисляется значение функции и наклон касательной к интегральной кривой в точке x_{k +1}. Затем находится среднее положение касательной из сравнения соответствующих наклонов в точках x _k и x_{k +1}, которое и будет использоваться для расчета точки y_k+1.

Данная схема является наиболее употребительной. Здесь в разложении функции в ряд Тейлора учитываются члены до h 4 включительно, т.е. погрешность на каждом шаге пропорциональна h 5 . Для практических вычислений используются следующие соотношения, обобщенные в данном случае на решение системы ОДУ:

, где i = 1… p, p – число уравнений в системе.

; k – номер точки, для которой осуществляется расчет;

К достоинствам метода следует отнести высокую точность вычислений. Схемы более высокого порядка точности практически не употребляются в силу своей громоздкости. Также немаловажно, что метод является явным, т.е. значение y_{k +1} вычисляется по ранее найденным значениям за известное заранее число действий.

Все представленные выше схемы допускают расчет с переменным шагом. Например, шаг можно уменьшить там, где функция быстро изменяется, и увеличить в обратном случае. Так, метод Рунге-Кутты-Мерсона позволяет оценивать погрешность на каждом шаге и, в зависимости от полученной оценки принимать решение об изменении шага. Автоматический выбор шага позволяет значительно сократить время вычислений.

Этот метод отличается от метода Рунге – Кутта четвертого порядка возможностью оценивать погрешность на каждом шаге и в зависимости от этого принимать решение об изменении шага. Один из вариантов формул:

;

Пусть задана максимальна погрешность . Если , h = h /2 , и n +1 цикл расчета повторяется (с точки x_n , y_n ) c новым шагом. Если же

Автоматический выбор шага позволяет значительно сократить время решения ОДУ.

Схема РКМ обобщается на системы ОДУ аналогично классической схеме Рунге – Кутта.

Метод основан на аппроксимации интерполяционными полиномами правых частей ОДУ.

Пусть с помощью любого из методов, рассмотренных выше, вычислено решение заданного дифференциального уравнения в точках x₁, x₂, x₃ (а в точке x₀ решение и так известно – поставлена задача Коши). Полученные значения функции обозначим как y₀, y₁, y₂, y₃, а значения правой части дифференциального уравнения как f ₀, f₁, f₂, f ₃, где f _k = f (x_k, y_k). Начиная с четвертой точки, на каждом шаге интегрирования дифференциального уравнения вычисления осуществляются по схеме

где P – прогноз решения; Е – вычисление f (x , y ); С – коррекция решения; m – количество итераций коррекции. Схемы такого типа называют «прогноз-коррекция»: это подразумевает сначала приблизительное вычисление решение по формуле низкого порядка, а затем уточнение с учетом полученной информации о поведении интегральной кривой.

Прогноз осуществляется по экстраполяционной формуле Адамса:

Коррекция осуществляется по интерполяционной формуле Адамса:

Количество итераций m ≤ p, где p – порядок используемого метода. В ходе каждой итерации решается нелинейное уравнение (11) относительно неизвестной y₄ (обычно методом простых итераций).

Иногда в методе Адамса используется схеме PECE на каждом шаге процесса интегрирования, т.е. осуществляется только одна коррекция. В силу сложности вычислений метод используется только в мощных программных пакетах численного анализа. Формулы метода также легко переносятся на решение систем ОДУ первого порядка.

Краевая задача – задача отыскания частного решения системы

(12) , 1 ≤ k ≤ p, на отрезке a ≤ x ≤ b, на котором дополнительные условия налагаются на значения функции u_k(x) более чем в одной точке этого отрезка. В качестве примера моно привести задачу нахождения статического прогиба u(x) нагруженной струны с закрепленными концами:

, a ≤ x ≤ b, u(a) = u(b) = 0. Здесь f (x) имеет смысл внешней изгибающей нагрузки на единицу длины струны, деленной на упругость струны.

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

Рассмотрим систему двух дифференциальных уравнений первого порядка

Сначала выберем некоторое значение , так что . Это уравнение мы можем решить и определить . Таким образом у нас появились два числа и β, которые будут определять задачу Коши , для системы (13), удовлетворять левому краевому условию, но не удовлетворять второму уравнению (14). Тем не менее, решим систему ОДУ с задачей Коши каким-либо из известных нам численных методов. Получим решение вида

(15) , зависящее от как от параметра.

Теперь мы должны каким-либо способом менять параметр до тех пор, пока не подберется такое значение, для которого будет выполнено условие (16) , т.е. правое краевое условие. Для этого случайным образом берут значения до тех пор, пока среди величин не окажется разных по знаку.

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

К уравнениям в частных производных приводят задачи газодинамики, теплопроводности, переноса излучения, электромагнитных полей, процессов переноса в газах, и др. Независимыми переменными в физических задачах обычно являются время t, координаты , скорости частиц . Пример – уравнение теплопроводности

где U – температура, – теплоемкость, – коэффициент теплопроводности, q – плотность источников тепла.

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

Рассмотрим в качестве примера одномерную задачу, близкую по смыслу к (17):

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

Введем равномерную прямоугольную сетку по x и t с шагом h и τ соответственно и заменить производные соответствующими разностными отношениями. Тогда

Здесь 1 ≤ k ≤ N-1 – число точек по координате x; 0 ≤ m ≤ M – число точек по координате t. Число неизвестных в (19) больше числа уравнений, недостающие уравнения выводятся из начальных и граничных условий:

; ; 0 ≤ m ≤ M.

Схема (19) содержит в каждом уравнении несколько неизвестных значений функции. Такие схемы называют неявными. Для фактического вычисления решения следует переписать схему так:

; . (20)

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

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

Операционное исчисление