5.2.4. Обыкновенные дифференциальные уравнения

Как известно, всякую задачу Коши для системы обыкновенных дифференциальных уравнений в нормальном виде можно представить системой уравнений первого порядка, поэтому ограничимся рассмот-рением задачи . Точное (аналитическое) решение удается найти лишь в исключительных случаях, поэтому широко распростра-нены приближенные (численные) методы. Простейший - способ Эйле-ра, использующий правую часть уравнения для проведения касатель-ной к решению: если t _i = ih (i = 0, 1, ...), то x _i = x(t _i) ≈ x _{i - 1} +

+ hf(t _{i - 1}, x _{i - 1}), что легко реализуется кодом (n задано, h выбрано):

i 0 . . n - 1 t _i i  h x _i x0 if i = 0

x _{i - 1} + h  f(t _{i - 1}, x _{i - 1}) otherwise .

Этот алгоритм плохо сходится даже для "хороших" задач, однако он легко улучшается, если перемещаться по направлению среднему из на-правлений в соседних узлах (модифицированный метод), полагая k _i = x _i + hf(t _i , x _i), x _{i + 1} = x _i + 0.5h(f(t _i , x _i) + f(t _{i + 1} , k _i)), i = 0, 1, ... . Рунге и Кутт довели модификацию до логического конца. Приведем

пример метода Рунге - Кутта четвертого порядка:

x _i p ← x0 if i = 0

otherwise

k1 ← f(t _{i - 1}, x _{i - 1})h

k2 ← f(t _i - ₁+ 0.5h, x _{i - 1} + 0.5k1)h

k3 ← f(t _{i - 1} + 0.5h, x _{i - 1} + 0.5k2)h

k4 ← f(t _{i - 1} + h, x _{i - 1} + k3)h

p ← x _{i - 1} + (k1 + 2 k2 + 2 k3 + k4)6 ^{- 1}

p .

MathCad имеет несколько встроенных функций для решения задачи Коши различными методами (Bulstoer, Rkadapt, rkfixed, Stiffb и др.), а также блок Given ... Odesolve для определения функции-решения.

Пример 1. Найдем решение уравнения y′ = на [0, 5] с начальным условием у(0) = 1 методом rkfixed с шагом h = 0.25 ( в 20 точках) и блоком Given ... Odesolve и сравним результат с точным ре-шением y(x) = (Ctrl + F7 вводит знак производной):

F(t, x) D(t, y) F(t, y ₀) r(t)

y₀ = 1 M rkfixed(y, 0, 5, 20, D) s M

Given

y′(x) = y(0) = 1 (клавиша "Ctrl+F7")

f odesolve(x,5) | r(s) - M | = 4.311  10 ^{- 6} | r(s) - | = 2.769  10 ^{- 5}

Здесь M = {t _i} - вектор значений аргументов (сетка), M = {y _i} - вектор значений решающей функции, значение функции на векторе есть вектор соответствующих значений функции, однако результат ра-боты Odesolve требует в этом случае векторизации (знак "→ ").

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

Пример 2. Зависимость величины горизонтального отклонения маятника единичной массы на тонкой нити длины r от времени мож-но описать нелинейным дифференциальным уравнением

x''(t) + + k x'(t) = 0,

где g – ускорение свободного падения, k – коэффициент сопротивле-ния среды. Пусть для определенности r = 10, k = 0.2 и x(0) = 7, x'(0) = 0. Для приведения модели к нормальной системе уравнений первого порядка введем новые переменные: обозначим x = y₀, x' = y₁ . В результате получаем систему уравнений первого порядка

или в векторном виде y' = D(t, y) с правой частью

D(t, y) = и с начальным условием y = . Те--перь достаточно обратиться к какой-нибудь встроенной функции, ре-шающей эту систему, например, положив z := rkfixed(y, 0, 20, 100, D). Матрица z будет иметь 101 значение приближенного решения задачи на отрезке [0, 20], первый столбец матрицы – точки отсчета времени, второй – соответствующие значения у₀ (или х), третий – значения у₁ (или x'). Задача решена. Можно строить график z^{‹ 1 ›} относительно

z^{‹ 0 ›} (зависимость x(t)) или z^{‹ 1 ›} относительно z^{‹ 2 ›} (фазовый портрет,

см. рис. 5).

амплитуда колебаний фазовый портрет

Рис. 5

Решение краевых задач существенно сложнее и часто сводится к

решению подходящих задач Коши (см. встроенные функции MathCad: bvalfit, sbval и др.). Блок Given ... Odesolve справляется и с некоторы-ми краевыми задачами. Многие задачи (особенно линейные) можно приближенно решить уже известными средствами: на равномерной сетке {x _i} c шагом h производные удобно заменять разностными со-отношениями (см. предыдущий пункт) и тогда задача сводится к сис-теме (вообще говоря нелинейных) алгебраических уравнений, которые можно решать методом Ньютона или (для линейных) прогонки (см. п. 2.1.1).

Пример 3. Рассмотрим линейную краевую задачу

p(x)y′′ + q(x)y′ + r(x)y = f(x) при условиях y(a) = u, y(b) = v.

Фиксируем n (число узлов сетки), положим h = , x _i = a + ih,

y _i = y(x _i), y′(x _i ) ≈ 0.5h ^{- 1}(y _{i + 1} - y _{i - 1}), y′′(x _i) ≈ h ^{- 2}(y _{i - 1}- 2y _i + y _{i + 1}) и уравнение (в этих терминах) становится чисто алгебраическим для каждого i, система таких уравнений линейная, 3 - диагональная, для ее решения можно использовать метод прогонки с матрицей А (три

столбца) и вектором правых частей b , которые имеют вид:

k 0 . . 2 A_{j, k} 0 if (j + k = 0) + (j + k = n)

otherwise

2p(x _{j + 1}) - hq(x _{j + 1}) if k = 0

2h ²r(x _{j + 1}) - 4p(x _{j + 1}) if k = 1

2p(x _{j + 1}) + hq(x _{j + 1}) otherwise

b _j r ← 2h ² f(x _{j + 1}), s ← 2p(x _{j + 1}) - hq(x _{j + 1})

r ← r - s  u if j = 0

r ← r - s  v if j = n - 2

r

Для конкретной задачи x (3x ² +1)y′′ + 2y′ - 6xy = 4 - 12x ² с условиями y(1) = y(2) = 0 положим F(x) 2  x - 0.75  (x ² + 1) - (точное решение)

n 7 a 1 b 2 u 0 v 0

p(x) x(3x ² + 1) q(x) 2 r(x) - 6x f(x) 4 - 12x ²

h := i 0 . . n x _i a + i  h z Prog(A, b) y _i if(i = 0, u, z _{i - 1})

Given

p(x)y′′(x) + q(x)y′(x) + r(x)y(x) = f(x)

y(1) = 0 y(2) = 0

g odesolve(x,2) |y - F(x)| = 7.106  10 ^{- 4} |g(x - F(x)| = 3.895  10 ^{- 4}