Разностные схемы для решения жестких обыкновенных дифференциальных уравнений в пространстве неопределенных коэффициентов (стр. 3 )

		0
			(39)
			(b21 ¹ 0).

Схема Хаммера–Холлингсворта — 2-го порядка, неявная, не L-устойчивая (q(–¥) = 1; пример неудачной схемы):

Явная схема Эйлера–Коши (точка O ):

0	0	0
1	1	0	(40)
	0	½

Две полуявные схемы 3-го порядка, не L-устойчивые, немонотонные, A-устойчивые:

1) точка B1 (схема Розенброка)

		0
			(41)

2) точка B2, для которой

1	1	0
1/3	–1/3	2/3	(42)
	1/4	3/4

Основная гипотеза, принимаемая здесь и в дальнейшем при анализе разностных схем: точкам в пространстве неопределенных коэффициентов {bkl}, которые близки между собой (в смысле расстояния, например, в евклидовой метрике), соответствуют схемы, близкие по своим свойствам (точности, устойчивости и т. п.).

1.4. Линейные многошаговые схемы
(методы типа Адамса)

1.4.1. Алгоритм и аппроксимация

Используя, как и в разделе 1.3, метод неопределенных коэффициентов и записывая линейную комбинацию вектор-функций v и f в некоторой последовательности равноотстоящих точек tn, tn+1, …, tn+k (t  = tn+kt – tn+k = const), для линейных многошаговых методов получим следующее выражение:

. (43)

Из-за однородности (43), коэффициент при искомом значении vn+k можно выбрать единичным:

аK = 1. (44)

Неопределенные коэффициенты ak, bk ( k = 0, …, K ) из разложения (43) в ряд Тейлора в точке tn+K (или tn) связаны условиями аппроксимации (на решениях (18)):

,  (45)

(обеспечивают 1-й порядок аппроксимации),

(46)

(обеспечивает 2-й порядок аппроксимации вместе с (45)),

…

(47)

(обеспечивает n-й порядок аппроксимации на решениях (18) вместе с предшествующими условиями).

Схемы первого порядка аппроксимации малоинтересны, поэтому в дальнейшем будем рассматривать схемы второго или более высокого порядка точности. Если с учетом нормировки (44) и условий аппроксимации второго порядка точности (45), (46) исключить, например,

,

,

,

то оставшиеся свободными коэффициенты (если позволяет выбранное K > 2), можно принять за линейное пространство размерности 2(K–1), например, с евклидовой метрикой. Каждой точке в этом пространстве будет соответствовать некоторая разностная схема второго порядка точности, а условия более высокого порядка аппроксимации ((47) при n = 3 и т. д.) после исключения в них aK, aK–1, bK, bK–1, являясь относительно оставшихся свободными коэффициентов ak, bk ( k = 0 , …, K–2 ) линейными уравнениями

,

образуют в этом пространстве соответствующую гиперплоскость схем третьего порядка аппроксимации. На пересечении двух таких гиперплоскостей ((47) при n = 3 и n = 4) будут расположены схемы с четвертым порядком аппроксимации и т. д., пока мы не исчерпаем все возможности для выбранного K, т. е. пока не найдем единственную схему с наиболее высоким порядком аппроксимации (точку в этом пространстве).

В частности, при K = 2 имеем:

a2 = 1, a1 = –1 – a0,

b2 = (1 + a0 + 2b0) / 2, b1 = (1 – 3a0 – 4b0) / 2, (48)

где a0 и b0 произвольны для схем 2-го порядка точности.

Условие 3-го порядка аппроксимации дает прямую в плоскости {a0, b0} с уравнением

5a0 + 12b0 + 1 = 0,

а единственной схемой 4-го порядка аппроксимации будет точка на этой прямой с координатами:

a0 = –1, b0 = 1/3. (49)

1.4.2. Устойчивость

Как уже отмечалось, одной аппроксимации недостаточно для сходимости решений (43) к решениям исходной дифференциальной задачи даже в линейном случае (7). Нужно еще обеспечить устойчивость разностных схем (43). Для одного линейного уравнения, полагая f = lv, из (43) получим

. (50)

Устойчивость разностных схем (50) будем исследовать на специальных решениях вида vn = qn. Тогда для определения q имеем многочлен устойчивости степени K и уравнение

,

корни которого qj, j = 1, 2, …, K (действительные или комплексные, в зависимости от s, ak, bk) должны удовлетворять условиям устойчивости

|q(s)| ≤ 1, j = 1, …, K. (51)

