Постановка задачи. Типы задач для ОДУ.

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

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

y ^/_k(x) = y _k+1(x) для 0 ≤ k ≤ p –2 (1)

y ^/_p-1(x) = f (x, y ₀, y ₁, … y _p–1), при этом y ₀(x) ≡ u (x). (2)

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

Предполагая тождественную замену y ₁(x) ≡ y (x) представим систему ОДУ в следующем виде:

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

Известно, что система p -го порядка имеет множество решений, которые в общем случае зависят от p параметров { C ₁, C ₂, … C _p}. Для определения значений этих параметров, т.е. для выделения единственного решения необходимо наложить p дополнительных условий на u _k(x). По способу задания условий различают три вида задач, для которых доказано существование и единственность решения. Это

1) Задача Коши. Задается координата u _k(x ₀) = u _k0 начальной точки интегральной кривой в (p +1) – мерном пространстве (k = 1… p). Решение при этом требуется найти на некотором отрезке x ₀ ≤ x ≤ x _max.

2) Краевая задача. Это задача отыскания частного решения системы ОДУ на отрезке a ≤ x ≤ b, в которой дополнительные условия налагаются на значения функции u _k(x) более чем в одной точке этого отрезка.

3) Задача на собственные значения. Кроме искомых функций и их производных в уравнение входят дополнительно m неизвестных параметров λ₁, λ₂, … λ_m, которые называются собственными значениями. Для единственности решения на некотором интервале необходимо задать p+m граничных условий.

В большинстве случаев необходимость численного решения систем ОДУ возникает в случае, когда аналитическое решение найти либо невозможно, либо нерационально, а приближенное решение (в виде набора интерполирующих функций) не дает требуемой точности. Численное решение системы ОДУ, в отличие от двух вышеприведенных случаев, никогда не покажет общего решения системы, так как даст только таблицу значений неизвестных функций, удовлетворяющих начальным условиям. По этим значениям можно построить графики данных функций или рассчитать для некоторого x > x ₀ соотвествующие u _k(x), что обычно и требуется в физических или инженерных задачах. При этом требуется, чтобы соблюдались условия корректно поставленной задачи.

 

Метод Эйлера

Рассмотрим задачу Коши для уравнения первого порядка:

u ^/(x) = f (x, u(x)), x ₀ ≤ x ≤ x _max, u (x ₀) = u ₀. (3)

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

(4)

x0

f(x)-точное решение

x

y0

x1=x0+h

x2

y

Поле интегральных кривых

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

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

.

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

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

(5),

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

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

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

(6),

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

Алгоритм

13. Построить сетку с шагом . Значение шага выбирается из соображений требуемой точности, учитывая, что порядок точности метода Эйлера . Левая граница отрезка, на котором строится сетка, задана начальными условиями задачи Коши.

14. Решение дифференциального уравнения ищется в виде сеточной функции. Значение известно из начальных условий, все следующие значения рассчитываются по формуле Эйлера:

,

где – правая часть дифференциального уравнения.

15. Для проверки сравнить значения приближенного решения со значениями точного решения в узлах сетки.

 

Модифицированный метод Эйлера

16. Построить сетку с шагом . Значение шага выбирается из соображений требуемой точности, учитывая, что порядок точности модифицированного метода Эйлера . Левая граница отрезка, на котором строится сетка, задана начальными условиями задачи Коши.

17. Решение дифференциального уравнения ищется в виде сеточной функции. Значение известно из начальных условий, все следующие значения рассчитываются по модифицированной формуле Эйлера (формула Эйлера с пересчетом):

,

где – правая часть дифференциального уравнения;

– вычисляется предварительно по формуле простого метода Эйлера.

18. Для проверки сравнить значения приближенного решения со значениями точного решения в узлах сетки.

 

Метод Рунге–Кутта четвертого порядка точности