дифф длинные

1. 7. Численные методы решения обыкновенных дифференциальных уравнений

Численные методы решения
обыкновенных
дифференциальных уравнений
7.
Уравнение F x, y, y , y ,..., y n 0 связывающее неизвестную функцию y(x),
независимую переменную x и производные y x , y x , ..., y n x неизвестной
функции, называется обыкновенным дифференциальным уравнением.
Порядок n старшей производной называется порядком дифференциального
уравнения.
В задаче Коши для дифференциального уравнения n-го порядка искомая
функция y(x) кроме самого дифференциального уравнения должна
удовлетворять начальным условиям
y x0 0 , y x0 1 , y x0 2 , ..., y n 1 x0 n 1.
Методы решения задач для дифференциальных уравнений можно разбить на
три типа: точные, приближенные и численные.

2.

Точными называют методы, с помощью которых решение
дифференциального уравнения можно выразить через известные функции
(элементарные функции или интегралы от элементарных функций). Точные
методы решения известны только для некоторых классов
дифференциальных уравнений (линейные дифференциальные уравнения,
уравнения с разделяющимися переменными и др.).
Приближенными называются методы, в которых решение находят как
предел последовательности функций, являющихся элементарными или
интегралами от элементарных функций. Например, метод разложения
искомой функции в ряд Тейлора является приближенным методом.
Численный метод решения дифференциального уравнения — это алгоритм
вычисления значений искомого решения y(x) на некотором дискретном
множестве значений аргумента x. При этом вычисляемые значения искомого
решения y(x) являются приближенными, но могут быть и точными.

3. 7.1. Задача Коши

Рассмотрим задачу Коши для обыкновенного дифференциального уравнения
первого порядка
(7.1)
y x f x, y ,
y x0 y 0
(7.2)
Требуется найти функцию y y x , которая удовлетворяет уравнению (7.1)
на интервале (x0, X) и начальному условию (7.2) в точке x0.
Приведем без доказательства теорему существования и единственности
решения задачи Коши.
Теорема 7.1. Пусть в области R x, y , x x0 a, y y0 b функция f(x, y)
непрерывна. Тогда на некотором отрезке x x0 d существует решение
уравнения (7.1), удовлетворяющее условию (7.2).
Если в области R функция f(x, y) удовлетворяет условию Липшица
f x, y1 f x, y2 k y1 y2 , k 0, k const

4.

то указанное решение единственно.
Произведем разбиение отрезка [x0, X] на n частей
xi x0 ih, h
X x0
.
n
(7.3)
Найдем приближенные значения решения y(x) в точках xi, i = 1, 2, …, n.
7.1.1. Метод Эйлера (метод ломаных)
Рассмотрим уравнение (7.1) в точках xi, i = 0, 1, …, n – 1 и заменим y xi
производную разностной формулой
y xi
y xi 1 y xi
h
(7.4)
Тогда получим рекуррентную формулу метода Эйлера для вычисления
приближенных значений y(xi + 1):
yi 1 yi hf xi , yi , i 0, 1, ..., n 1. (7.5)

5.

Здесь через yi обозначены приближенные значения y(xi), т.е. yi ≈ y(xi),
i = 1, 2, …, n.
На рис. 7.1 дана геометрическая иллюстрация метода Эйлера (7.5).
Уравнение касательной к графику решения y(x) в точке (x0, y0) имеет вид
(7.6)
y y f x , y x x ,
0
0
0
0
так как . Интегральная кривая y(x) на отрезке [x0, x1] заменяется отрезком
касательной (7.6), соединяющей точку (x0, y0) с точкой (x1, y1), где
y1 y0 f x0 , y0 x1 x0 (см. рис. 7.1). Точка (x1, y1) уже не лежит на
интегральной кривой y = y(x), удовлетворяющей начальному условию (7.2).
При i = 1 формула (7.5) дает точку (x2, y2), которая определяется с помощью
касательной y y1 f x1 , y1 x x1 , проведенной в точке (x1, y1) к
интегральной кривой y(x), удовлетворяющей уравнению (7.1) и начальному
условию y x1 y1.
Таким образом, с каждым шагом i метод Эйлера (7.5) дает точки (xi, yi),
которые, вообще говоря, удаляются от интегральной кривой, соответствующей
точному решению задачи Коши (7.1), (7.2). Вместо интегральной кривой метод
Эйлера дает ломаную, изображенную на рис. 7.1, поэтому его часто называют
методом ломаных.