Если в плоскости {Re s, Im s} область устойчивости, определяемая неравенствами (51), включает всю левую полуплоскость (т. е. (51) выполняется для всех s с Re s < 0), имеем A-устойчивую схему. Если эта область включает в себя заштрихованную на рис. 1.11а или 1.11б часть плоскости {Re s, Im s }, имеем жестко-устойчивую схему. Для действительных s и qj(s) условия устойчивости схемы будут иметь вид –1 ≤ qj(s) ≤ 1. Если это условие выполняется для всех s: –¥ < s ≤ 0, то схема A-устойчива. Если

|q(s)| = 0, j = 1, 2, …, K, (52)

то схема будет L-устойчивой.

1.4.3. Примеры линейных многошаговых схем

Для случая K = 2 с учетом (48), (50) имеем

vn+2 = a1vn+1 + a0vn, q2 – a1q – a0 = 0.

Отсюда

, ,

где

, .

Рис. 1.13. Линейные многошаговые схемы в плоскости неопределённых коэффициентов ( K = 2 )

Условия |q1(s, a0, b0)|≤1, |q2(s, a0, b0)|≤1 ограничивают область A-устойчивых схем. На рис. 1.13 эта область в плоскости свободных коэффициентов {a0, b0} заштрихована (четырехугольник A1, A2, A3, A4). Сплошная прямая B1B2 — множество схем третьего порядка аппроксимации — не пересекается с множеством A-устойчивых схем. Точка B2 на этой прямой — единственная в этом случае схема четвертого порядка аппроксимации (49). И в общем случае доказано (теорема Далквиста), что линейных многошаговых A-устойчивых схем с порядком аппроксимации выше второго не существует.

Точка A6 с координатами

a0 = 1/3, b0 = 0 (53)

— единственная в этом случае L-устойчивая схема (схема Кёртисса–Хиршфельдера). Область жестко-устойчивых схем, очевидно, содержится внутри себя множество A-устойчивых схем. В частности, схема, соответствующая точке A6, является жестко-устойчивой.

Точка A1 с коэффициентами

a0 = b0 = 0, (54)

т. е. одношаговая линейная схема (метод трапеций), из числа A-устойчивых схем является оптимальной в том смысле, что она среди других A-устойчивых схем при K = 2 ближе всего расположена и к L-устойчивой схеме (точке A6), и к схемам 3–го порядка (наиболее точная из A-устойчивых схем).

Отметим еще схему 3-го порядка точности, соответствующую точке B1. Она является жестко-устойчивой и наиболее близка к множеству A-устойчивых схем. Для этой схемы

a0 = 0, b0 = –1/

Явным схемам на рис. 1.13 соответствует штриховая прямая. В частности, точка A5 с коэффициентами

a0 = 0, b0 = –1/2 (56)

— явная схема Адамса. Все явные схемы, как видно, не являются A-устойчивыми. А в целом, семейство линейных многошаговых схем существенно беднее схем типа Рунге–Кутты.

1.5. Схемы для продолженных систем
(схемы Обрешкова)

1.5.1. Алгоритм и аппроксимация

Лучшие черты методов Рунге–Кутты и линейных многошаговых методов сочетают схемы для продолженной системы (19) или линейные многошаговые схемы с использованием 2-й (схемы Обрешкова) и более высоких производных:

vtt(t, v) = ft = ¶f/¶t + fv f. (57)

Добавляя в (43) линейную комбинацию вторых производных vtt с дополнительными неопределенными коэффициентами сk, получим для этих схем (aK = 1)

,

(58)

.

Из начальных условий в (18) и из (57) видно, что vtt(0, v0) фактически можно считать известным.

Условия аппроксимации, как и в случае линейных многошаговых методов, являются линейными (относительно неопределенных коэффициентов ak, bk, ck) уравнениями:

,  (59)

(обеспечивают 1-й порядок аппроксимации),

(60)

(обеспечивает 2-й порядок аппроксимации вместе с (59)),

…

(61)

(обеспечивает n-й порядок аппроксимации на решениях (18) вместе с предшествующими условиями).

В случае неявных схем (bK ¹ 0, cK ¹ 0) нелинейная относительно vn+k система (58) с ограничениями (59) – (61):

решается каким-либо итерационным методом, например,

s = 0,1, …,  , (62)

B = ¶R/¶vn+k = E – tbkfv – t2ck[¶(¶f/¶t)/¶v + fvfv + C].

В (62) С — матрица, столбцами которой являются векторы

