Схема Эйлера.
Схема Эйлера является простейшей одношаговой разностной схемой. В ней для расчета интеграла I
i,i+1
используется только значение функции f(t, u) в начале отрезка [t
i
, t
i+1
]:
I
i,i+1
≈ ∆t
i
·f(t
i
, u
i
). (1.3)
Здесь ∆t
i
= (t
i+1
- t
i
) – шаг изменения независимой переменной.
Подставляя выражение (1.3) в уравнение (1.2), получим разно- стную схему Эйлера:
)
,
(
1
i
i
i
i
i
u
t
f
t
u
u
⋅
∆
+
=
+
. (1.4)
Определение интеграла по значению функции в одной точке от- резка является достаточно грубым, поэтому схема Эйлера имеет только первый порядок аппроксимации.
Исправленный метод Эйлера.
Точность расчета интеграла I
i,i+1
можно повысить, добавив еще одну расчетную точку. Сначала по- лучим с помощью схемы Эйлера приближенное значение решения
)
,
(
1
i
i
i
i
i
u
t
f
t
u
u
⋅
∆
+
=
+
в конечной точке рассматриваемого отрез- ка. Используя это значение, получим приближенное значение
12 функции
)
,
(
1 1
+
+
i
i
u
t
f
в конце отрезка. Приближенное значение интеграла получается по формуле трапеций:
)))
,
(
,
(
)
,
(
(
5 0
))
,
(
)
,
(
(
5 0
1 1
1 1
,
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
u
t
f
t
u
t
f
u
t
f
t
u
t
f
u
t
f
t
I
∆
+
+
∆
=
+
∆
≈
+
+
+
+
. (1.5)
Подставляя выражение (1.5) в уравнение (1.2), получим разно- стную схему исправленного метода Эйлера:
)))
,
(
,
(
)
,
(
(
5 0
1 1
i
i
i
i
i
i
i
i
i
i
u
t
f
t
u
t
f
u
t
f
t
u
u
∆
+
+
∆
+
=
+
+
. (1.6)
Схема имеет второй порядок аппроксимации.
Модифицированный метод Эйлера.
После определения по схеме
Эйлера приближенного значения решения
1
+
i
u
для расчета инте- грала можно использовать метод прямоугольников:
)))
,
(
5 0
,
5 0
(
(
)
2
/
)
(
,
2
/
)
((
(
1 1
1
,
i
i
i
i
i
i
i
i
i
i
i
i
i
i
u
t
f
t
u
t
t
f
t
u
u
t
t
f
t
I
∆
+
∆
+
∆
=
+
+
∆
≈
+
+
+
. (1.7)
Подставляя выражение (1.7) в уравнение (1.2), получим разно- стную схему модифицированного метода Эйлера:
)))
,
(
5 0
,
5 0
(
(
1
i
i
i
i
i
i
i
i
i
u
t
f
t
u
t
t
f
t
u
u
∆
+
∆
+
∆
+
=
+
. (1.8)
Как и для исправленного метода Эйлера, схема модифициро- ванного метода Эйлера имеет второй порядок аппроксимации.
Метод Рунге-Кутта.
В общем случае интеграл I
i,i+1
определяется с помощью n точек на интервале [t
i
, t
i+1
]. Первая точка совпадает t
i
, а остальные располагаются оптимальным образом с точки зрения получения при заданном n наивысшего порядка аппроксимации.
Наиболее широко используется схема с 4-мя точками, в которой расчет интеграла I
i,i+1
производится по формуле Симпсона:
13 6
/
)
2 2
(
4 3
2 1
1
fffftuuiii+
+
+
∆
+
=
+
, (1.9) где
f1
=
f(
ti
,
ui
),
f2
=
f(
ti
+∆t i
/2,
ui
+ ∆t i
f1
/2),
f3
=
f(
ti
+∆t i
/2,
ui
+ ∆t i
f2
/2),
f4
=
f(
ti
+∆t i
,
ui
+ ∆t i
f3
).
Эта схема называется схемой Рунге-Кутта и имеет 4-й порядок аппроксимации.
Все рассмотренные схемы являются явными и обладают услов- ной устойчивостью.
1.1.3. Многошаговые методы решения задачи Коши для ОДУ Для оценки интеграла
Ii,i+1 можно использовать информацию об изменении функции
f(t, u) в предыдущие моменты времени. На ос- нове значений
f(t, u) в точках
ti,
ti-1
, …
ti-k+1
построим интерполяци- онный полином, совпадающий в этих точках с
f(ti, ui),
f(ti-1
, ui-1
), …
f(ti-k+1
, ui-k+1
). Проэкстраполируем полином на отрезок [
ti, ti+1] и най- дем аналитическое выражение для интеграла
Ii,i+1. Подставив его в уравнение (1.2), получим разностную схему вида
)
,
(
1 1
1 1
+
−
+
−
=
+
∑
β
∆
+
=
niniknniiutftuu , (1.10) где β
n
– постоянные коэффициенты, зависящие от числа использо- ванных для интерполяции точек
k, ∆
t – шаг
изменения независимой переменнной1
Разностные схемы вида (1.10) называются явными схемами
Адамса. Порядок аппроксимации явных схем Адамса равен числу использованных для интерполяции точек
k.
Аналогично можно получить неявные схемы Адамса. Для этого необходимо, чтобы при построении интерполяционного полинома
1
При построении многошаговых разностных схем с целью упрощения получаемых формул обычно рассматривается случай постоянного шага изменения независимой переменной.
14 рассматривалась также точка
ti+1
. В общем виде неявные схемы
Адамса отличаются от явных только начальной точкой при сумми- ровании:
)
,
(
1 1
0 1
+
−
+
−
=
+
∑
β
∆
+
=
niniknniiutftuu. (1.11)
Порядок аппроксимации неявных схем Адамса равен
k + 1.
При сравнении одношаговых и многошаговых разностных схем можно отметить, что:
1) с точки зрения условий устойчивости явные многошаговые разностные схемы не имеют преимуществ по сравнению с яв- ными одношаговыми схемами. Неявные многошаговые схемы абсолютно устойчивы только, если их порядок аппроксимации не выше второго;
2) одношаговые схемы требуют
больших вычислительных затрат на один шаг, т.к. промежуточные значения функции
f(t, u) не используются на следующем шаге, а пересчитываются заново.
В многошаговых разностных схемах при любом порядке ап- проксимации требуется пересчитывать только значение в самой удаленной точке
f(ti-k+1, u i-k+1), заменяя его на значение
f(ti+1, ui+1);
3) в то же время многошаговые схемы не являются самостартую- щими, т.к. требуют для начала расчета определения значений
u0
,
u1
, …
uk-1 1
. Это затрудняет изменение шага независимой пе- ременной в процессе расчета.
В качестве примера рассмотрим явные и неявные схемы Адамса
2-го и 4-го порядков.
Явная схема Адамса 2-го порядка. Построим полином первого порядка А(t) =
а0
+
а1
t, совпадающий с функцией
f(t, u) в точках
ti-1
и
ti. Коэффициенты полинома определим из решения системы уравнений
1
Стартовые значения обычно находят с помощью одношаговых разност- ных схем того же порядка аппроксимации.
15
)
,
(
;
)
,
(
1 0
1 1
1 1
1 0
i
i
i
i
i
i
i
i
f
u
t
f
t
a
a
f
u
t
f
t
a
a
≡
=
+
≡
=
+
−
−
−
−
(1.12)
Получим
t
f
f
a
t
a
f
a
i
i
i
i
∆
−
=
−
=
−
/
)
(
,
1 1
1 0
. Тогда интеграл I
i,i+1
будет равен
2
/
)
3
(
2
/
)
(
)
(
1 2
2 1
1 0
1 0
1
,
1
t
f
f
t
t
a
t
a
dt
t
a
a
I
i
i
i
i
t
t
i
i
i
i
∆
−
=
−
+
∆
=
+
=
−
+
+
∫
+
(1.13)
Явная схема Адамса 2-го порядка будет иметь вид
2
/
))
,
(
)
,
(
3
(
1 1
1
−
−
+
−
∆
+
=
i
i
i
i
i
i
u
t
f
u
t
f
t
u
u
. (1.14)
Неявная схема Адамса 2-го порядка.
Аналогично случаю явной схемы Адамса построим полином первого порядка А(t) = а
0
+ а
1
t, совпадающий с функцией f(t, u) в точках t
i
и t
i+1
. Коэффициенты полинома определим из решения системы уравнений
)
,
(
;
)
,
(
1 1
1 1
1 0
1 0
+
+
+
+
≡
=
+
≡
=
+
i
i
i
i
i
i
i
i
f
u
t
f
t
a
a
f
u
t
f
t
a
a
(1.15)
Получим
t
f
f
a
t
a
f
a
i
i
i
i
∆
−
=
−
=
+
/
)
(
,
1 1
1 0
. Тогда интеграл I
i,i+1
будет равен
2
/
)
(
2
/
)
(
)
(
1 2
2 1
1 0
1 0
1
,
1
t
f
f
t
t
a
t
a
dt
t
a
a
I
i
i
i
i
t
t
i
i
i
i
∆
+
=
−
+
∆
=
+
=
+
+
+
∫
+
(1.16)
Неявная схема Адамса 2-го порядка будет иметь вид
16 2
/
))
,
(
)
,
(
(
1 1
1
+
+
+
+
∆
+
=
i
i
i
i
i
i
u
t
f
u
t
f
t
u
u
. (1.17)
Схемы Адамса 4-го порядка аппроксимации.
Используя при оп- ределении интерполяционного полинома значения функции f(t, u) в четырех точках, можно аналогично получить явную:
24
/
))
,
(
9
)
,
(
37
)
,
(
59
)
,
(
55
(
3 3
2 2
1 1
1
−
−
−
−
−
−
+
−
+
−
∆
+
=
i
i
i
i
i
i
i
i
i
i
u
t
f
u
t
f
u
t
f
u
t
f
t
u
u
(1.18) и неявную схему Адамса 4-го порядка:
24
/
))
,
(
)
,
(
5
)
,
(
19
)
,
(
9
(
2 2
1 1
1 1
1
−
−
−
−
+
+
+
+
−
+
∆
+
=
i
i
i
i
i
i
i
i
i
i
u
t
f
u
t
f
u
t
f
u
t
f
t
u
u
. (1.19)
Особенности реализации неявных многошаговых разностных
схем.
Перепишем неявную схему (1.11) в виде:
)
,
(
)
,
(
1 1
1 0
1
i
n
i
i
i
i
i
i
u
u
g
u
t
f
t
u
u
K
+
−
+
+
+
+
β
∆
+
=
. (1.20)
В выражении (1.20) функция g
i
является известной. Для реше- ния полученного нелинейного относительно u
i+1
алгебраического уравнения используют два метода:
− метод простой итерации:
− метод Ньютона.
В методе простой итерации для искомого значения u
i+1
сначала определяется начальное приближение
0 1
+
i
u
, которое обычно полу- чают с помощью явной многошаговой схемы. Уточненные значе- ния u
i+1
находятся итерационно:
)
,
(
)
,
(
1 1
1 1
0 1
i
n
i
i
s
i
i
i
s
i
u
u
g
u
t
f
t
u
u
K
+
−
−
+
+
+
+
β
∆
+
=
. (1.21)
Здесь s – номер итерации.
Рассмотренная процедура называется методом прогноза- коррекции или методом предиктор-корректор. При реализации ме-
17 тода на каждом шаге можно выполнять либо фиксированное число итераций, либо проводить их до достижения заданной точности определения u
i+1
, т.е пока разница между значениями
1 1
−
+
s
i
u
и
s
i
u
1
+
не станет меньше заданной погрешности.
В ранних реализациях метода обычно использовали подход с выполнением итераций до совпадения соседних итерационных приближений искомой величины с машинной точностью. Однако даже в случае полной сходимости, полученное значение является только приближенным значением u
i+1
, а не точным значением функции Т(t
i+1
)
, поэтому такая точность определения u
i+1
является избыточной.
В настоящее время в методе прогноза-коррекции используют фиксированное количество итераций, обычно не более двух. Такой подход эффективен для большинства задач, но обладает одним крупным недостатком. При фиксированном числе итераций m ме- тод прогноза-коррекции является явным методом, т.к. конечное наилучшее приближение искомой функции
m
i
u
1
+
можно выразить по явной формуле через известные величины. Поэтому метод уже не обладает абсолютной устойчивостью.
В случае сильной нелинейности функции f(t, u) или при реше- нии так называемых жестких
1
систем уравнений метод прогноза- коррекции не работает и используется метод Ньютона.
Перепишем выражение (1.20) в виде уравнения F(u) = 0
)
,
(
)
,
(
)
(
0 1
1 0
i
n
i
i
i
i
u
u
g
u
t
f
t
u
u
u
F
K
+
−
+
−
β
∆
−
−
≡
=
, (1.22) тогда итерационный процесс по методу Ньютона организуется сле- дующим образом
1
Строгое определение жесткости требует достаточно детального описа- ния. Будем считать жесткой систему уравнений, моделирующую процес- сы с сильно отличающимися характерными временами. При решении та- ких систем с целью определения поведения медленно протекающих про- цессов приходится использовать шаги интегрирования, намного превы- шающие характерные времена для других рассматриваемых процессов.
18
( ) ( )
sisisisiuFuFuu1 1
1 1
1
+
+
+
+
+
′
−
=
. (1.23)
Здесь
uutftuFi∂
∂
β
∆
−
=
′
+
)
,
(
1
)
(
1 0
Начальное приближение
0 1
+
iu, как и методе прогноза-корреции, определяется с помощью явной мно- гошаговой схемы.
Сходимость метода Ньютона и количество необходимых итера- ций зависит от начального приближения в отличие от метода про- гноза-коррекции, для которого сходимость итераций определяется величиной шага изменения независимой переменной. Точность на- чального приближения при этом зависит от шага изменения
t.
Для реализации метода Ньютона требуется гораздо больше уси- лий, чем в случае метода прогноза-коррекции, но для жестких сис- тем уравнений ему практически нет альтернативы.
Скачать 1.6 Mb.