Численные методы решения обыкновенных дифференциальных уравнений
Уравнение
,
связывающее неизвестную функцию y (x), независимую переменную x и производные 
неизвестной функции, называется обыкновенным дифференциальным уравнением. Порядок n старшей производной называется порядком дифференциального уравнения.
В задаче Коши для дифференциального уравнения n-го порядка искомая функция y (x) кроме самого дифференциального уравнения должна удовлетворять начальным условиям
Методы решения задач для дифференциальных уравнений можно разбить на три типа: точные, приближенные и численные [9].
Точными называют методы, с помощью которых решение дифференциального уравнения можно выразить через известные функции (элементарные функции или интегралы от элементарных функций). Точные методы решения известны только для некоторых классов дифференциальных уравнений (линейные дифференциальные уравнения, уравнения с разделяющимися переменными и др.).
Приближенными называются методы, в которых решение находят как предел последовательности функций, являющихся элементарными или интегралами от элементарных функций. Например, метод разложения искомой функции в ряд Тейлора является приближенным методом.
Численный метод решения дифференциального уравнения — это алгоритм вычисления значений искомого решения y (x) на некотором дискретном множестве значений аргумента x. При этом вычисляемые значения искомого решения y (x) являются приближенными, но могут быть и точными.
Задача Коши
Рассмотрим задачу Коши для обыкновенного дифференциального уравнения первого порядка
(7.1)
. (7.2)
Требуется найти функцию 
, которая удовлетворяет уравнению (7.1) на интервале (x 0, X) и начальному условию (7.2) в точке x 0.
Приведем без доказательства теорему существования и единственности решения задачи Коши.
Теорема 7.1. Пусть в области 
функция f (x, y) непрерывна. Тогда на некотором отрезке 
существует решение уравнения (7.1), удовлетворяющее условию (7.2).
Если в области R функция f (x, y) удовлетворяет условию Липшица
,
то указанное решение единственно.
Произведем разбиение отрезка [ x 0, X ] на n частей
(7.3)
Найдем приближенные значения решения y (x) в точках xi, i = 1, 2, …, n.
Метод Эйлера (метод ломаных)
Рассмотрим уравнение (7.1) в точках xi, i = 0, 1, …, n – 1 и заменим производную 
разностной формулой
. (7.4)
Тогда получим рекуррентную формулу метода Эйлера для вычисления приближенных значений y (xi + 1):
(7.5)
Здесь через yi обозначены приближенные значения y (xi), т.е. yi ≈ y (xi),
i = 1, 2, …, n.
На рис. 7.1 дана геометрическая иллюстрация метода Эйлера (7.5). Уравнение касательной к графику решения y (x) в точке (x 0, y 0) имеет вид
(7.6)
так как 
. Интегральная кривая y (x) на отрезке [ x 0, x 1] заменяется отрезком касательной (7.6), соединяющей точку (x 0, y 0) с точкой (x 1, y 1), где 
(см. рис. 7.1). Точка (x 1, y 1) уже не лежит на интегральной кривой y = y (x), удовлетворяющей начальному условию (7.2).
При i = 1 формула (7.5) дает точку (x 2, y 2), которая определяется с помощью касательной 
проведенной в точке (x 1, y 1) к интегральной кривой y (x), удовлетворяющей уравнению (7.1) и начальному условию 
.
Таким образом, с каждым шагом i метод Эйлера (7.5) дает точки (xi, yi), которые, вообще говоря, удаляются от интегральной кривой, соответствующей точному решению задачи Коши (7.1), (7.2). Вместо интегральной кривой метод Эйлера дает ломаную, изображенную на рис. 7.1, поэтому его часто называют методом ломаных.
Рис.7.1.
Формулу (7.5) можно получить и другим способом. Рассмотрим разложение искомого решения y (x) в ряд Тейлора в окрестности точки x 0:
Ограничившись двумя слагаемыми и учитывая, что 
, при x = x 1 получим (7.5).
Мы также получили следующий результат — погрешность вычисления значения y 1 есть величина порядка O (h 2), а погрешность значения yn — величина порядка O (h). Из-за большой погрешности метод Эйлера применяется редко.
Пример 7.1. Решить методом Эйлера на отрезке [1, 2] задачу Коши
Решение в программе Excel. Точным решением данной задачи является функция 
.
Выберем шаг h = 0,1.
В диапазоне A 1: B 12 заполняем ячейки, как показано в таблице 7.1.
В ячейке D 2 записываем начальное значение y (0) = 1, а в ячейке C 2 — формулу =(1+C2^2)*0,5/A2, которая соответствует выражению f (x 0, y 0). Затем выделим ячейку C 2 и маркером заполнения протянем вниз до C 11.
В ячейку D 3 введем формулу = D 2+0,1* C 2 и протянем ячейку D 3 маркером заполнения вниз до D 12.
В столбце E 2: E 12 вычислим значения точного решения. Для этого в ячейку E 2 введем формулу =TAN(LN(КОРЕНЬ(A2))) и протянем ячейку E 2 маркером заполнения вниз до E 12.
В столбце F 2: F 12 вычислим абсолютные ошибки с помощью формулы =ABS(C2-D2), введенной в ячейку F 2 и скопированной в остальные ячейки.
Как видим, значения, полученные методом Эйлера, приближенно равны точным значениям. Абсолютная погрешность каждого значения равна в среднем ε ≈ 0,01. При этом видно, что погрешность yi +1 увеличивается вместе с номером i.
Таблица 7.1
| A | B | C | D | E | F | |
| 1 | i | xi | f(xi,yi) | yi | y(x) | Ошибки |
| 2 | 0 | 1 | 0,5 | 0 | 0 | 0 |
| 3 | 1 | 1,1 | 0,455682 | 0,05 | 0,047691 | 0,002309 |
| 4 | 2 | 1,2 | 0,420472 | 0,095568 | 0,091414 | 0,004154 |
| 5 | 3 | 1,3 | 0,391899 | 0,137615 | 0,13194 | 0,005676 |
| 6 | 4 | 1,4 | 0,368307 | 0,176805 | 0,169842 | 0,006964 |
| 7 | 5 | 1,5 | 0,348547 | 0,213636 | 0,205556 | 0,00808 |
| 8 | 6 | 1,6 | 0,331796 | 0,248491 | 0,239426 | 0,009065 |
| 9 | 7 | 1,7 | 0,317452 | 0,28167 | 0,27172 | 0,00995 |
| 10 | 8 | 1,8 | 0,305064 | 0,313416 | 0,302658 | 0,010758 |
| 11 | 9 | 1,9 | 0,294285 | 0,343922 | 0,332418 | 0,011503 |
| 12 | 10 | 2 |
| 0,37335 | 0,36115 | 0,0122 |
Создадим в программе Excel пользовательские функции для решения задачи Коши (7.1), (7.2) методом Эйлера, используя встроенный язык Visual Basic. Правую часть уравнения (7.1) для нашего примера определяет функция fxy(x, y):
Function fxy(ByVal x, ByVal y): fxy = (1 + y ^ 2) / (2 * x): End Function
Метод Эйлера реализует функция Eiler(x0, y0, x, n):
Function Eiler(ByVal x0, ByVal y0, ByVal x, ByVal n)
h = (x - x0) / n
For i = 1 To n
y1 = y0 + h * fxy(x0, y0): x0 = x0 + h: y0 = y1
Next i: Eiler = y1: End Function
Ключевое слово ByVal в списке параметров указывает на то, что каждый параметр будет передаваться по значению.
Здесь исходными данными являются начальные значения x0, y0, значение переменной x и числовой параметр n, который означает число частей, на которые разбивается отрезок [x0, x], т.е. метод Эйлера применяется с шагом h = (x - x0) / n. Результатом является приближенное значение решения в точке x.
Чтобы ввести эти функции в Excel, выполним команду «Сервис – макрос – Редактор Visual Basic» и в открывшемся окне выбираем пункт меню «Insert – Module» и вводим тексты программ-функций (см. рис. 7.2). Теперь эти функции можно использовать в формулах, как и обычные функции программы Excel.
Переходим в Excel, и построим таблицу решения примера 7.1 с помощью созданных функций. В таблице 7.2 приведены результаты вычислений, а в таблице 7.3 показаны формулы. Уменьшение погрешности по сравнению с таблицей 7.1 объясняется тем, что метод Эйлера в таблице 7.2 применяется с шагом h = 0,01.
Отметим, что правая часть уравнения вызывается в описании функции Eiler. Если понадобится решить задачу Коши с другой правой частью, то достаточно переопределить функцию fxy(x, y).
Таблица 7.2
| A | B | C | D | |
| 1 | xi | yi | y(x) | Ошибки |
| 2 | 1 | 0 | 0 | 0 |
| 3 | 1,1 | 0,047914 | 0,047691 | 0,000223 |
| 4 | 1,2 | 0,091817 | 0,091414 | 0,000403 |
| 5 | 1,3 | 0,132492 | 0,13194 | 0,000552 |
| 6 | 1,4 | 0,17052 | 0,169842 | 0,000678 |
| 7 | 1,5 | 0,206345 | 0,205556 | 0,000788 |
| 8 | 1,6 | 0,240311 | 0,239426 | 0,000885 |
| 9 | 1,7 | 0,272693 | 0,27172 | 0,000973 |
| 10 | 1,8 | 0,30371 | 0,302658 | 0,001052 |
| 11 | 1,9 | 0,333545 | 0,332418 | 0,001126 |
| 12 | 2 | 0,362345 | 0,36115 | 0,001195 |
Таблица 7.3
| A | B | C | D | |
| 1 | xi | yi | y(x) | Ошибки |
| 2 | 1 | 0 | =TAN(LN(КОРЕНЬ(A2))) | =ABS(B2-C2) |
| 3 | 1,1 | =Eiler(A2;B2;A3;10) | =TAN(LN(КОРЕНЬ(A3))) | = ABS (B3-C3) |
| 4 | 1,2 | =Eiler(A3;B3;A4;10) | =TAN(LN(КОРЕНЬ(A4))) | = ABS (B4-C4) |
| 5 | 1,3 | =Eiler(A4;B4;A5;10) | =TAN(LN(КОРЕНЬ(A5))) | = ABS (B5-C5) |
| 6 | 1,4 | =Eiler(A5;B5;A6;10) | =TAN(LN(КОРЕНЬ(A6))) | = ABS (B6-C6) |
| 7 | 1,5 | =Eiler(A6;B6;A7;10) | =TAN(LN(КОРЕНЬ(A7))) | = ABS (B7-C7) |
| 8 | 1,6 | =Eiler(A7;B7;A8;10) | =TAN(LN(КОРЕНЬ(A8))) | = ABS (B8-C8) |
| 9 | 1,7 | =Eiler(A8;B8;A9;10) | =TAN(LN(КОРЕНЬ(A9))) | = ABS (B9-C9) |
| 10 | 1,8 | =Eiler(A9;B9;A10;10) | =TAN(LN(КОРЕНЬ(A10))) | = ABS (B10-C10) |
| 11 | 1,9 | =Eiler(A10;B10;A11;10) | =TAN(LN(КОРЕНЬ(A11))) | = ABS (B11-C11) |
| 12 | 2 | =Eiler(A11;B11;A12;10) | =TAN(LN(КОРЕНЬ(A12))) | = ABS (B12-C12) |
Рис.7.2. Окно редактора Visual Basic
Создадим в программе Mathcad функции для решения задачи Коши методом Эйлера и вычислим решение в точке x = 2 с шагом h = 0,1 (n = 10) и шагом h = 0,01 (n = 100):
|
|
|
|
|
|
|
|
Как видим, результаты расчетов в программе Mathcad совпадают со значениями, полученными в Excel.
Метод Эйлера с уточнением
Для повышения точности метода Эйлера применяют следующий прием. Сначала находят приближенное значение решения по методу Эйлера
,
а затем уточняют его по формуле
Этот метод называется методом Эйлера-Коши, и рекуррентные соотношения для его реализации могут быть записаны в виде:
(7.7)
Метод Эйлера-Коши имеет погрешность порядка O (h 2). Геометрическая иллюстрация метода Эйлера-Коши дается на рис. 7.2. Очередное приближение метода Эйлера-Коши соответствует точке пересечения диагоналей параллелограмма, построенного на двух звеньях ломаной метода Эйлера.
Метод Эйлера-Коши является одним из частных случаев методов Рунге-Кутта, которые рассмотрены в следующем пункте.
Рис.7.2
Методы Рунге-Кутта
Рассмотрим наиболее популярный метод решения задачи Коши — метод Рунге-Кутта. Этот метод позволяет строить формулы расчета приближенного решения практически любого порядка точности.
Выведем формулы метода Рунге-Кутта второго порядка точности. Для этого решение представим куском ряда Тейлора, отбрасывая члены с порядком выше второго. Тогда приближенное значение искомой функции в точке x 1 можно записать в виде:
. (7.8)
Вторую производную 
можно выразить через производную функции f (x, y), однако в методе Рунге-Кутта вместо производной используют разность
,
соответственно подбирая значения параметров 
. Тогда (7.8) можно переписать в виде
, (7.9)
где α, β, γ и δ — некоторые параметры. Рассматривая правую часть (7.9) как функцию аргумента h, разложим её по степеням h:
,
и выберем параметры α, β, γ и δ так, чтобы это разложение было близко к (7.8). Отсюда следует, что
α + β = 1, αγ = 0,5, αδ = 0,5 f (x 0, y 0).
С помощью этих трех уравнений выразим β, γ и δ через параметр α, тогда получим
(7.10)
Теперь, если вместо (x 0, y 0) в (7.10) подставить (x 1, y 1), то получим формулу для вычисления y 2 — приближенного значения искомой функции в точке x 2.
В общем случае метод Рунге-Кутта применяется на произвольном разбиении отрезка [ x 0, X ] на n частей, т.е. с переменным шагом
(7.11)
Параметр α выбирают равным 1 или 0,5. Запишем окончательно расчетные формулы метода Рунге-Кутта второго порядка с переменным шагом для α = 1:
(7.12)
и α = 0,5:
(7.13)
Дадим геометрическую интерпретацию методов Рунге-Кутта (7.12), (7.13).
На рисунке 7.3 а) иллюстрируется формула (7.12). Сначала по методу Эйлера вычисляется приближенное значение в точке xi + hi /2. В этой точке определяется касательная к интегральной кривой, параллельно которой через точку (xi, yi) проводится прямая до точки пересечения с прямой x = xi + 1. Ордината точки пересечения yi + 1 принимается за приближенное значение искомого решения в точке xi + 1.
Рисунок 7.3 б) интерпретирует формулу (7.13). Методом Эйлера вычисляется приближенное значение 
в точке xi +1. В этой точке тангенс угла наклона касательной к интегральной кривой равен выражению 
. Через точку (xi, yi) проводится прямая, тангенс угла наклона которой определяется как полусумма угловых коэффициентов касательных, проведенных через точки (xi, yi) и 
. За приближенное значение искомого решения в точке xi +1 принимается ордината точки пересечения этой прямой с прямой x = xi +1.
Отметим, что метод (7.13) совпадает с методом Эйлера-Коши (7.7).
Рис.7.3.
Теорема 7.1. Если правая часть f (x, y) уравнения (7.1) непрерывна и ограничена вместе со своими производными до второго порядка включительно, то решение, полученное по формулам (7.12), (7.13), равномерно сходится к решению задачи (7.1), (7.2) с погрешностью 
.
Приведем наиболее употребительные формулы метода Рунге-Кутта, формулы четвертого порядка точности:
(7.14)
Пример 7.2. Решить методом Рунге-Кутта (7.12) задачу Коши из примера 7.1:
Решение в программе Excel. В таблице 7.4 приведены результаты применения метода Рунге-Кутта второго порядка (7.12) с шагом h = 0,1.
Сравнивая таблицы 7.2 и 7.4 можно увидеть, что метод Рунге-Кутта второго порядка (7.12) с шагом h = 0,1 дает лучшие результаты, чем метод Эйлера с шагом 0,01.
Таблица 7.4
| A | B | C | D | E | F | |
| 1 | xi | f(xi,yi) | f(xi+0,5h, yi+0,5 f(xi,yi)) | yi | y(x) | Ошибки |
| 2 | 1 | 0,5 | 0,476488 | 0 | 0 | 0 |
| 3 | 1,1 | 0,455577 | 0,436939 | 0,047649 | 0,047691 | 4,24E-05 |
| 4 | 1,2 | 0,420143 | 0,405049 | 0,091343 | 0,091414 | 7,14E-05 |
| 5 | 1,3 | 0,391301 | 0,378861 | 0,131848 | 0,13194 | 9,22E-05 |
| 6 | 1,4 | 0,367432 | 0,357029 | 0,169734 | 0,169842 | 0,000108 |
| 7 | 1,5 | 0,347401 | 0,338594 | 0,205437 | 0,205556 | 0,00012 |
| 8 | 1,6 | 0,330395 | 0,322861 | 0,239296 | 0,239426 | 0,00013 |
| 9 | 1,7 | 0,315811 | 0,309309 | 0,271582 | 0,27172 | 0,000138 |
| 10 | 1,8 | 0,303198 | 0,297545 | 0,302513 | 0,302658 | 0,000145 |
| 11 | 1,9 | 0,292211 | 0,287263 | 0,332268 | 0,332418 | 0,000151 |
| 12 | 2 |
|
| 0,360994 | 0,36115 | 0,000156 |
Создадим функцию для вычисления решения задачи Коши (7.1), (7.2) методом Рунге-Кутта четвертого порядка (7.14):
Function Runge_Kutta4(ByVal x0, ByVal y0, ByVal x, ByVal n)
h = (x - x0) / n
For i = 1 To n
k1 = fxy(x0, y0)
k2 = fxy(x0 + h / 2, y0 + h * k1 / 2)
k3 = fxy(x0 + h / 2, y0 + h * k2 / 2)
k4 = fxy(x0 + h, y0 + h * k3)
y1 = y0 + h * (k1 + 2 * k2 + 2 * k3 + k4) / 6
x0 = x0 + h: y0 = y1
Next i
Runge_Kutta4 = y0
End Function
В таблице 7.5 приведены результаты решения задачи Коши методом Рунге-Кутта четвертого порядка (7.14) с параметром n = 10, т.е. на каждом этапе вычисления yi применялся шаг h = 0,01. Поэтому абсолютные ошибки очень малы.
Таблица 7.5
| xi | yi | y(x) | Ошибки |
| 1 | 0,000000000000000 | 0,000000000000000 | 0,000000000000000 |
| 1,1 | 0,047691197731806 | 0,047691197726655 | 0,000000000005151 |
| 1,2 | 0,091414144750546 | 0,091414144742135 | 0,000000000008412 |
| 1,3 | 0,131939841952911 | 0,131939841942306 | 0,000000000010605 |
| 1,4 | 0,169841513601824 | 0,169841513589657 | 0,000000000012167 |
| 1,5 | 0,205556457698103 | 0,205556457684762 | 0,000000000013341 |
| 1,6 | 0,239425622368332 | 0,239425622354065 | 0,000000000014267 |
| 1,7 | 0,271719843610883 | 0,271719843595851 | 0,000000000015032 |
| 1,8 | 0,302657774396418 | 0,302657774380728 | 0,000000000015690 |
| 1,9 | 0,332418460630980 | 0,332418460614704 | 0,000000000016276 |
| 2 | 0,361150365759415 | 0,361150365742600 | 0,000000000016814 |
Создадим в программе Mathcad функции для решения задачи Коши методом Рунге-Кутта и вычислим решение в точке x = 2 с шагом h = 0,01 (n = 100):
|
|
|
|