(¶(fv)/¶v1)f, (¶(fv)/¶v2)f, …, (¶(fv)/¶vK)f. (63)

1.5.2. Устойчивость

Так как свободных параметров у этого класса разностных схем больше, появляется возможность при меньшем числе точек сеточного шаблона (т. е. при меньших K ) строить схемы более высокого порядка аппроксимации, чем в случае линейных многошаговых схем. В частности, при K = 1 (одношаговые, как методы Рунге–Кутты, схемы) имеем для схем второго порядка аппроксимации

a1 = 1, a0 = –1, b0 = 1/2+c0+c1, b1 = 1/2–c0–c1, (64)

c0, c1 — произвольны. Если в дополнение к (64)

c1 = c0 – 1/6, (65)

то имеем однопараметрическое семейство 3-го порядка, а при

c0 = 1/12, c1 = –1/12 (66)

— единственную на данном шаблоне схему 4-го порядка.

Для тестового уравнения с f = lv имеем

vt = f = lv, vtt = ft = lv = l2v, …,

,

т. е. то же, что и в случае линейных многошаговых схем. В частности, при K = 1, как и в методах Рунге–Кутты, получаем геометрическую прогрессию:

vn+1 = qvn, где q = –a0/a1 = [1 + s(1/2 + c0 + c1) +

+ s2c0] / [1 – s(1/2 – c0 – c1) – s2c1]. (67)

Из условия устойчивости |q(s, c0 ,c1)| ≤ 1, требуя его выполнения для всех значений s: –¥ ≤ Re(s) ≤ 0, в плоскости свободных параметров {c0, c1} можно получить все множество A-устойчивых схем (заштриховано на рис. 1.14; вертикальная штриховка — монотонные схемы, для которых, в случае действительных значений s 0 ≤ q ≤ 1, горизонтальная штриховка — устойчивые, но не монотонные схемы, для которых 1 ≤ q < 0).

Определяя из (67) |q(s, c0, c1)| = –c0/c1, и приравнивая его нулю, получим, что множество L-устойчивых схем расположено на прямой с0 = 0 (исключая точку с1 = 0).

Схемам 3-го порядка аппроксимации (65) на рис. 1.14 соответствует прямая B1B2B3, единственной схеме четвертого порядка (66) — точка B3, расположенная на границе A-устойчивых схем.

Точка O с координатами

с0 = 0, с1 = 0 (68)

соответствует известной схеме «трапеций».

Рис. 1.14. Схемы Обрешкова в плоскости неопределённых коэффициентов ( K = 1 )

В целом данный класс схем заметно богаче класса линейных многошаговых методов. Он содержит целую полупрямую B1B2 A-устойчивых схем третьего порядка аппроксимации, в том числе L-устойчивую схему — точку B2 с координатами

c0 = 0, c1 = –1/6, (69)

чего нет даже в методах Рунге–Кутты при K = 2.

1.6. Контрольные вопросы

1.6.1. Общие вопросы к лабораторным работам 1–3

·  Жесткие системы ОДУ. Поведение решений. Примеры.

·  Жесткие системы ОДУ. Расположение собственных значений матрицы (Якоби). Признаки жёсткости.

·  Требования к численным методам решения жёстких систем.

·  Сходимость разностных схем в линейном приближении.
А- и L-устойчивые, жёстко устойчивые, монотонные схемы.

·  Построение схем методом неопределённых коэффициентов. Графическая интерпретация свойств различных методов.

·  Какие достоинства и недостатки имеют неявные схемы? схемы высокого порядка аппроксимации?

·  В каких случаях и зачем при реализации методов решения систем ОДУ используется матрица Якоби правой части?

·  Достоинства и недостатки трёх изучаемых типов схем.

1.6.2. Схемы Рунге–Кутты (работа №1)

·  Алгоритм явных, полуявных и неявных методов Рунге–Кутты.

·  Условия аппроксимации. Сколько свободных коэффициентов имеют неявные (полуявные) схемы Рунге–Кутты 2-го (3-го) порядка при K=1, 2 (K — число членов в выражении для схемы)?

·  Операторный метод исследования устойчивости. Какие методы (по признакам устойчивости) есть среди полуявных и неявных схем 2-го и 3-го порядка?

1.6.3. Линейные многошаговые схемы (работа №2)

·  Условия аппроксимации. Сколько свободных коэффициентов имеют явные (неявные) линейные многошаговые схемы 2-го, 3-го и 4-го порядка при K=1, 2, 3 (K — число точек шаблона)?

·  Исследование устойчивости. Многочлен устойчивости.