6.

y3
y2
y1
y0
x0
x1
x2
x3
x
Рис.7.1.
Формулу (7.5) можно получить и другим способом. Рассмотрим разложение
искомого решения y(x) в ряд Тейлора в окрестности точки x0:
y x y x0 y x0 x x0
y x0
x x0 2 y x0 x x0 3 ...
2!
3!
Ограничившись двумя слагаемыми и учитывая, что y x0 f x0 , y0 , при x = x1
получим (7.5).
Мы также получили следующий результат — погрешность вычисления
значения y1 есть величина порядка O(h2), а погрешность значения yn —
величина порядка O(h). Из-за большой погрешности метод Эйлера
применяется редко.

7. 7.1.2. Метод Эйлера с уточнением

Для повышения точности метода Эйлера применяют следующий прием.
Сначала находят приближенное значение решения по методу Эйлера
~
yi 1 yi hf xi , yi
а затем уточняют его по формуле
f xi , yi f xi 1 , ~
yi 1
yi 1 yi h
.
2
Этот метод называется методом Эйлера-Коши, и рекуррентные соотношения
для его реализации могут быть записаны в виде:
yi 1 yi h
f xi , yi f xi h, yi hf xi , yi
, i 0, 1, ..., n 1. (7.7)
2
Метод Эйлера-Коши имеет погрешность порядка O(h2). Геометрическая
иллюстрация метода Эйлера-Коши дается на рис. 7.2. Очередное
приближение метода Эйлера-Коши соответствует точке пересечения
диагоналей параллелограмма, построенного на двух звеньях ломаной метода
Эйлера.

8.

Метод Эйлера-Коши является одним из частных случаев методов Рунге-Кутта,
которые рассмотрены в следующем пункте.
yi+1
yi
x
xi
xi+1
Рис. 7.2
xi+1

9. 7.1.3. Методы Рунге-Кутта

Рассмотрим наиболее популярный метод решения задачи Коши — метод
Рунге-Кутта. Этот метод позволяет строить формулы расчета приближенного
решения практически любого порядка точности.
Выведем формулы метода Рунге-Кутта второго порядка точности. Для этого
решение представим куском ряда Тейлора, отбрасывая члены с порядком
выше второго. Тогда приближенное значение искомой функции в точке x1
можно записать в виде:
y x0
y x0 2
2
x1 x0 y0 f x0 , y0 h
y1 y x0 y x0 x1 x0
h (7.8)
2
2
Вторую производную можно выразить через производную функции f(x, y),
однако в методе Рунге-Кутта вместо производной используют разность
d
f ~
x, ~
y f x, y
y x
f x, y
dx
x
соответственно подбирая значения параметров . Тогда (7.8) можно
переписать в виде
y1 y0 h βf x0 , y0 αf x0 γh, y0 δh
(7.9)

10.

где α, β, γ и δ — некоторые параметры. Рассматривая правую часть (7.9) как
функцию аргумента h, разложим её по степеням h:
y1 y0 α β hf x0 , y0 αh 2 γf x x0 , y0 δf y x0 , y0
и выберем параметры α, β, γ и δ так, чтобы это разложение было близко к
(7.8). Отсюда следует, что
α + β = 1, αγ = 0,5, αδ = 0,5 f(x0, y0).
С помощью этих трех уравнений выразим β, γ и δ через параметр α, тогда
получим
h
h
y1 y0 h 1 - α f x0 , y0 αf x0
, y0
f x0 , y 0 ,
2α
2α
(7.10)
0 α 1.
Теперь, если вместо (x0, y0) в (7.10) подставить (x1, y1), то получим формулу
для вычисления y2 — приближенного значения искомой функции в точке x2.
В общем случае метод Рунге-Кутта применяется на произвольном разбиении
отрезка [x0, X] на n частей, т.е. с переменным шагом
(7.11)
x , x , ... x ; h x x , x X .
0
1
n
i
i 1
i
n

11.

