Численное решение обыкновенных дифференциальных уравнений

Краткое введение. Дифференциальное уравнение первого порядка,разрешенное относительно производной, имеет вид

y ^' = f (x,y). (1)

Решением дифференциального уравнения (1) называется функция φ (x), подстановка которой в уравнение обращает его в тождество: φ ^' (x) ₌ f (x, φ (x)). График решения y = φ (x) называется интегральной кривой.

Задача Коши для дифференциального уравнения (1) состоит в том, чтобы найти решение уравнения (1), удовлетворяющему начальному условию

(2)

пару чисел (x₀, y₀ ) называют начальными данными. Решение задачи Коши называется частным решением уравнения (1) при условии (2).

Численное решение задачи Коши (1) - (2) состоит в том, чтобы получить искомое решение φ (x) в виде таблицы его приближенных решений для заданных значений аргумента x на некотором отрезке [ a, b ]:

X₀ = a, x ₁, x ₂,..., x _m = b (3)

Точки (3) называются узловыми точками, а множество этих точек называется сеткой на отрезке [ a, b ]. Будем использовать равномерную сетку с шагом h:

h = (b - a) / m; x _i - x _{i - 1} = h или x _i = x ₀ + ih (i = 1,..., m).

Приближенные значения численного решения задачи Коши в узловых точках x _i обозначим через y _i; таким образом,

(i = 1,..., m).

Для любого численного метода решения задачи (1) - (2) начальное условие (2) выполняется точно, т. е. .

Величина погрешности численного метода решения задачи Коши на сетке отрезка

[ a, b ] оценивается величиной

,

Говорят, что численный метод имеет p - й порядок точности по шагу h на сетке, если расстояние d можно представить в виде степенной функции от h:

, p > 0,

где c - некоторая положительная постоянная, зависящая от правой части уравнения (1) и от рассматриваемого метода. В данном случае очевидно, что когда шаг h стремится к нулю, погрешность d также стремится к нулю.

Метод Эйлера. Простейшим численным методом решения задачи Коши (1) - (2) является метод Эйлера.

Угловой коэффициент касательной к интегральной кривой в точке P₀ (x₀, y₀) есть

y ^'₀ = f (x₀,y₀ ).

Найдем ординату y₁ касательной, соответствующей абсциссе x ₁ = x ₀ + h. Так как уравнение касательной к кривой в точке P ₀ имеет вид y - y₀ = y ^' (x - x₀ ), то

y₁ = y₀ + h f (x₀, y₀ ).

Угловой коэффициент в точке P₁ (x₁ ,y₁) также находится из данного дифференциального уравнения y^'₁ = f(x₁,y₁). На следующем шаге получаем новую точку P₂ (x₂ ,y₂), причем

x₂ = x₁ + h, y₂ = y₁ + hf(x₁ ,y₁).

Продолжая вычисления в соответствии с намеченной схемой, получим формулы Эйлера для m приближенных значений решения задачи Коши с начальными данными (x₀, y₀) на сетке отрезка [ a, b ] с шагом h:

x_i = x_{i -1} + h, y_i = y_{i - 1} + hf(x_{i - 1} , y_{i - 1} ) (i = 1,2,..., m) (4)

Графической иллюстрацией приближенного решения является ломаная, соединяющая последовательно точки P₀ , P₁, P₂,...,P_m, которую называют ломаной Эйлера.

Погрешность метода Эйлера можно оценить неравенством

или представить в виде

, где .

Это означает, что метод Эйлера имеет первый порядок точности. В частности, при уменьшении шага h в 10 раз погрешность уменьшится примерно в 10 раз.

Практическую оценку погрешности решения, найденного на сетке с шагом h/2, в точке производят с помощью приближенного равенства - правила Рунге:

, (5)

где p - порядок точности численного метода. Таким образом, оценка полученного результата по формуле (5) вынуждает проводить вычисления дважды: один раз с шагом h, другой - с шагом h/2.

Методы Рунге - Кутта. Численные методы решения задачи Коши

,

на равномерной сетке { x₀ = a, x₁ , x₂ ,...,x_m = b}отрезка[ a, b ]с шагом

h = (b -a ) / m являются методами Рунге - Кутта, если, начиная с данных (x_0,y₀ ), решение ведется по следующим рекуррентным формулам

(6)

Метод называют методом Рунге - Кутта порядка p,если он имеет p - й порядок точности по шагу h на сетке. Порядок точности p достигается с помощью формул (6) при определенных значениях коэффициентов c_j и d_j(j = 1,2,...,p); c₁ всегда полагают равным нулю.

Метод Рунге - Кутта второго порядка называют методом Эйлера - Коши, если p = 2,

c₁ = 0, c₂ = 1, d₁ = d₂ = 1/2. Алгоритм Эйлера - Коши получается из формул (6):

 

x_i =x_i-1 + h, y_i = y_i-1 + Δy_i-1, Δy_i-1 = (1/2)[ k₁^{[i -1]} + k₂^{[i -1]}] (i = 1,..., m), (7)

k₁^{[ i - 1]} = hf (x_i-1,y_i-1), k₂^{[ i - 1 ]} = hf (x_i-1 + h, y_i-1 + hf (x_i-1 ,y_i-1))

 

Для практической оценки погрешности решения можно применять правило Рунге, полагая в формуле (5) р = 2.

Метод Рунге - Кутта четвертого порядка называют классическим методом Рунге - Кутта, если p = 4, c₁ = 0, c₂ = c₃ = 1/2, c₄ = 1, d₁ = d₄ = 1/6, d₂ = d₃ =1/3.

Из рекуррентных формул (6)получим алгоритм решения задачи Коши классическим методом Рунге - Кутта:

 

x _I = x _i - ₁ + h, y _i = y _{i - 1} + Δy _{i – 1} (i = 1,2,..., m),

Δy_i-1 = 1/6 [ k₁^{[ i - 1]} + 2 k₂^{[ i - 1]} + 2k₃^{[ i - 1]} + k₄^{[ i - 1]} ],

k₁^{[ i - 1]} = h f (x_{i - 1,} y_{i -1}),

k₂^{[ i - 1 ]} = h f(x_{i - 1} + (1/2) h, y _{i - 1} + (1/2)k₁^{[ i - 1 ]}),

k₃^{[ i - 1]} = h f (x_i - ₁ + (1/2)h, y _{i - 1} + (1/2)k₂^{[ i - 1 ]}),

k₄^{[ i - 1 ]} = h f(x_{i - 1} + h, y _{i - 1} + k₃^{[ i - 1 ]}),

Графиком приближенного решения является ломаная, последовательно соединяющая точки

P_i(x_i, y_i) ( i = 0, 1, 2,..., m ). С увеличением порядка численного метода звенья ломаной приближаются к ломаной, образованной хордами интегральной кривой y = φ(x), последовательно соединяющими точки (x_i, φ(x_i)) на интегральной кривой.

Правило Рунге (5) практической оценки погрешности решения для численного метода четвертого порядка имеет вид