1.6.4. Схемы для продолженных систем (работа №3)

·  Возможности построения A- и L-устойчивых схем Рунге‑Кутты, линейных многошаговых схем и схем Обрешкова. Влияние явности и порядка точности на эти возможности.

2. Примеры жёстких систем ОДУ

2.1. Модель Филда–Нойса «орегонатор»

Простейшая модель периодической химической реакции Белоусова–Жаботинского состоит из трех уравнений:

На то, что система жесткая, указывают большие различия в константах скоростей реакций — есть процессы «быстрые», и есть «медленные».

Так как переменные системы — концентрации (HBrO2, Br– и Ce(IV) соответственно), то начальные условия для системы следует выбирать положительными. Как правило, их выбирают достаточно близкими к 0. Конечное время интегрирования системы Tk = 800. О системе см., например, [1–3].

2.2. Уравнение Ван-дер-Поля

Типичным примером жесткой задачи малой размерности является уравнение Ван-дер-Поля [1, 2, 4, 5]. Его возможно записать в виде системы

(1)

или в виде

, (2)

(представление Льенара). Считаем, что параметр a — большой. В расчетах рассмотреть два случая: a = 103 и a = 106. Для тестов обычно полагают Конечное время интегрирования системы, записанной в виде (1), Tk = 20.

Периодические решения жестких систем ОДУ иногда называют релаксационными автоколебаниями [4, 5].

Дополнительный вопрос: укажите преобразование, переводящее представление (1) в представление Льенара (2).

2.3. Система Ван-дер-Поля и траектории-утки

Рассмотрим неавтономную систему Ван-дер-Поля:

Как и в предыдущей задаче считаем, что a = 103 и a = 106, Рассмотреть численно случаи 0 < A < 1 и Tk = 200.

О траекториях-утках в системе Ван-дер-Поля см. [5] (строгое исследование) и [6] (популярное изложение).

2.4. Суточные колебания озона в атмосфере

Рассмотрим простейшую математическую модель колебаний концентрации озона в атмосфере [7]. Она описывается следующей неавтономной системой ОДУ:

В данной модели уравнения описывают изменение концентрации атомарного кислорода, молекулярного кислорода и озона соответственно. Считается, что изменения концентрации молекулярного кислорода невелики. Начальные значения для задачи таковы: (см–3), (см–3), (см–3), а значения констант скоростей химических реакций , .

Две другие химические реакции зависят от локальной освещенности участка земной поверхности и приближаются следующим выражением:

где с–1, с3 = 22,62, с4 = 7,601. Значения констант скоростей обращаются в ноль «ночью», резко возрастают «на рассвете», достигают максимума «в полдень» и падают до ноля «на закате». Конечное время интегрирования Tk = 172 800 c (двое суток). Данная система является жесткой «ночью» и умеренно жесткой «в светлое время суток».

2.5. Уравнение Бонгоффера–Ван-дер-Поля

Рассмотрим еще один пример жесткой задачи малой размерности, имеющей периодическое решение [5, 8]:

Здесь a = 103 и a = 106, Уравнение описывает протекание тока через клеточную мембрану. Постоянная компонента тока c в безразмерной записи системы такова, что 0 < c <1, b > 0. Tk = 20.

2.6. Сингулярно-возмущенная система — модель двухлампового генератора Фрюгауфа

Система более высокой размерности, имеющая решение в виде релаксационного цикла, приведена в [4] (см. также [8]). Она имеет вид:

Здесь α > 0 порядка 1, функция Tk = 20, ε = 10–3, 10–6.

2.7. Простейшая модель гликолиза

Простейшая модель гликолиза описывается уравнениями следующего вида [8]:

предложенными Дж. Хиггинсом. В системе β = 10, α = 100, 200, 400, 1000. Начальные условия для системы: y1(0) = 1, y2(0) = 0,001, Tk = 50. Решение этой системы — релаксационные автоколебания (жесткий предельный цикл).

2.8. Модель химических реакций Робертсона

Один из первых и самых популярных примеров «жесткой» системы ОДУ принадлежит Робертсону (1966) и имеет вид, типичный для моделей химической кинетики — в правой части системы стоят полиномы второй степени от концентраций (сравните с орегонатором). Система Робертсона имеет вид [1]:

Начальные условия для системы таковы: y1(0) = 1, y2(0) = 0, y3(0) = 0. Рассматриваются следующие величины отрезка интегрирования: Tk = 40 (в работе Робертсона рассматривался именно такой отрезок интегрирования), Tk = 100, 1000, …, 1011. О свойствах задачи см. в [1].

