Санкт-Петербургский государственный университет
Кафедра моделирования электромеханических и
компьютерных систем
Рыжиков Александр Вячеславович
Выпускная квалификационная работа бакалавра
Исследование явных методов решения задачи
Коши основанных на разложении
Лагранжа – Бюрмана
Направление 010900
Прикладные математика и физика
Научный руководитель,
кандидат физ.-мат. наук,
доцент
Кривовичев Г. В.
Санкт-Петербург
2016
Содержание
Введение . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3
Глава 1. Обзор литературы . . . . . . . . . . . . . . . . . . . . . .
5
1.1
Метод рядов Тейлора . . . . . . . . . . . . . . . . . . . . . .
5
1.2
Явные методы Рунге – Кутты . . . . . . . . . . . . . . . . . .
6
1.3
ЯМРК на основе разложений Лагранжа – Бюрмана . . . . .
8
1.4
Устойчивость и жесткие задачи . . . . . . . . . . . . . . . .
10
1.5
Вывод . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
11
Глава 2. ЯМРК на основе разложений Лагранжа – Бюрмана
13
Глава 3. Исследование устойчивости . . . . . . . . . . . . . . .
17
Заключение . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
24
2
Введение
Обыкновенные дифференциальные уравнения (ОДУ) имеют важное
значение при моделировании многих прикладных задач в разных областях
знания. Однако получить точное решение задач можно лишь для некоторых классов ОДУ, для которых функция, входящая в правую часть, имеет простой вид. Для подавляющего числа задач, имеющих практическое
значение, получить решение в аналитическом виде невозможно. Поэтому
широкое распространение получили различные численные методы.
При построении и исследовании численных методов, мы сталкиваемся с проблемами устойчивости и точности метода. Проблема устойчивости
особенно усложняется в случае жестких систем, характерных для задач
химической кинетики, гидродинамики, теории электрических цепей. Известно, что для таких систем неявные методы дают лучший результат, чем
явные методы. Но неявные методы имеют дополнительные сложности в
реализации. В связи с этим актуально построение и исследование явных
методов типа Рунге – Кутты с расширенной областью устойчивости.
В [1] предлагается подход к построению явных методов Рунге – Кутты
на основе разложений Лагранжа – Бюрмана. За счет выбора функции, по
которой идет разложение в формуле Лагранжа – Бюрмана, ее параметров
и параметров метода Рунге – Кутты можно влиять на свойства таких методов.
Работа состоит из введения, трех глав и заключения.
Первая глава является обзорной. В ней вводятся основные понятия
и определения: задача Коши, метод рядов Тейлора, явные методы Рунге – Кутты (ЯМРК), формула Лагранжа – Бюрмана и ЯМРК на ее основе, функция устойчивости, область устойчивости, жесткие задачи, схемы
неявных методов, с которыми сравниваются исследуемые методы. В конце
cформулированы цель и задачи работы.
Вторая глава посвящена явным методам Рунге – Кутты на основе разложений Лагранжа – Бюрмана.
3
Третья глава посвящена исследованию устойчивости. Получены графики областей устойчивости. Приведено сравнение значений площадей
областей устойчивости. Решены тестовые задачи. Приведено сравнение с
неявными методами.
В заключении приведены основные результаты работы.
4
Глава 1. Обзор литературы
Будем рассматривать задачу Коши (1)–(2) для скалярного уравнения, но все результаты можно обобщить на случай систем.
du
= f (x, u),
dx
x ∈ [a, b],
u(a) = u0 ,
(1)
(2)
где x – независимая переменная, f (x, u) – скалярная функция. Полагаем,
что функция f (x, u) определена и непрерывно дифференцируема необходимое число раз на R2 .
Далее речь пойдет об одношаговых численных методах. Это методы,
которые последовательно дают приближения uj ≈ u(xj ) на основе известного приближения uj−1 ≈ u(xj−1 ) на сетке:
Ω = {xj , j = 0, N |x0 = a, xj = xj−1 + hj , j = 1, N − 1, xN = b}.
Будем считать, что hj = h = const, т.е. сетка является равномерной.
1.1. Метод рядов Тейлора
Одним из подходов при построении расчетных схем различных порядков точности для решения уравнения (1) является использование разложения решения рассматриваемой задачи Коши по формуле Тейлора [1]:
h2 00
hp (p)
u(xj+1 ) ≈ u(xj ) + hu (xj ) + u (xj ) + · · · + u (xj ).
2
p!
0
(3)
Производные, входящие в правую часть (3) могут быть фактически
найдены:
u0 (xj ) = f (xj , uj ),
u00 (xj ) = {fx0 + f fu0 }|xj ,
00
00
00
u000 (xj ) = {fxx
+ 2f fxu
+ f 2 fuu
+ (fx0 + f fu0 )fu0 }|xj ,
.........
5
(4)
Например, полагая в (3) p = 1 и p = 2, заменяя u(xj ) на uj с учетом
(4) получим расчетные схемы [3]:
uj+1 = uj + hf (xj , uj ) (явный метод Эйлера),
uj+1 = uj + h[f (xj , uj ) + h2 (fx0 + f fu0 )].
Как видно, такие схемы не требуют вычисления дополнительных
условий и позволяют легко менять шаг интегрирования [4]. Однако эти
методы имеют свои недостатки. Так из теории аппроксимации функций
многочленами известно, что такие приближения приводят к колебаниям
интерполирующей функции в окрестности разрывов исходной функции и
в областях больших градиентов решения [1], а практическое применение
ограничено задачами, для которых легко вычисляются полные производные высшего порядка.
1.2. Явные методы Рунге – Кутты
Определим общую схему построения явных методов Рунге – Кутты.
Пусть значение u(xj ) известно. Тогда значение решения в узле xj+1
вычисляется по формуле:
u(xj+1 ) ≈ z1 = u(xj ) + ∆uh,j .
(5)
Формула для вычисления ∆uh,j зависит от количества этапов используемого метода Рунге – Кутты. В общем случае m -этапного метода сначала
последовательно вычисляют величины, которые иногда называют корректирующими функциями:
k1 = hf (xj , u(xj )),
kq = hf (xj + cq h, u(xj ) +
q−1
P
aqν kν (h)),
q = 2, ..., m,
(6)
ν=1
где a21 , a31 , ..., am1 , am2 , ..., am,m−1 , b1 , ..., bm , c1 , ..., cm ∈ R – параметры метода.
6
Затем приращение решения ∆uh,j находят в виде линейной комбинации:
∆uh,j =
m
X
bi ki .
(7)
i=1
ЯМРК (5)–(7) имеет порядок точности p, если разложения Тейлора для точного решения u(xj+1 ) и полученного приближенного к нему z1
совпадают до члена hp включительно.
Введем в рассмотрение функцию:
Ψ(h) = u(xj + h) − u(xj ) −
m
X
bi ki
i=1
– методическая погрешность одношагового метода.
Для того, чтобы получить метод конкретного порядка точности, нужно выбрать его параметры так, чтобы выполнялись следующие равенства:
Ψ0 (0) = Ψ00 (0) = · · · = Ψ(p) (0) = 0,
(8)
причем Ψ(p+1) (0) 6= 0.
Фиксируя число этапов m ЯМРК (5)–(7) мы фактически фиксируем
количество параметров, подлежащих определению, а задавая порядок p метода, с учетом выполнения (8) формируем систему нелинейных алгебраических уравнений, связывающую параметры m-этапного ЯМРК p-порядка
– условия порядка.
Любое частное решение системы условий порядка дает расчетную
схему. При решении этой системы, когда m > 1, мы сталкиваемся с тем,
что неизвестных больше, чем уравнений. Поэтому один или несколько параметров полагаются в качестве входных. Обычно входные параметры подбираются так, чтобы получить удобные для вычислений расчетные формулы,
и такие, которые обеспечивают подавление того или иного члена в главном
члене погрешности [2].
Например, классический метод Рунге – Кутты 4-го порядка [3] имеет
вид:
h
u(xj+1 ) = u(xj ) + (k1 + 2k2 + 2k3 + k4 ),
6
7
hk1
h
,
k1 = f (xj , uj ), k2 = f xj + , uj +
2
2
h
hk2
k3 = f xj + , uj +
, k4 = f (xj + h, uj + hk3 ).
2
2
ЯМРК так же как и методы Тейлора не требуют вычисления до
полнительных условий и позволяют легко менять шаг интегрирования. К
тому же здесь не вычисляются производные от функции f (x, u), поэтому
область применения таких методов больше [4]. Однако они не пригодны
для решения жестких задач, которые часто встречаются в приложениях.
1.3. ЯМРК на основе разложений
Лагранжа – Бюрмана
Как отмечается в [1], функцию можно разложить и в степенной ряд
общего вида, т.е. по степеням некоторой другой функции ϕ(x). Тогда за
счет выбора функции ϕ(x) можно уменьшить амплитуду колебаний численного решения. Разложить функцию в ряд по степеням другой функции
можно с помощью формулы Лагранжа – Бюрмана [5].
Пусть u(x) и ϕ(x) – достаточно гладкие функции и ϕ(0) = 0, ϕ0 (0) 6=
0, т.е. точка x = 0 является простым нулем функции ϕ(x). Тогда формулу
Лагранжа – Бюрмана можно записать в виде:
u(x) ≈ u(x0 ) +
m
X
ϕk (x − x0 )bk ,
(9)
k=1
1 n dk−1 h 0 x − x0 k io
u (x)
, k = 1, 2, . . . .
(10)
bk = lim
x→x0 k! dxk−1
ϕ(x − x0 )
С помощью формулы (9)–(10) можно разложить функцию u(x) в ряд
по степеням функции ϕ(x − x0 ). Как видно при ϕ(x) = x формула Лагранжа – Бюрмана переходит в формулу Тейлора. Стоит отметить, что сложность вычисления пределов в (10) нелинейно возрастает с ростом k [1]. В
частности, для b1 , b2 имеем:
u0 (x0 )
b1 = 0
,
ϕ (0)
8
(11)
u0 (x0 )ϕ00 (0)
u00 (x0 )
.
b2 = 0 2 −
2ϕ (0)
2ϕ0 (0)3
(12)
Далее предлагается строить ЯМРК на основе разложений Лагранжа – Бюрмана. Применять разложение Лагранжа – Бюрмана для построения разностных схем впервые предложил Е. В. Ворожцов [1].
Общая схема ЯМРК на основе разложений (9)–(10) имеет вид:
u(xj+1 ) ≈ z1 = u(xj ) + ∆uh,j ,
∆uh,j =
m
X
(13)
Ai ki ,
(14)
i=1
k1 = ϕ(h)f (xj , u(xj )),
kq = ϕ(h)f (xj + cq h, u(xj ) +
q−1
P
aqν kν (h)),
q = 2, ..., m.
(15)
ν=1
где a21 , a31 , ..., am1 , am2 , ..., am,m−1 , b1 , ..., bm , c1 , ..., cm ∈ R – параметры метода.
ЯМРК (13)–(15) имеет порядок точности p, если разложения Лагранжа – Бюрмана для точного решения u(xj+1 ) и полученного приближенного
к нему z1 совпадают до члена ϕp (h) включительно.
Одна из проблем возникающая при построении методов Рунге –
Кутты на основе разложений Лагранжа – Бюрмана связана с выбором
функции ϕ(x). В работе были рассмотрены функции:
x − x
j
β h,
ϕ(x − xj ) = th
h
x − x
j
ϕ(x − xj ) = arctg
β h,
h
где β – положительная постоянная.
(16)
(17)
Функция (16) была предложена в работе [1], а функция (17) предложена в этой работе. Эти функции рассматриваются с тем, что за счет
выбора параметра β можно сделать так, чтобы функции сильно менялись
на маленьких промежутках и почти не изменялись в остальной области
определения, что важно для решения жестких задач.
9
На свойства метода можно влиять за счет выбора:
1. функции ϕ(x);
2. параметров метода Рунге – Кутты.
1.4. Устойчивость и жесткие задачи
Функцию устойчивости можно интерпретировать как численное решение (19) для одного шага тестовой задачи Далквиста (18) [6]:
u0 (x) = λu(x),
u(x0 ) = 1,
u(xj ) = R(z)u(xj−1 ).
λ ∈ C.
(18)
(19)
Множество S = {Z ∈ C; |R(z)| 6 1} – называют областью устойчивости [6].
Дать определение «Что такое жесткие задачи?» строго математически не представляется возможным. Однако для систем с постоянной матрицей n × n, можно определять «жесткость» задачи по собственным значениям матрицы, так, если
1. Reλk < 0 для любого k = 1, 2, . . . , n,
2. Число жесткости g =
max |Re(λk )|
min |Re(λk )|
1,
то говорят, что задача жесткая [3]. Можно считать, что число жесткости
велико, если оно больше или равно 10 [4].
Для жестких систем характерно наличие двух участков решения: пограничный слой (область, где компоненты решения характеризуется быстрыми изменениями) и квазистационарный режим (медленные изменения).
Однако пограничный слой может отсутствовать при специальном выборе
начальных условий.
Рассмотрим систему [7]:
du
= Ju,
dt
10
(20)
где t ∈ [0, 0.5] – время, J – матрица второго порядка:
!
−1000 999
J=
.
1
−2
С двумя видами начальных условий:
u1 (0) = u2 (0) = 1,
u1 (0) = −1,
u2 (0) = 1.
(21)
(22)
Легко убедится, что собственные числа матрицы J удовлетворяют
двум поставленным условиям. Данная задача является тестовой и имеет
аналитическое решение:
u1 (t) = 0.999 u1 (0) − u2 (0) e−1001t + 0.001u1 (0) + 0.999u2 (0) e−t ,
u2 (t) = −0.001 u1 (0) − u2 (0) e−1001t + 0.001u1 (0) + 0.999u2 (0) e−t .
Видно, что для задачи (20), (22) внутри пограничного слоя компонента решения u1 (t) меняется заметно быстрее, чем u2 (t). Поэтому иногда
u1 (t) называют «быстрой» компонентой, а u2 (t) – «медленной» [7]. В случае
выбора начальных условий (21) – пограничный слой отсутствует, и задача
(20), (21) является «умеренно жесткой».
Ниже представлены схемы неявного метода Эйлера и неявного метода трапеций [3], с которыми сравнивались полученные явные методы
Рунге – Кутты на основе разложений Лагранжа – Бюрмана:
uj+1 = uj + hf (xj+1 , uj+1 ),
h
uj+1 = uj + (f (xj , uj ) + f (xj+1 , uj+1 )).
2
1.5. Вывод
В литературе представлен новый подход к построению явных методов Рунге – Кутты на основе разложений Лагранжа – Бюрмана. Такие методы обладают рядом преимуществ по сравнению с классическими методами Рунге – Кутты. Например, за счет выбора функции, по которой идет
11
разложение в формуле Лагранжа – Бюрмана, ее параметров и параметров
метода можно влиять на свойства метода в целом. Данные методы были
предложены для решения жестких задач Коши. К сожалению, к настоящему моменту нельзя утверждать, что свойства таких методов являются
в полной мере исследованными.
Таким образом, цель работы заключается в исследовании явных
методов Рунге – Кутты, построенных на основе разложений Лагранжа –
Бюрмана. Для достижения цели необходимо решить следующие задачи:
1. Построение методов с малым числом этапов (до 3).
2. Исследование устойчивости в зависимости от значений параметра.
3. Решение тестовых задач.
12
Глава 2. ЯМРК на основе разложений
Лагранжа – Бюрмана
Методы 1–2 этапов представлены в [1], где описан алгоритм их построения. Например, одноэтапный метод можно получить, если положить
m = 1 в формуле Лагранжа – Бюрмана, тогда из (15) с учетом (11) и (1)
получим:
u(xj+1 ) = u(xj ) +
ϕ(h)
f
x
,
u(x
)
j
j .
ϕ0 (0)
(23)
Для краткости метод (23) будем называть методом LB1.
Полагая m = 2 в (14) из (13), (15) получим схему двухэтапного метода:
u(xj+1 ) = u(xj ) + b1 k1 + b2 k2 ,
(24)
k1 = ϕ(h)f (xj , u(xj )),
(25)
k2 = ϕ(h)f (xj + c1 h, u(xj ) + a21 k1 ).
Условия порядка выведены в [1], так если функция ϕ(x) – нечетная,
условия порядка имеют вид:
A1 + A2 =
1
,
ϕ0 (0)
c1 A2 h =
ϕ(h)
,
2ϕ0 (0)2
A2 a21 =
1
.
2ϕ0 (0)2
(26)
Как можно видеть система условий порядка – это система из 3 уравнений с 4 неизвестными, поэтому один параметр полагается произвольным.
Очевидно, что в качесте входного параметра надо брать A2 , зная его мы
легко найдем остальные. Следуя [1], из принципа минимизации главного
члена погрешности, выберем A2 =
3
4ϕ0 (0) .
Тогда разрешая систему (26) по-
лучим:
A2 =
3
,
4ϕ0 (0)
A1 =
1
,
4ϕ0 (0)
c1 =
2ϕ(h)
,
3hϕ0 (0)
a21 =
2
.
3ϕ0 (0)
(27)
Для краткости метод (24), (25), (27) будем называть методом LB2.
Подробнее остановимся на построении трехэтапного метода для автономной системы. Положим в (14) m = 3, тогда приращение ∆uh,j в (13)
13
с учетом (15) вычисляется по формулам:
∆uh,j = A1 k1 + A2 k2 + A3 k3 ,
(28)
k1 = ϕ(h)f (u), k2 = ϕ(h)f (u + a21 k1 ), k3 = ϕ(h)f (u + a31 k1 + a32 k2 ). (29)
Выпишем выражение для приращения по формуле Лагранжа –
Бюрмана:
∆uh = b1 ϕ(h) + b2 ϕ2 (h) + b3 ϕ3 (h) + o(ϕ3 (h)).
Коэффициенты имеют вид:
(30)
u0 (x)
b1 = 0 ,
ϕ (0)
b2 =
u0 (x)ϕ00 (0)
u00 (x)
−
,
2ϕ0 (0)2
2ϕ0 (0)3
u0 (x)ϕ00 (0)2 u00 (x)ϕ00 (0)
u000 (x)
u0 (x)ϕ000 (0)
b3 =
−
+ 0 3−
.
2ϕ0 (0)5
2ϕ0 (0)4
6ϕ (0)
6ϕ0 (0)4
Сравним величину (28) с величиной (30). Для этого разложим k2 и
k3 из (29) по формуле Лагража – Бюрмана и найдем разность δu = ∆uh −
∆uh,j .
Подставим k1 в k2 и разложим f по формуле Лагранжа – Бюрмана:
k2 = ϕ(h)(f + ϕ(a21 ϕ(h)f (u))C1 + ϕ2 (a21 ϕ(h)f (u))C2 ) + o(ϕ2 (h)),
здесь C1 =
fu
ϕ0 (0) ,
C2 =
fuu
2ϕ0 (0)2 .
Затем ϕ(a21 ϕ(h)f (u)) разложим по формуле Тейлора по степеням
ϕ(h) в окрестности нуля, т.к. h 1 и учтем, что ϕ(0) = 0:
ϕ00 (0) 2
ϕ(a21 ϕ(h)f (u)) = ϕ (0)ϕ(h)a21 f +
ϕ (h)a221 f 2 + o(ϕ2 (h)).
2
0
Окончательно для k2 имеем следующее разложение:
ϕ00 (0) 2
k2 = ϕ(h)(f + (ϕ (0)ϕ(h)a21 f +
ϕ (h)a221 f 2 )C1 +
2
0
ϕ0 (0)2 ϕ2 (h)a221 f 2 C2 ) + o(ϕ3 (h)).
Аналогично делаем для k3 .
14
В результате получено следующее выражение:
δu = r1 ϕ(h) + r2 ϕ2 (h) + r3 ϕ3 (h) + o(ϕ3 (h)).
(31)
Здесь
1
r1 =
− A1 − A2 − A3 f,
ϕ0 (0)
r2 = r2,0 f + r2,1 f fu ,
r3 = r3,0 f + r3,1 f 2 fuu + r3,2 f fu + r3,3 f fu2 ,
r2,0
r3,0 =
ϕ00 (0)
= − 0 3,
2ϕ (0)
ϕ00 (0)2 ϕ(3) (0)
−
2ϕ0 (0)5 6ϕ0 (0)4
r2,1 =
r3,1 =
1
− A2 a21 − A3 a31 − A3 a32 ,
2ϕ0 (0)2
1
1
1
1
2
2
−
A
a
−
A
a
−
A3 a232 −A3 a31 a32 ,
2
3
21
31
0
3
6ϕ (0) 2
2
2
ϕ00 (0)
1
r3,2 = − 0 4 , r3,3 = 0 3 − A3 a21 a32 .
2ϕ (0)
6ϕ (0)
Для того чтобы получить метод третьего порядка точности надо занулить функции r1 , r2 , r3 из (31). Как можно видеть, занулить член r3,0
нельзя (третья производная не обнуляется). Поэтому построенный трехэтапный метод будет иметь второй порядок точности. Условия порядка
имеют вид:
A1 + A2 + A3 =
1
ϕ0 (0)
A2 a21 + A3 a31 + A3 a32 =
,
(32)
1
2ϕ0 (0)2
,
(33)
1
1
1
1
A2 a221 + A3 a231 + A3 a232 + A3 a31 a32 = 0 3 ,
(34)
2
2
2
6ϕ (0)
1
A3 a21 a32 = 0 3 .
(35)
6ϕ (0)
Трехэтапный метод для краткости будем называть методом LB3.
Как можно видеть условия порядка для метода LB3 представляют
собой систему из 4 уравнений с 6 неизвестными, поэтому в качестве входных надо выбрать два параметра. Из (35) следует, что A3 , a21 , a32 не равны
нулю. Выберем входными параметрами a21 и a32 , тогда из (35) сразу найдем A3 , а уравнение (34) будет квадратным только относительно a31 . Затем
15
подставим A3 в (32), (33), (34). Далее из (33) выразим A2 через a31 и подставим в (34). Разрешая (34) относительно a31 получим два корня. Затем
последовательно найдем A2 и A1 . Таким образом, система условий порядка
метода LB3 имеет два решения:
1
,
6a21 a32 ϕ0 (0)3
q
1
a21 − 2a32 ± a221 + 8a21 a32 − 12a221 a32 ϕ0 (0) ,
a31 =
2
a31
1
1
1
−
−
,
A
=
− A2 − A3 .
A2 =
1
2a21 ϕ0 (0)2 6a221 a32 ϕ0 (0)3 6a221 ϕ0 (0)3
ϕ0 (0)
A3 =
16
Глава 3. Исследование устойчивости
Для построенных методов 1–3 этапов проведено исследование на
устойчивость. Задача заключалась в том, чтобы посмотреть как выбор функции ϕ(x), по которой идет разложение в формуле Лагранжа –
Бюрмана, и параметр β влияют на область устойчивости.
Применяя методы 1–3 этапов к решению задачи Далквиста (18) на
одном шаге, получаем соответственно следующие виды функций устойчивости:
R(z) = 1 + γz,
1
R(z) = 1 + γz + γ 2 z 2 ,
2
1
1
R(z) = 1 + γz + γ 2 z 2 + γ 3 z 3 ,
2
6
где z = λh, γ =
ϕj (h)
hϕ0j (0) .
Для функции (17) γ =
arctg(β)
,
β
а для функции (16) γ =
th(β)
β .
Построены графики областей устойчивости на комплексной плоскости (рис. 1). По графикам можно сделать следующие выводы:
1. Области устойчивости увеличиваются при увеличении значения параметра β.
2. Область устойчивости метода LB1 с функцией th(β) больше, чем с
функцией arctg(β) при одном и том же значении параметра. Аналогичная ситуация для методов LB2 и LB3.
3. Области устойчивости увеличиваются при увеличении этапов метода.
17
-35
-30
-25
-20
-15
-10
30
30
20
20
10
10
-50
-5
-40
-30
-20
-10
-10
-10
-20
-20
-30
-30
Β
5
10
15
20
25
б)
а)
30
40
Β
20
20
10
-50
-35
-30
-25
-20
-15
-10
-40
-30
-20
-10
-5
-10
-20
-20
-40
5
10
15
20
25
-30
г)
в)
60
40
40
20
20
-60
-40
-30
-20
-50
-40
-30
-20
-10
-10
-20
-20
Β
5
10
15
20
25
-40
-60
-40
е)
д)
Рис. 2. Изменение области устойчивости в зависимости от значения параметра.
a) LB1 c функцией arctg(β), б) LB1 c функцией th(β), в) LB2 c функцией arctg(β),
г) LB2 c функцией th(β), д) LB3 c функцией arctg(β), е) LB3 c функцией th(β).
18
В табл. 1 представлены результаты вычислений площадей областей
устойчивости для методов 1–3 этапов при различном выборе функции ϕ и
в зависимости от значения параметра β.
Табл. 1. Значение площади устойчивости в зависимости от значения параметра β.
arctg(β)
th(β)
β
LB1
LB2
LB3
LB1
LB2
LB3
2
10.2375
19.1662
29.7215
15.5300
25.2670
39.1917
3
18.1246
33.8339
52.5623
28.6074
53.4235
82.8488
4
28.6029
53.5732
82.9452
50.2842
94.1213
145.9296
5
41.7012
77.9514
120.7036
78.5380
146.5652 227.7175
6
57.2537
106.8258 165.6500 113.1759 211.5690 327.6717
7
75.2708
140.9964 218.6308 153.9961 287.6746 446.3889
8
96.1558
179.9344 278.7758 201.2440 376.0992 583.0342
9
119.2364 223.1580 346.3325 254.4439 475.4141 737.9864
10 145.2145 270.9186 420.7657 314.2157 587.2724 909.8922
Используя полученные данные, построены графики (рис. 2). По графикам видно, что площади областей устойчивости методов LB1, LB2 и
LB3 неограниченно возрастают при увеличении значения параметра β. При
этом оказалось, что для случая, когда в качестве функции ϕ берем гиперболический тангенс, площади областей устойчивости больше (рис. 1, 2).
SH ΒL
SH ΒL
140
120
100
80
60
40
20
250
400
200
300
150
200
100
th(Β)
arctg(Β)
50
Β
4
6
LB1
8
10
100
Β
4
6
LB2
8
10
4
6
8
10
β
LB3
Рис. 2. Сравнение значений площадей областей устойчивости.
Были рассмотрены две тестовые задачи (20), (21) и (20), (22) с той целью, чтобы на практике посмотреть как выбор параметра и сетки влияют
на устойчивость. Задача заключалась в определении значения параметра
19
β, после которого метод LB1 устойчив при данном разбиении сетки. Параметр подбирался вручную, основываясь на графическом представлении
решения, как показано на рис. 3.
1.10
1.00
1.05
0.98
0.96
0.95
u1, u2
u1, u2
1.00
0.90
0.94
0.85
0.92
0.80
0.90
0.75
0.0
0.1
0.2
0.3
0.4
0.5
0.0
0.1
0.2
0.3
0.4
0.5
t
t
б)устойчивое
а) неустойчивое
Рис. 3. Графики решения задачи (20), (22) по схеме LB1 c функцией arctg(β)
Расчеты проводились на равномерных сетках различных разбиений.
Результаты занесены в табл. 2, где n – число узлов сетки.
Табл. 2. Значение параметра β, после которого метод LB1 будет устойчив.
к задаче (20), (21) к задаче (20), (22)
arctg(β)
n
th(β)
arctg(β)
th(β)
устойчив при β >
100
2.6
2.1
8.8
6
125
2
1.7
6.8
4.8
150
1.5
1.3
5.5
4
175
1.1
1
4.6
3.4
200
0.8
0.7
3.9
3
225
0.4
0.4
3.4
2.6
Как видно из табл. 2 устойчивость улучшается при увеличении значения параметра, что подтверждает теоретические результаты.
Таким образом, за счет увеличения значения параметра можно существенно улучшить свойства устойчивости, задавая область устойчивости
20
сколь угодно большой площади. Однако при больших значениях параметра существенно ухудшается точность численного решения. Эту ситуацию
можно проиллюстрировать при решении простейшей тестовой задачи:
dy
= λy,
dt
y(0) = 1,
t ∈ [0, 10],
λ ∈ R,
для которой известно точное решение y(t) = eλt .
Численное решение будем искать в виде y(tj ) = q j y0 .
Для метода LB1:
q =1+
λh
arctg(β).
β
Для метода LB2:
λh
λ2 h2
q =1+
arctg(β) +
arctg2 (β).
β
2β
а)
б)
Рис.4. Графики решения задачи (36). Численное решение получено по методу LB1.
a)λ = −0.1, β = 3, б)λ = −0.1, β = 0.1.
Как можно видеть (см. рис. 4 а)), при значении β на порядок больше |λ|, погрешность численного решения весьма существенна. Это связано
с тем, что решение полученное по разностной схеме стремится к нулю с
меньшей скоростью, чем точное решение. Если брать параметр величиной
одного порядка с |λ|, то ситуацию можно исправить, при этом оказалось,
что шаг сетки ни на что не влияет (см. рис. 4 б)).
При решении жесткой тестовой задачи (20), (22) было проведено
сравнение методов 1-го и 2-го порядков точности с неявными методами
21
таких же порядков (неявным методом Эйлера и неявным методом трапеций), на одних и тех же сетках, построенных с постоянным шагом.
Погрешности численного решения вычислялись по формуле:
∆uk = kuk −
u∗k k2
N
−1
i0.5
X
1
∗ 2
=
(ukj − ukj ) (tj+1 − tj ) ,
tN − t0 j=0
h
k = 1, 2,
где u∗k – точное решение.
Результаты представлены в табл. 3, 4
Табл. 3. Погрешности численных решений для методов первого порядка точности.
Сетка 500 узлов
Сетка 1000 узлов
Метод LB1 Метод BEM
Метод LB1 Метод BEM
∆u1
3.54 × 10−2
1.77 × 10−2
1.34 × 10−2
9.98 × 10−3
∆u2
−5
−4
−5
−5
9.98 × 10
1.01 × 10
4.53 × 10
quick
5.09 × 10
slow
Здесь
Метод LB1 – метод первого порядка точности на основе разложений
Лагрнажа – Бюрмана с функцией (17) и β = 0.01;
Метод BEM (backward Euler method) – неявный метод Эйлера.
quick – быстрая компонента решения, slow – медленная.
Табл. 4. Погрешности численных решений для методов второго порядка точности.
Сетка 500 узлов
Сетка 1000 узлов
Метод LB2
Метод TR
Метод LB2
Метод TR
∆u1
3.54 × 10−2
4 × 10−3
1.34 × 10−2
1.34 × 10−4
∆u2
−4
−5
−5
1.96 × 10
−5
9.99 × 10
4.53 × 10
9.43 × 10
quick
slow
Метод LB2 – метод второго порядка точности на основе разложений
Лагрнажа – Бюрмана с функцией (17) и β = 0.01;
Метод TR (trapezoidal rule) – неявный метод трапеций.
Как можно видеть (см. табл. 3, 4), для быстрой компоненты решения
неявные методы позволяют получить более точные результаты, в то время
22
как для медленной компоненты погрешности примерно равны. Несмотря
на то, что у явных методов точность оказалась хуже, погрешность все равно является достаточно малой, при этом явные методы являются менее
затратными в реализации.
Таким образом, в главе 3 получены следующие результаты:
1. Проведен анализ устойчивости методов LB1, LB2, LB3.
2. Решены тестовые задачи. Приведены сравнения с результатами, полученными неявными методами.
Из полученных результатов можно сделать вывод: при увеличении
числа этапов методов Рунге – Кутты на основе разложений Лагранжа –
Бюрмана улучшается их устойчивость, но при этом ухудшается точность
полученных результатов.
23
Заключение
Основные результаты работы:
1. Построены явные методы Рунге – Кутты на основе разложений
Лагранжа – Бюрмана 1–3 этапов. Трехэтапный метод построен впервые.
2. Проведено исследование устойчивости явных методов Рунге – Кутты
1–3 этапов, построенных на основе разложений Лагранжа – Бюрмана.
Показано, что область устойчивости увеличивается при увеличении
значения параметра β.
3. При решении тестовых задач проведено сравнение с неявными методами. Показано, что явные методы Рунге – Кутты на основе разложений Лагранжа – Бюрмана могут применяться для решения жестких
задач.
24
Список литературы
1. Ворожцов Е. В. Построение явных разностных схем для обыкновенных дифференциальных уравнений с помощью разложений ЛагранжаБюрмана // Вычислительные методы и программирование, 2010. Т. 11,
с. 198–209
2. Арушанян О. Б., Залеткин C. Ф. Численное решение обыкновенных
дифференциальных уравнений на Фортране. М.: Изд-во МГУ, 1990.
336 с.
3. Вержбицкий В. М. Основы численных методов: Учебник для вузов. М.:
Высшая школа, 2002. 840 с.
4. Холл Дж., Уатт Дж. Современные численные методы решения обыкновенных дифференциальных уравнений. Пер. с англ. под ред. Горбунова
А. Д. / М.: Мир, 1979. 312 c.
5. Лаврентьев М. А., Шабат Б. В. Методы теории функций комплексного
переменного. М: Наука, 1973. 749 с.
6. Хайрер Э., Ваннер Г. Решение обыкновенных дифференциальных уравнений. Жесткие и дифференциально-алгебраические задачи. Пер. с англ. под ред. Филиппова С. C. / М.: Мир, 1999. 685 с.
7. Ракитский Ю. В., Устинов С. М., Черноруцкий И. Г. Численные методы
решения жестких систем. М.: Наука, 1979. 208 c.
25
Отзывы:
Авторизуйтесь, чтобы оставить отзыв