ВУЗ: Не указан
Категория: Не указан
Дисциплина: Не указана
Добавлен: 25.01.2021
Просмотров: 426
Скачиваний: 1
16
4. ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ ОБЫКНОВЕННЫХ
ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ
Обыкновенным дифференциальным уравнением порядка
k
называется
уравнение
,
0
))
(
),...,
(
'
),
(
,
(
)
(
=
x
y
x
y
x
y
x
F
k
(4.1)
которое связывает независимую переменную
x
, искомую функцию
)
(
x
y
и ее
производные
).
(
),...,
(
'
)
(
x
y
x
y
k
Решение дифференциального уравнения (4.1) заключается в нахождении
функций
)
(
x
y
, которые удовлетворяют этому уравнению для всех значений
x
в
определенном бесконечном или конечном замкнутом или открытом интервале
( )
b
a
,
.
Общее решение обыкновенного дифференциального уравнения порядка
k
имеет вид
),
,...,
,
,
(
2
1
k
C
C
C
x
y
y
=
(4.2)
где
k
C
C
C
,...,
,
2
1
- произвольные постоянные, частный выбор которых дает
частное решение.
Обыкновенные дифференциальные уравнения часто встречаются в
различных прикладных задачах. Однако лишь очень немногие из этих
уравнений могут быть решены точно. В этих случаях используют методы
численного анализа. Следует отметить, что численные методы позволяют
определить лишь частные решения, при этом приближенное значение искомой
функции
)
(
x
y
вычисляется на конечном множестве точек
[ ]
b
a
x
n
,
Î
. Решение
получается в виде таблицы
)
(
n
n
x
y
y
=
. Поэтому такие методы известны как
дискретные.
4.1. Численные методы решения задачи Коши для обыкновенных
дифференциальных уравнений первого порядка
Отыскание частного решения уравнения (4.1), удовлетворяющего
начальным условиям
,
)
(
,...,
'
)
(
'
,
)
(
)
1
(
0
)
1
(
0
0
-
-
=
=
=
k
k
y
a
y
y
a
y
y
a
y
(4.3)
составляет задачу Коши. По
k
начальным условиям (4.3) определяются
k
неизвестных констант
k
C
C
C
,...,
,
2
1
.
Рассмотрим методы численного решения задачи Коши для обыкновенных
дифференциальных уравнений первого порядка.
Пусть требуется найти решение дифференциального уравнения
),
,
(
'
y
x
f
y
=
[ ]
,
,
b
a
x
Î
(4.4)
удовлетворяющее начальному условию
17
.
0
)
(
y
a
y
=
(4.5)
Выберем на отрезке
[ ]
b
a
,
конечное множество точек
n
x
(
n
= 0, 1, ... ,
N
):
b
x
x
x
a
N
=
<
<
<
=
...
1
0
,
в которых и будем искать приближенные значения функции
)
(
n
n
x
y
y
=
.
Узловые точки
n
x
будем считать равноотстоящими (равномерная сетка) с
шагом
n
n
x
x
h
-
=
+
1
.
4.1.1.Разложение решения в ряд Тейлора. Формула Эйлера
Рассмотрим более общую задачу. Предположим, что приближенное
решение в точке
n
x
нам известно. Попытаемся найти его в следующей точке
сетки
1
+
n
x
. Наиболее простой способ его построения основан на разложении в
ряд Тейлора (в предположении дифференцируемости функции
f(x,y)
). Разложим
)
(
)
(
1
h
x
y
x
y
n
n
+
=
+
в ряд Тейлора в окрестности точки
n
x
:
...
!
3
)
(
''
'
!
2
)
(
''
)
(
'
)
(
)
(
3
2
1
+
+
+
+
=
+
h
x
y
h
x
y
h
x
y
x
y
x
y
n
n
n
n
n
,
или
),
,
,
(
)
(
)
(
1
h
y
x
h
x
y
x
y
n
n
n
n
D
+
=
+
где
...
!
2
)
(
''
)
(
'
)
,
,
(
+
+
=
D
h
x
y
x
y
h
y
x
.
Если этот ряд оборвать, оставив в нем
m
слагаемых и заменить
)
(
n
x
y
приближенным значением
n
y
, получим
),
(
)
,
,
(
1
1
+
+
+
+
=
m
n
n
n
n
h
O
h
y
x
h
y
y
j
(
n
= 0, 1, ... ,
N-
1), (4.6)
где
1
)
(
!
)
(
...
!
2
)
(
''
)
(
'
)
,
,
(
-
+
+
+
=
m
m
h
k
x
y
h
x
y
x
y
h
y
x
j
. (4.7)
Входящие в правую часть производные могут быть вычислены:
)
,
(
)
(
'
y
x
f
x
y
=
,
)
,
(
)
,
(
)
,
(
)
,
(
)
,
(
)
,
(
)
(
''
'
'
y
x
f
y
x
f
y
x
f
dx
dy
y
y
x
f
x
y
x
f
dx
y
x
df
x
y
y
x
+
=
¶
¶
+
¶
¶
=
=
и т.д.
Тогда при
m
= 1 имеем
)
(
)
,
(
2
1
h
O
y
x
hf
y
y
n
n
n
n
+
+
=
+
(
формула Эйлера
).
При
m
= 2:
)
(
)]
,
(
)
,
(
)
,
(
(
2
)
,
(
[
3
'
'
1
h
O
y
x
f
y
x
f
y
x
f
h
y
x
f
h
y
y
n
n
n
n
y
n
n
x
n
n
n
n
+
+
+
+
=
+
18
и т.д.
К сожалению, применение этих методов ограничено лишь теми задачами,
для которых легко вычисляются полные производные высокого порядка.
Поэтому на практике широко используется лишь формула Эйлера, хотя ее
точность и не очень высока. Метод Эйлера легко обобщается на случай систем
дифференциальных уравнений 1-го порядка:
),
,...,
,
,
(
2
1
1
'
1
k
y
y
y
x
f
y
=
),
,...,
,
,
(
2
1
2
'
2
k
y
y
y
x
f
y
=
. . . . . . . . . . . . . . . . . . . .
)
,...,
,
,
(
2
1
'
k
k
k
y
y
y
x
f
y
=
с начальными условиями
;
)
(
10
0
1
y
x
y
=
;
)
(
20
0
2
y
x
y
=
... ,
0
0
)
(
k
k
y
x
y
=
.
В этом случае
))
(
),...,
(
),
(
,
(
)
(
)
(
2
1
1
n
k
n
n
n
j
n
j
n
j
x
y
x
y
x
y
x
hf
x
y
x
y
+
=
+
,
(
j
= 1, 2, ... ,
k
;
n
= 0, 1, 2, ... ,
N
-1) .
4.1.2. Метод Рунге-Кутта
Для решения задачи Коши (4.4)-(4.5) воспользуемся формулами (4.6)-(4.7).
Идея метода Рунге-Кутта заключается в построении такого выражения для
)
,
,
(
h
y
x
j
, которое не содержит производных от функции
f
(
x,y
).
Рассмотрим построение формулы Рунге-Кутта для
m
= 2. В этом случае
)]
,
(
)
,
(
)
,
(
(
2
)
,
(
)
,
,
(
'
'
y
x
f
y
x
f
y
x
f
h
y
x
f
h
y
x
y
x
+
+
=
j
. (4.8)
В методе Рунге-Кутта предлагается искать
)
,
,
(
h
y
x
j
в виде:
))
,
(
,
(
)
,
(
)
,
,
(
21
2
2
1
y
x
f
hb
y
ha
x
f
C
y
x
f
C
h
y
x
+
+
+
=
j
, (4.9)
где
21
2
2
1
,
,
,
b
a
C
C
- константы, которые надо определить. Для этого разложим
))
,
(
,
(
21
2
y
x
f
hb
y
ha
x
f
+
+
в формуле (4.9) в ряд Тейлора в окрестности точки
(
x,y
):
=
+
+
+
=
)]
,
(
)
,
(
)
,
(
)
,
(
[
)
,
(
)
,
,
(
21
'
2
'
2
1
y
x
f
hb
y
x
f
ha
y
x
f
y
x
f
C
y
x
f
C
h
y
x
y
x
j
)
,
(
)
,
(
)
,
(
)
,
(
)
(
'
21
2
'
2
2
2
1
y
x
f
y
x
hf
b
C
y
x
hf
a
C
y
x
f
C
C
y
x
+
+
+
=
(4.10)
Сравнивая (4.8) и (4.10), можно получить условия, которым должны
удовлетворять коэффициенты
21
2
2
1
,
,
,
b
a
C
C
:
19
1
2
1
=
+
C
C
,
2
1
2
2
=
a
C
,
2
1
21
2
=
b
C
. (4.11)
Мы получили систему из 3-х уравнений относительно 4-х неизвестных. Эта
система имеет бесконечное множество решений, каждое из которых дает
формулу, имеющую точность порядка
3
h
. Для практического использования
следует выбрать такое решение системы (4.11), при котором формула (4.9)
приобретает наиболее удобный для вычислений вид. Обычно берут
2
1
2
1
=
=
C
C
. Тогда
1
21
2
=
=
b
a
. Подставив эти значения в (4.9) и использовав
полученное выражение в (4.6), приходим к следующей формуле, которая носит
название
двухэтапной схемы Рунге-Кутта
:
)
(
]
[
2
3
2
1
1
h
O
K
K
h
y
y
n
n
+
+
+
=
+
,
где
)
,
(
1
n
n
y
x
f
K
=
,
))
,
(
,
(
2
n
n
n
n
y
x
hf
y
h
x
f
K
+
+
=
.
Чтобы рассчитать
1
+
n
y
по этой формуле, необходимо вычислять функцию
)
,
(
y
x
f
дважды.
Формализм Рунге-Кутта можно обощить, используя
m
вычислений
функции
f(x,y)
на каждом шаге интегрирования. В этом случае получим
m
-
этапный метод Рунге
-
Кутта
:
)
(
)
,
,
(
1
1
+
+
+
+
=
m
n
n
n
n
h
O
h
y
x
h
y
y
j
,
å
=
=
m
r
r
r
K
C
h
y
x
1
)
,
,
(
j
,
)
,
(
1
y
x
f
K
=
,
å
-
=
+
+
=
1
1
)
,
(
r
s
s
rs
r
r
K
b
h
y
ha
x
f
K
,
r
= 2, 3, ... ,
m
.
Наиболее широко используется
четырехэтапная схема метода Рунге-
Кутта
(
m
=4):
)
(
]
2
2
[
6
5
4
3
2
1
1
h
O
K
K
K
K
h
y
y
n
n
+
+
+
+
+
=
+
,
)
,
(
1
n
n
y
x
f
K
=
,
)
2
,
2
(
1
2
K
h
y
h
x
f
K
n
n
+
+
=
,
20
)
2
,
2
(
2
3
K
h
y
h
x
f
K
n
n
+
+
=
,
)
,
(
3
4
hK
y
h
x
f
K
n
n
+
+
=
.
4.1.3. Метод Адамса
Недостатком метода Рунге-Кутта является то, что для получения решения
уравнения в одной точке приходится вычислять правую часть уравнения (4.4) в
нескольких точках. Если правая часть сложна, это приводит к большой
вычислительной работе. Однако существуют методы, применение которых на
каждом шаге требует только однократного вычисления правой части.
Рассмотрим задачу Коши:
)
,
(
'
y
x
f
y
=
,
[ ]
b
a
x
,
Î
,
0
)
(
y
a
y
=
.
Проинтегрируем это уравнение в пределах отрезка
[
]
x
+
x
x
,
, где
[ ]
b
a
x
x
,
,
Î
+
x
:
ò
+
=
-
+
x
x
x
x
dt
t
y
t
f
x
y
x
y
))
(
,
(
)
(
)
(
. (4.12)
Будем считать, что приближенные значения
y(x)
в точках
n
x
x
x
,...,
,
1
0
уже каким-
либо образом вычислены (
n
n
y
x
y
y
x
y
y
x
y
=
=
=
)
(
,...,
)
(
,
)
(
1
1
0
0
). Следовательно, в
этих точках известны и значения производной
)
,
(
'
i
i
i
y
x
f
y
=
(
i
= 0, 1, ... ,
n
).
Заменим в формуле (4.12)
))
(
,
(
t
y
t
f
интерполяционным многочленом, значения
которого совпадают со значениями функции
)
,
(
i
i
y
x
f
в точках
n
x
x
x
,...,
,
1
0
, и
положим в ней
n
x
x
=
и
h
=
x
. Тогда
ò
+
=
-
+
1
)
(
1
n
n
x
x
n
n
n
dx
x
p
y
y
.
Если в качестве
)
(
x
p
n
взять интерполяционный полином Ньютона для
интерполирования назад, то после несложных, но громоздких преобразований
получим:
+
+
-
+
-
+
=
-
-
-
-
+
)
2
(
12
5
)
(
2
'
2
'
1
'
'
1
'
'
1
n
n
n
n
n
n
n
n
y
y
y
h
y
y
h
hy
y
y
...
)
3
3
(
8
3
'
3
'
2
'
1
'
+
-
+
-
+
-
-
-
n
n
n
n
y
y
y
y
h
.