2.9. Модель дифференциации растительной ткани

Данный пример из [1] — типичный случай биохимической модели «умеренной» размерности (современные модели, например, фотосинтеза включают сотни уравнений подобного типа). Хотя данная модель является умеренно жесткой, тем не менее, ее лучше решать с помощью методов, предназначенных для решения ЖС ОДУ.

Начальные значения всех переменных системы равны 0, кроме y1(0) = 1 и y8(0) = 0,0057. Длина отрезка интегрирования Tk = 421,8122.

2.10. Задача E5

Еще одна модель химической реакции из [1], получившая свое название Е5 в более ранних публикациях.

Начальные условия: y1(0) = 1,76·10–3, а все остальные переменные равны 0. Значения коэффициентов модели следующие: A = 7,89·10–10, B = 1,1·107, C = 1,13·103, M = 106. Первоначально задача ставилась на отрезке Tk = 1000, но впоследствии было обнаружено, что она обладает нетривиальными свойствами вплоть до времени Tk = 1013 (см. [1]).

Обратите особое внимание, что в процессе расчетов приходится иметь дело с очень малыми концентрациями реагентов (малы значения y2, y3 и y4).

2.11. Уравнение Релея

Уравнение Релея во многом похоже на уравнение Ван-дер-Поля [8]. Рассматривается задача вида

.

Решить задачу, записав уравнение Релея в виде системы ОДУ. Начальные условия: x(0) = 0, μ = 1000, Tk = 1000.

2.12. Экогенетические модели

Рассмотрим пример системы уравнений, которая описывает изменения численности популяций двух видов и эволюцию некого генетического признака α. Система имеет вид:

Параметры задачи таковы: ε ≤ 0,01, 0 ≤ x0 ≤ 3, 0 ≤ y0 ≤ 15, α0 = 0,01, Tk = 1500. Наличие малого параметра в третьем уравнении системы показывает, что генетический признак меняется медленнее, чем численность популяций. Решение системы — релаксационные колебания.

Наряду с данной задачей в [9] описана более интересная: численность двух популяций зависит от их взаимодействия и двух медленно меняющихся генетических признаков.

Параметры задачи таковы: ε ≤ 0,01, 0 ≤ x0 ≤ 40, 0 ≤ y0 ≤ 40, α10 = 0, α20 = 10, Tk = 2000. Рассмотреть также модификацию предыдущей системы [9]:

Параметры задачи: ε ≤ 0,01, 0 ≤ x0 ≤ 40, 0 ≤ y0 ≤ 40, α10 = 0, α20 = 10, Tk = 2000.

Список литературы

1.  Хайрер Э., Ваннер Г. Решение обыкновенных дифференциальных уравнений. Жесткие и дифференциально–алгебраические задачи. — М.: Мир, 1999. — 685 с.

2.  Хайрер Э., Нерсетт С., Ваннер Г. Решение обыкновенных дифференциальных уравнений. Нежесткие задачи. — М.: Мир, 1990. — 512 с.

3.  Колебания и бегущие волны в химических системах. / Под ред. Р. Филда, М. Бургер. — М.: Мир, 1988. — 720 с.

4.  Мищенко Е. Ф., Розов Н. Х. Дифференциальные уравнения с малым параметром и релаксационные колебания. — М.: Наука, 1975. — 248 с.

5.  Мищенко Е. Ф., Колесов Ю. С., Колесов А. Ю., Розов Н. Х. Периодические движения и бифуркационные процессы в сингулярно-возмущенных системах — М.: Наука. Физматлит, 1995. — 336 с.

6.  Малинецкий Г. Г. Хаос. Структуры. Вычислительный эксперимент. — М.: Наука, 1998.

7.  Каханер Д., Моулер К., Нэш С.. Численные методы и программное обеспечение. — М.: Мир, 1998. — 575 с.

8.  Ланда П. С. Нелинейные колебания и волны. — М.: Наука. Физматлит, 1997. — 496 с.

9.  Кондрашов А. С, Хибник А. И. Экогенетические модели как быстро-медленные системы. // Исследования по математической биологии. — Пущино, 1996. — с. 88–123.

10.  Ракитский Ю. В., Устинов С. М., Черноруцкий И. Г. Численные методы решения жестких систем. — М.: Наука, 1979.

11.  Холл Дж., Уатт Дж. (ред.). Современные численные методы решения обыкновенных дифференциальных уравнений. — М.: Мир, 1979.