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

Простейшие алгоритмы численного интегрировании. Идеи численных методов проиллюстрируем на примере интегрирования уравнения первого порядка где х, t— зависимая и независимая переменные; fix, /) — известная (линейная или нелинейная) функция.

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

Уравнению (1) на интервале /₀-./1 соответствует интеграл (отсюда понятие интегрирование дифференциального уравнения)

Допустим, что:

=> начальное значение x (to) = хо известно;

=> функцияfix, t) на интервале /о…t изменяется незначительно.

В этом случае интеграл (2) с учетом (1) может быть представлен в следующих двух формах, известных соответственно как прямая (явная) и обратная (неявная) формулы Эйлера [8]:

где *i ~ x (t); h = t- fo; *о,, v'i — значение производных искомой функции x (t) при /₍) и (1).

С помощью рис. 1 можно дать геометрическую трактовку соотношениям.

(3) и (4). При получении прямой формулы Эйлера (3) искомая функция .*(/) аппроксимируется на интервале (шаге) интегрирования прямой линией, совпадающей с касательной к этой функции в точке jc₀= х (/о) (рис. , а), а при получении образной формулы Эйлера (4) — касательной в точке *(fi), выходящей из точки Хо = x (to) (рис. 1, б). В обоих случаях значение искомой функции x (t) определяется с некоторой ошибкой ? = x (t) — Х].

Представив (3) в виде соотношения (xj — Xo)/h = *'₀, а (4) — в виде соотношения (Х| - xo)/h = х и используя линейную комбинацию производных box’o + Ьх = (xi -xo)/h, можно получить при Ь₀ = Ь = ½ известную формулу трапеций

Рис. 1. Геометрическая трактовка прямой (а) и обратной (б) формулы Эйлера При численном интегрировании промежуток времени /₀…/_т, в пределах которого требуется получить решение л*(/) уравнения (1), разбивают на малые интервалы времени At = h, называемые шагом интегрирования. Полагают, что начальное значение x (to) = л*о известно и решение, притом единственное, существует на всем интересующем интервале времени. Соотношения (3)-(5) записывают в форме алгоритма.

позволяющего последовательно (задаваясь к = 0, 1,2,…) вычислить значения*), *2,д:з,… от начальной расчетной точки /] до конечной t_m с шагом h. Результаты численного интегрирования представляют собой таблицу соответствующих значений t_n wx_n (п = 0, 1,2, 3,…).