Параметр α выбирают равным 1 или 0,5. Запишем окончательно расчетные
формулы метода Рунге-Кутта второго порядка с переменным шагом для α = 1:
h
h
yi 1 yi hi f xi i , yi i f xi , yi , i 0, 1, ..., n 1.
2
2
(7.12)
и α = 0,5:
(7.13)
hi
yi 1 yi f xi , yi f xi hi , yi hi f xi , yi , i 0, 1, ..., n 1.
2
Дадим геометрическую интерпретацию методов Рунге-Кутта (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. В этой точке тангенс угла наклона
касательной к интегральной кривой равен выражению

12.

f xi hi , yi hi f xi , yi . Через точку (xi, yi) проводится прямая, тангенс угла
наклона которой определяется как полусумма угловых коэффициентов
касательных, проведенных через точки (xi, yi) и xi 1 , yi hi f xi , yi . За
приближенное значение искомого решения в точке xi+1 принимается ордината
точки пересечения этой прямой с прямой x = xi+1.
Отметим, что метод (7.13) совпадает с методом Эйлера-Коши (7.7).
y
y
y(xi+1)
y(xi+1)
yi+1
yi+1
yi
x
xi
xi + hi/2
xi+1
yi
x
xi
а)
xi+1
б)
Рис.7.3.

13.

Теорема 7.1. Если правая часть f(x, y) уравнения (7.1) непрерывна и
ограничена вместе со своими производными до второго порядка
включительно, то решение, полученное по формулам (7.12), (7.13),
равномерно сходится к решению задачи (7.1), (7.2) с погрешностью O max hi2 .
Приведем наиболее употребительные формулы метода Рунге-Кутта,
формулы четвертого порядка точности:
h
yi 1 yi k1 2k 2 2k 3 k 4 ,
6
h
h
k1 f xi , yi , k 2 f xi , yi k1 ,
2
2
h
h
k 3 f xi , yi k 2 , k 4 f xi h, yi hk 3 .
2
2
(7.14)

14. 7.1.4. Метод Рунге-Кутта с автоматическим выбором шага. Правило Рунге оценки погрешности

Для метода Рунге-Кутта применимо правило Рунге для оценки погрешности.
Пусть y x; h приближенное значение решения в точке x, полученное по
формулам (7.12), (7.13) или (7.14) с шагом h, а p — порядок точности
соответствующей формулы. Тогда погрешность R h значения y x; h можно
оценить, используя приближенное значение y x;2h решения в точке x,
полученное с шагом 2h
R h
y x; h y x;2h
2 p 1
(7.15)
где p = 2 для формул (7.12), (7.13) и p = 4 для (7.14). Уточненное решение
запишется в виде
y x; h y x;2h
~
y x; h y x; h
(7.16)
2p 1

15.

В алгоритмах с автоматическим выбором шага предварительно задают
погрешность в виде положительного параметра ε, и на каждом этапе
вычисления следующего значения подбирают шаг h такой, что выполняется
неравенство
y x; h y x;2h
(7.17)
ε
p
2 1
7.1.5. Задача Коши для системы
дифференциальных уравнений
Метод Рунге-Кутта применим и к задаче Коши для системы m
дифференциальных уравнений с m неизвестными функциями
y1 f1 x, y1 ,..., ym ,
y f x, y ,..., y ,
2
2
1
m
...
ym f m x, y1 ,..., ym ,
x x0 , X
(7.18)

16.

y1 x0 y1, 0 ,
y2 x0 y2, 0 , ...,
ym x0 ym, 0 .
(7.19)
Приведем для задачи (7.18), (7.19) расчетные формулы метода Рунге-Кутта
четвёртого порядка. Пусть требуется найти систему m функций
удовлетворяющих в интервале (x0, X) дифференциальным уравнениям (7.18),
а в точке x0 — начальным условиям (7.19). Предположим, что отрезок [x0, X]
разбит на N частей
xi x0 ih , h
X x0
.
N
Тогда каждую l-ую функцию yl(x) можно приближенно вычислять в точках xi+1
по формулам Рунге-Кутта
K l ,1 f l xi , y1,i , y2,i , ..., ym ,i , l 1, 2, ..., m;
h
h
h
h
K l , 2 f l xi , y1,i K1,1 , y2,i K 2,1 , ..., ym ,i K m,1 , l 1, 2, ..., m;
2
2
2
2
h
h
h
h
(7.20)
K l ,3 f l xi , y1,i K1, 2 , y2,i K 2, 2 , ..., ym ,i K m , 2 , l 1, 2, ..., m;
2
2
2
2
K l , 4 f l xi h, y1,i hK1,3 , y2,i hK 2,3 , ..., ym,i hK m ,3 , l 1, 2, ..., m;
yl ,i 1 yl ,i
h
K l ,1 2 K l ,2 2 K l ,3 K l ,4 , l 1, 2, ..., m.
6