Алгоритм (6), позволяющий непосредственно определить решение в следующей точке /*+_ь называется явным алгоритмом. В неявных алгоритмах Эйлера (7) и трапеций (8) используется производная в следующей точке (/*+))> значение которой неизвестно. Возникает задача определения се приближенного значения, или предсказания. Одна из возможностей ее решения состоит в применении на каждом шаге прямой формулы Эйлера (6). На основе найденного предсказания можно рассчитать значение производной х' из (1) и использовать его при коррекции, которую можно выполнить по неявным формулам (7), (8). Повторным расчетом по формуле коррекции можно уточнить рассогласование начальных значений * и Хотя непосредственная необходимость в таких повторных расчетах отсутствует, в результате нескольких итераций (как показано в приведенных ниже примерах) можно.

получить более точное значение х". Итерации желательны, поскольку выполняются сравнительно просто. Кроме того, коррекция помогает в фиксации ошибок и управлении размером шага И.

Примеры. Для иллюстрации алгоритмов численного интегрирования решим неоднородное линейное дифференциальное уравнение х' = х + г, приняв шаг интегрирования Л = 0,025; x (t₀) = jc₀ = 1 (*о=0). Для оценки погрешности численного интегрирования (решения) воспользуемся точным решением дифференциального уравнения, имеющим следующий вид: л: = Зе¹ -1² — 2t — 2. Погрешность е определим как разность точного и полученного (численного) результатов.

Явный алгоритм Эйлера: jc*+i = ** + h (x_k + /*").

Шаг 1 (к = 0): = *₀ + Л (х₀ + Af) = 1 + 0,025(1 + О²) = 1,025.

Шаг 2 (к = 1): л:₂ = *, + Л (*, + г,²) = 1,025 + 0,025(1,025 + 0,025²) = 1,506 406 ит. д.

Неявный алгоритм Эйлера. На каждом шаге алгоритма будем использовать прямую формулу Эйлера, v*+iⁿ = х* + А (х* + f*²) для предсказания и обратную формулу Эйлера х&.i = x_k + h (x_k+ + 1²) для коррекции, используя при этом три итерации. Третья итерация является результатом и используется на следующем шаге алгоритма.

Шаг 1 (к = 0).

Предсказание: *|^П = x₀+h (xo + fo²) = 1 + 0,025(1+0²) = 1,025.

Итерация 1: *,⁽¹⁾ = *₀ + А (х,^п + /,²) = 1 + 0,025(1,025 + 0,025²) = = 1,256 406.

Итерация 2: лг,⁽²⁾ = хо + h (*,^(,) + /,²) = 1 + 0,025(1,256 406 + 0, 025²) = = 1,256 566.

Итерация 3: лг,⁽³⁾ = лг₀ + h (*,⁽²⁾ + /,²) = 1 + 0,025(1,256 566 + 0, 025²) = = 1,256 570.

Шаг 2(k= 1).

Предсказание: х₂^а = x_x+h (x_x + t_x²) = 1,256 570 + 0,025(1,256 570 + + 0,025²) = 1,51 314.

Итерация 1: *₂^(,) = х_х + h (х₂ + t₂) = 1,256 570 + 0,025(1,51 314 + + 0,05²)= 1,520 023.

Итерация 2: jc₂^<2) = Х + h (jc₂⁽¹⁾ + /₂²) = 1,256 570 + 0,025(1,520 023 + + 0,05²)= 1,520 195.

Итерация 3: л:₂⁽³⁾ = х_х + И (х₂^{2) + /₂²) = 1,256 570 + 0,025(1,520 195 + + 0,05²)= 1,520 199 ит. д.

Неявный алгоритм трапеций. На каждом шаге алгоритма предсказания будем рассчитывать по прямой формуле Эйлера = х_к + И (х_к + /д²), а для коррекции — использовать три итерации по.

формуле трапеций = ** + 0,5Л (дг* + t² + jc*+i + /*+1²) — Третья итерация является результатом и используется на следующем шаге алгоритма.

Шаг 1 (к = 0).

Предсказание: «Yiⁿ = x₀+h (xo + to²) = 1 + 0,025(1 + О²) = 1,025.

Итерация 1: *|^(,) = х₀ + 0,5/?(.y₀ + t₀² + х" + /,²) =.

= 1 + 0,50,025(1 + 0²+ 1,025 + 0,025²) = 1,235 203.

Итерация 2: *|^<2) = лг₀ + 0,5/?(jc₀ + /о² + *|^(,) + t²) =

= 1 +0,5−0,025(1 +0²+ 1,253 203 + 0,025²)= 1,253 243.

Итерация 3:, Vi⁽³⁾ = хо + 0,5Л (дго + to + Xi^(l) + t²) =

= 1 + 0,5−0,025(1 + 0²+ 1,253 243+ 0,025²) = 1,253 244 и т. д.

Результаты расчетов сведены в табл. 1.

Приведенные примеры показывают, что:

=> абсолютная (и относительная) ошибка возрастает с каждым шагом численного интегрирования. При этом ошибки явного и неявного алгоритма Эйлера имеют разные знаки (см. рис. 1). Неявный алгоритм трапеций, представляющий собой комбинацию прямой и обратной формул Эйлера, обладает меньшей погрешностью (этому ниже будет дано объяснение);

Таблица 1.

=> объелинение явных и неявных формул приводит к метолу предсказания и коррекции. При этом коррекцию можно рассматривать как способ итеративного решения в общем случае нелинейных алгебраических уравнений (относительно лг*+|) на каждом шаге неявного алгоритма. Неявные алгоритмы более трудоемки с точки зрения вычислений, однако обладают большей устойчивостью, о которой пойдет речь ниже.