17.

Здесь через yl,i мы обозначили приближенное значение функции yl(x) в точке xi.
Обращаем внимание на порядок вычислений по формулам (7.20). На каждом
шаге сначала вычисляются коэффициенты Kl,i в следующем порядке
K1,1, K2,1, …, Km,1,
K1,2, K2,2, …, Km,2,
K1,3, K2,3, …, Km,3,
K1,4, K2,4, …, Km,4,
и лишь затем приближенные значения функций y1,i+1, y2,i+1, …, ym,i+1.
7.1.6. Задача Коши для дифференциальных
уравнений высших порядков
Задача Коши для дифференциального уравнения n-го порядка
y ( n ) f x, y, y ..., y ( n 1) ,
x x0 , X
(7.21)
y x0 y0 , y x0 y1,0 , ..., y ( n 1) x0 yn 1,0 (7.22)

18.

легко сводится к задаче Коши для системы дифференциальных уравнений с
помощью замены переменных
z0 y, z1 y , ..., z n 1 y ( n 1) . (7.23)
Учитывая (7.23) из уравнения (7.21) получим систему дифференциальных
уравнений
z 0 z1 ,
z1 z 2 ,
...
z z
n 1,
n 2
z n 1 f x, z 0 , z1 ,..., z n 1 .
(7.24)
Начальные условия (7.22) для функций zl переписываются в виде
z0 x0 y0 ,
z1 x0 y1, 0 , ...,
z n 1 x0 yn 1, 0
(7.25)
Запишем для полученной системы метод Рунге-Кутта
z l ,i 1 z l ,i
h
K l ,1 2K l ,2 2 K l ,3 K l ,4 , i 0, 1, ..., N ; l 0, 1, ..., n 1 (7.26)
6

19.

Для вычисления коэффициентов Kl,1, Kl,2, Kl,3 и Kl,4 имеем следующие
формулы:
K 0,1 z1,i ,
K 1,1 z 2,i ,
...
K n 1,1 f xi , z 0,i , z1,i ,..., z n 1,i ,
h
K 1,1 ,
2
h
K1, 2 z 2,i K 2,1 ,
2
...
K 0, 2 z1,i
h
h
h
h
K n 1, 2 f xi , z 0,i K 0,1 , z1,i K 1,1 ,..., z n 1,i K n 1,1 ,
2
2
2
2

20.

h
K1, 2 ,
2
h
K1,3 z 2,i K 2, 2 ,
2
...
K 0,3 z1,i
h
h
h
h
K n 1,3 f xi , z 0,i K 0, 2 , z1,i K 1, 2 ,..., z n 1,i K n 1, 2 ,
2
2
2
2
K 0, 4 z1,i hK1,3 ,
K1, 4 z 2,i hK 2,3 ,
...
K n 1, 4 f xi h, z0,i hK 0, 2 , z1,i hK1, 2 ,..., z n 1,i hK n 1, 2 .

21. 7.2. Краевые задачи для обыкновенных дифференциальных уравнений

Примерами краевых задач для обыкновенных дифференциальных уравнений
являются следующие задачи:
1) Найти функцию u(x), которая удовлетворяет на отрезке [a, b] уравнению
(7.27)
u x f x
а на концах отрезка условиям
u(a) = u(b) = 0.
(7.28)
Задача (7.27), (7.28) имеет следующее содержание. Между точками
x = a и x = b натянута упругая струна, находящаяся под действием внешней
изгибающей нагрузки. f(x) — величина нагрузки, а u(x) — прогиб струны в
безразмерных единицах.
2) Двухточечная краевая задача для линейного дифференциального
уравнения второго порядка
u x p x u x q x u x f x
(7.29)
0u a 1u a A,
0u b 1u b B.
(7.30)

22. 7.2.1. Метод прогонки

Рассмотрим частный случай задачи (7.29), (7.30)
u x p x u x f x , x 0, X , (7.31)
u 0 a, u X b.
(7.32)
Введем сетку: xi ih, i 0, 1, ..., N ; h X / N . Обозначим через ui
приближенные значения u(x) в узлах сетки. Рассмотрим уравнение (7.31) во
внутренних узлах и заменим производную второго порядка разностной
формулой (5.11):
u xi
ui 1 2ui ui 1
h2
(7.33)
Тогда из (7.31), (7.32) получим для определения ui систему линейных
уравнений
ui 1 2ui ui 1
pi ui f i , i 1, 2, ..., N 1,
2
h
(7.34)
u0 a, u N b.

23.

Система (7.34) при p(x) ≥ 0 имеет решение.
Система (7.34) представляет собой систему линейных уравнений с
трехдиагональной матрицей
2 p1h 2 u1 u2 h 2 f1 a, i 1,
2
2
ui 1 2 pi h ui ui 1 h f i , i 2, 3,..., N 2,
2
2
u
2
p
h
u
h
f N 1 b, i N 1,
N 1
N 1
N 2
(7.35)
и для её решения применим метод прогонки, который фактически является
методом исключения неизвестных Гаусса.
1. Прямой ход прогонки. Запишем первое уравнение (7.35) в виде
1
a h 2 f1
(7.36)
u1 α1u2 β1 , α1
2 p h , β 2 p h .
2
1
2
1
1
Подставив (7.36) во второе уравнение (7.35) и упростив выражения, увидим,
что можно для ui получить формулы
β i 1 h 2 f i
1
ui α i ui 1 β i , α i
, βi
, i 2,..., N 2.
2
2
2 pi h α i 1
2 pi h α i 1
(7.37)

24.

Из последнего уравнения (7.35), учитывая (7.37) при i N 2 , получим
b β N 2 h 2 f N 1
u N 1 β N 1 , α N 1 0, β N 1
.
2
2 p N 1h α N 2
(7.38)
2. Обратный ход прогонки. После вычисления прогоночных коэффициентов
можно найти значения решения задачи. Формула (7.38) дает значение u N 1 :
u N 1 β N 1.
(7.39)
Остальные значения вычисляем в обратном порядке
ui α i ui 1 βi , i N 2, N 3, ..., 2, 1.
(7.40)

25. 7.2.2. Метод стрельбы (метод пристрелки)

Рассмотрим метод пристрелки на примере решения линейной краевой задачи
(7.31), (7.32)
u x p x u x f x , x 0, X ,
u 0 a, u X b.
Если известны частное решение u x неоднородного уравнения u x p x u x f x
и два линейно независимых решения u1 x , u2 x однородного уравнения
u x p x u x 0 , то общее решение неоднородного дифференциального
уравнения u x p x u x f x можно записать в виде
u x C1u1 x C2u2 x .
Постоянные C1, C2 можно определить из краевых условий (7.32).
В методе пристрелки используется следующий способ.
Сначала находят частное решение u x неоднородного уравнения u x p x u x f x
удовлетворяющее условию u 0 a , и частное решение u1 x однородного
уравнения u x p x u x 0, удовлетворяющее условию u 0 0 .
Затем общее решение u x неоднородного уравнения u x p x u x f x ,
удовлетворяющее условию u 0 a записывают в виде u x Cu1 x . И тогда
остается найти найти постоянную C из условия u X Cu1 X b .

26.

Приведем сеточный аналог метода пристрелки. Пусть краевая задача
приведена к системе уравнений (7.34)
ui 1 2ui ui 1
pi ui f i , i 1, 2, ..., N 1,
2
h
u0 a, u N b.
(7.41)
Найдем частные решения неоднородной (7.42) и однородной (7.43) систем
уравнений:
ui0 1 2ui0 ui0 1
pi ui0 f i , i 1, 2, ..., N 1,
2
h
(7.42)
u 0 a.
0
ui1 1 2ui1 ui1 1
1
p
u
0, i 1, 2, ..., N 1,
i
i
2
h
u1 0.
0
(7.43)
Выбирая произвольные значения для u1 , u1 (при этом должно быть u1 0 ),
находим из (7.42) и (7.43) формулы для вычисления частных решений:
0
1
0
0
u0 a, u1 a h,
0
2
0
0
2
ui 1 2 h pi ui ui 1 h f i , i 1, 2, ..., N 1.
1
(7.44)

27.

1
1
u0 0, u1 h,
1
2
1
1
u
2
h
p
u
u
i 1, 2, ..., N 1.
i
1
i
i
i
1i ,
(7.45)
Найдем C из условия u N C u N b и запишем решение:
0
1
b u N0
C
u1N
ui ui0 C ui1 , i 0, 1, ..., N .
(7.46)
(7.47)
Алгоритм метода пристрелки заключается в том, чтобы выбрать шаг h и
выполнить последовательно вычисления по формулам (7.44) — (7.47).