авторефераты диссертаций БЕСПЛАТНАЯ БИБЛИОТЕКА РОССИИ

КОНФЕРЕНЦИИ, КНИГИ, ПОСОБИЯ, НАУЧНЫЕ ИЗДАНИЯ

<< ГЛАВНАЯ
АГРОИНЖЕНЕРИЯ
АСТРОНОМИЯ
БЕЗОПАСНОСТЬ
БИОЛОГИЯ
ЗЕМЛЯ
ИНФОРМАТИКА
ИСКУССТВОВЕДЕНИЕ
ИСТОРИЯ
КУЛЬТУРОЛОГИЯ
МАШИНОСТРОЕНИЕ
МЕДИЦИНА
МЕТАЛЛУРГИЯ
МЕХАНИКА
ПЕДАГОГИКА
ПОЛИТИКА
ПРИБОРОСТРОЕНИЕ
ПРОДОВОЛЬСТВИЕ
ПСИХОЛОГИЯ
РАДИОТЕХНИКА
СЕЛЬСКОЕ ХОЗЯЙСТВО
СОЦИОЛОГИЯ
СТРОИТЕЛЬСТВО
ТЕХНИЧЕСКИЕ НАУКИ
ТРАНСПОРТ
ФАРМАЦЕВТИКА
ФИЗИКА
ФИЗИОЛОГИЯ
ФИЛОЛОГИЯ
ФИЛОСОФИЯ
ХИМИЯ
ЭКОНОМИКА
ЭЛЕКТРОТЕХНИКА
ЭНЕРГЕТИКА
ЮРИСПРУДЕНЦИЯ
ЯЗЫКОЗНАНИЕ
РАЗНОЕ
КОНТАКТЫ


Pages:     | 1 | 2 || 4 | 5 |   ...   | 9 |

«Е.А. Новиков, Ю.В. Шорников КОМПЬЮТЕРНОЕ МОДЕЛИРОВАНИЕ ЖЕСТКИХ ГИБРИДНЫХ СИСТЕМ Е.А. Новиков, Ю.В. Шорников КОМПЬЮТЕРНОЕ ...»

-- [ Страница 3 ] --

Будем предполагать, что ее коэффициенты удовлетворяют соотноше ниям g3 (3 2 ) = 32 2 (2 3 2 ) / 6, (4.4) 232 p3 = g, 2 p2 + 3 p3 = 1 / 2, p1 + p2 + p3 = 1, (4.5) при которых она имеет второй порядок точности, а ее локальная ошиб ка n,2 имеет вид n,2 = (1 6 g )h3 f 2 f / 6 + O (h 4 ), 2 = 21, 3 = 31 + 32. Тогда для контроля точности (4.3) при условиях (4.4), (4.5) можно применять неравенство || k2 k1 || 6 | 2 | / | 1 6 g |. (4.6) Получим неравенство для контроля устойчивости. Для этого с ис пользованием (3.33) составим линейную комбинацию векторов ki, 1 i 3, и hf ( yn+1 ) с неопределенными коэффициентами ci, 1 i 4, т. е.

4 ci ki + c4 hf ( yn+1 ) = ci hf n + ( 2c2 + 3c3 + c4 )h2 f n fn + i = i = ( ) + 0,5 2 c2 + 3 c3 + c4 h3 f n f n2 + ( 232 c3 + 0,5c4 )h3 f n2 f n + O (h 4 ).

2 Отсюда, полагая c4 свободным параметром, имеем, что при l = первое и второе условия (4.1) будут выполнены, если ci, 1 i 4, удовлетворяют соответственно следующим системам алгебраических уравнений:

Г л а в а 4. АЛГОРИТМЫ ИНТЕГРИРОВАНИЯ С КОНТРОЛЕМ УСТОЙЧИВОСТИ 1 c 1 1 2 c = c 1, 0 2 3 2 4 1 / 0 c 232 (4.7) 1 c1 1 1 2 3 c2 = c4 1.

0 3 c 2 Обозначим через bi и bi, 1 i 4, соответствующие решения (4.7).

Тогда для контроля устойчивости (4.3) можно использовать нера венство 3 max biki + b4 hf ( yn +1 ) 1 j N i =1 j 232b3 + 0,5b D, (4.8) 2b2 + 3b3 + b 3 bi ki + b4 hf ( yn +1 ) i =1 j где положительная постоянная D связана с размером области устой чивости схемы (4.3). Например, если собственные числа матрицы Яко би задачи (2.1) вещественные и отрицательные, а метод (4.3) использу ется с набором коэффициентов 1 1 g=, 2 = 21 =, 31 = 32 =, 16 3 7 1 3 3 =, p1 =, p2 =, p3 =, 9 7 8 при которых длина интервала устойчивости равна шести, то D 6.

На рис. 4.1–4.4 слева приведены линии уровня | Q3,2 ( z ) | = s при значениях s, равных числам 1;

0,7;

0,3;

0,1;

справа – объемные изо бражения областей устойчивости. На рис. 4.1 показана область устой чивости численной формулы (4.3).

4.1. Схемы второго порядка точности Рис. 4.1. Область устойчивости метода (3.40), (3.53) Из (4.4) и (4.5) следует, что при выборе коэффициентов численной схемы (4.3) имеется некоторый произвол, воспользовавшись которым, можно упростить неравенство (4.8). Выберем для определенности в формулах (4.4), (4.5) значение параметра g равным 1/16 и положим 2 = 3 и 32 = 1 / 3. Тогда коэффициенты схемы (4.3) однозначно оп ределяются из (4.4), (4.5) и имеют вид 2 1 1 15 2 = 3 = 21 =, 31 = 32 =, p1 =, p2 =, p3 =, (4.9) 3 3 4 32 при этом длина интервала устойчивости (4.3), (4.9) равна шести.

При получении неравенства для контроля устойчивости применим к разности (k3 k2 ) формулу Тейлора с остаточным членом в форме Лагранжа первого порядка, т. е.

2 1 f ( y n) 1 ( k2 k1 ), k3 k2 = h f yn + k1 + k2 f yn + k1 = h y 3 3 где y n вычислено в некоторой окрестности решения y (tn ). Учитывая, что k2 k1 = 2h 2 f n f n / 3 + O(h3 ), для контроля устойчивости схемы (4.3) с параметрами (4.9) можно применять следующее неравенство:

3 max ( k3 k2 )i / (k2 k1 )i 6. (4.10) 1i N Из многочисленных расчетов задач умеренной жесткости следует, что неравенство (4.10) не хуже (4.8) с точки зрения эффективности ал горитма интегрирования. Поэтому далее при изучении других числен ных схем будем поступать по аналогии с получением (4.10).

Г л а в а 4. АЛГОРИТМЫ ИНТЕГРИРОВАНИЯ С КОНТРОЛЕМ УСТОЙЧИВОСТИ 4.2. СХЕМА ТРЕТЬЕГО ПОРЯДКА ТОЧНОСТИ Рассмотрим четырехстадийную схему 2 1 1 yn+1 = yn + k1 + k2 + 72 g k3 + 6 k4, 6 3 72 g k1 = hf ( yn ), k2 = hf ( yn + 0,5k1 ), (4.11) k3 = hf ( yn + (0,5 12 g ) k1 + 12 gk2 ), k4 = hf ( yn + k3 ) и вспомогательный метод вида zn +1 = yn + k1 / 6 + 2k2 / 3 + k4 / 6, (4.12) где g – постоянная. С использованием (3.33) нетрудно убедиться в том, что порядок точности численных формул (4.11) и (4.12) равен со ответственно трем и двум, а их локальные ошибки n,3 и n,2 имеют вид 1 24 g 4 h f f + O(h 4 ), n,3 = h f f + O(h5 ).

n, 2 = 12 Учитывая, что yn+1 zn+1 = g 1 (k3 k2 ) / 72, неравенство для контроля точности схемы (4.11) можно записать в виде || k3 k2 || 144 g / (1 24 g ).

Теперь получим неравенство для контроля устойчивости схемы (4.11). Для этого применим к разности (k3 k2 ) формулу Тейлора с остаточным членом в форме Лагранжа первого порядка. В результате получим f ( y n) k3 k2 = 12 gh (k2 k1 ), y где y n вычислено в некоторой окрестности решения y (tn ). Учитывая, что k2 k1 = h 2 f n f n / 2 + O (h3 ), для контроля устойчивости схемы (4.11) можно использовать неравенство типа (4.10), т. е.

max (k3 k2 )i / (k2 k1 )i D / | 12 g |, (4.13) 1i N 4.3. Схемы четвертого и пятого порядков точности где постоянная D связана с размером области устойчивости схемы.

Известно, что максимальная длина интервала устойчивости четырех стадийной схемы третьего порядка точности равна 6,0273, причем она достигается при g = 0,018456. Выберем данное значение g и поло жим D = 6 в формуле (4.13).

4.3. СХЕМЫ ЧЕТВЕРТОГО И ПЯТОГО ПОРЯДКОВ ТОЧНОСТИ Как показывают расчеты, контроль устойчивости численной схемы не только приводит к повышению эффективности алгоритма интегри рования, но и позволяет вычислить решение ряда задач, которые не удается рассчитать алгоритмом без контроля устойчивости. Поэтому возникает вопрос о том, как организовать контроль устойчивости из вестных схем, хорошо зарекомендовавших себя при решении практи ческих задач, не внося существенных изменений в вычислительный процесс. К наиболее известным и удачным явным формулам типа Рун ге–Кутты относятся пятистадийная схема четвертого порядка точности 1 2 yn+1 = yn + k1 + k4 + k5, k1 = hf ( yn ), 6 3 1 k2 = hf yn + k1, k3 = hf yn + k1 + k2, (4.14) 3 3 1 1 k4 = hf yn + k1 + k3, k5 = hf yn + k1 k3 + 2k4, 8 2 построенная Р. Мерсоном, и шестистадийная схема пятого порядка точности, полученная Э. Фельбергом. Одна из причин их успеха за ключается, по-видимому, в экономичном способе оценки ошибки, на основе которой осуществляется контроль точности вычислений. Кроме того, области устойчивости обоих методов достаточно велики как по вещественной, так и по мнимой оси комплексной плоскости z = h.

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

Г л а в а 4. АЛГОРИТМЫ ИНТЕГРИРОВАНИЯ С КОНТРОЛЕМ УСТОЙЧИВОСТИ Сначала рассмотрим метод Мерсона, для контроля точности кото рого можно применять неравенство || 2k1 9k3 + 8k4 k5 || /30 55/ 4. (4.15) Функция устойчивости Q5,4 ( z ) численной формулы (4.14) имеет вид Q5,4 ( z ) = 1 + z + z 2 / 2 + z 3 / 6 + z 4 24 + z 5 / 144, z = h.

На рис. 4.2 показана область устойчивости метода Мерсона.

Рис. 4.2. Область устойчивости метода Мерсона Из (4.14) видно, что для данной численной схемы выполняется ус ловие 2 = 3, где 2 = 21, 3 = 31 + 32. Поэтому неравенство для контроля устойчивости получим по аналогии с (4.10). Применяя к раз ности (k3 k2 ) формулу Тейлора с остаточным членом в форме Ла гранжа первого порядка, будем иметь 1 f ( y n) k3 k2 = h ( k2 k1 ), y где вектор y n вычислен в некоторой окрестности решения y (tn ). Учи тывая, что k2 k1 = h 2 f n f n / 3 + O(h3 ), для контроля устойчивости (4.14) можно использовать неравенство 6 max (k3 k2 )i / (k2 k1 )i 3,5, (4.16) 1i N 4.3. Схемы четвертого и пятого порядков точности где числу 3,5 примерно равна длина интервала устойчивости числен ной формулы (4.14). Отметим, что по мнимой оси область устойчиво сти метода (4.14) также ограничена числом 3,5.

Теперь рассмотрим шестистадийный метод Рунге–Кутты–Фель берга, который имеет вид yn+1 = yn + pli ki, (4.17) i = где 1 3 k1 = hf ( yn ), k2 = hf yn + k1, k3 = hf yn + k1 + k2, 4 1932 k4 = hf yn + k1 k2 + k3, (4.18) 2197 439 3680 k5 = hf yn + k1 8k2 + k3 k4, 216 8 3544 k6 = hf yn k1 + 2k2 k3 + k4 k5.

27 2565 В (4.17) первый индекс при записи коэффициентов pli, 1 i 6, озна чает порядок точности метода. При значениях коэффициентов 16 p51 =, p52 = 0, p53 =, 135 (4.19) 28561 9 p54 =, p55 =, p56 = 56430 50 формула (4.17) совпадает с методом Фельберга и имеет пятый порядок точности. С использованием (3.33) получим, что ее локальная ошибка n,5 представима в виде n,5 = 17 h6 f 5 f / 18720 + O(h6 ), где элемен тарный дифференциал вычисляется на точном решении y (tn ).

Функция устойчивости Q6,5 ( z ), z = h, метода Фельберга имеет вид Q6,5 ( z ) = 1 + z + z 2 / 2 + z 3 / 6 + z 4 / 24 + z 5 / 120 + z 6 / 2080, Г л а в а 4. АЛГОРИТМЫ ИНТЕГРИРОВАНИЯ С КОНТРОЛЕМ УСТОЙЧИВОСТИ а область устойчивости показана на рис. 4.3. Известны коэффициенты 25 p41 =, p42 = 0, p43 =, 16 (4.20) 2197 p44 =, p45 =, p46 = 0, 4104 при которых (4.17) имеет четвертый порядок, а ее локальная ошибка имеет вид n,4 = h5 f 4 f / 780 + O (h5 ). (4.21) Рис. 4.3. Область устойчивости метода Фельберга Функция устойчивости Q5,4 ( z ), z = h, метода (4.18), (4.20) имеет вид Q5,4 ( z ) = 1 + z + z 2 / 2 + z 3 / 6 + z 4 / 24 + z 5 / 104, а область устойчивости показана на рис. 4.4.

Рис. 4.4. Область устойчивости метода (4.18), (4.20) 4.3. Схемы четвертого и пятого порядков точности Соотношение (4.21) получено с применением (3.33). Так как в каж дой точке имеем два приближения к решению, вычисленных методами четвертого и пятого порядков, то для контроля точности можно ис пользовать оценку ошибки n,5 вида n,5 = 17 ( p5i p4i )ki / 24, где i = ki, 1 i 6, определены формулами (4.18). В результате для контроля точности вычислений будем применять следующее неравенство:

( p5i p4i )ki / 24, 17 (4.22) i = где – требуемая точность расчетов. Значения параметров i, 2 i 6, определяемых по формуле (3.20), имеют вид 2 = 1 / 4, 3 = 3 / 8, 4 = 12 / 13, 5 = 1, 6 = 1 / 2. (4.23) Из (4.23) видно, что приращение k5 вычисляется в точке tn+1, и по этому нет необходимости в дополнительном неравенстве для выбора величины шага.

Теперь перейдем к построению неравенства для контроля устойчи вости. Отметим, что так как для схемы (4.17) не выполняется условие 2 = 3, то нельзя построить неравенство типа (4.10). Поэтому при получении оценки максимального собственного числа матрицы Якоби ограничимся рассуждениями на линейной задаче с постоянными ко эффициентами y = Ay + b. (4.24) Запишем ki, 1 i 3, применительно к задаче (4.24), т. е.

k1 = h( Ayn + b), k2 = k1 + (hA) k1 / 4, k3 = k1 + 3(hA)k1 / 8 + 9(hA)2 k1 / 128.

Составляя из данных стадий соотношения типа (4.1), получим k2 k1 = (hA)k1 / 4, 32k3 / 9 16k2 / 3 + 16k1 / 9 = (hA) 2 k1 / 4.

В результате для контроля устойчивости (3.3) можно применять нера венство max (32k3 48k2 + 16k1 )i / (k2 k1 )i / 9 3,6. (4.25) 1i N Г л а в а 4. АЛГОРИТМЫ ИНТЕГРИРОВАНИЯ С КОНТРОЛЕМ УСТОЙЧИВОСТИ В данном неравенстве используется число 3,6, которое примерно равно длине интервала устойчивости схемы (4.17) с коэффициентами (4.19). Отметим, что по мнимой оси область устойчивости также огра ничена числом 3,6.

Заметим, так как в случае нелинейной задачи в представлениях ki, 1 i 3, в виде рядов Тейлора имеются дополнительные слагаемые, т. е.

k2 k1 = h 2 f n f n / 4 + h3 f n f n2 / 32 + O (h 4 ), (4.26) (32k3 48k2 + 16k1 ) / 9 = h3 f n 2 f n / 4 + 4h3 f n f n2 / 3 + O(h 4 ), то они будут вносить дополнительные искажения в оценку максималь ного собственного числа. В этом смысле (4.25) грубее (4.16). Однако результаты расчетов показывают, что использование неравенства (4.25) совместно с выбором шага по формуле вида (2.48) приводит к повышению эффективности по сравнению с алгоритмом без контроля устойчивости.

Глава АЛГОРИТМЫ ИНТЕГРИРОВАНИЯ ПЕРЕМЕННОГО ПОРЯДКА И ШАГА Н а основе явных схем типа Рунге–Кутты с расширенными об ластями устойчивости разработаны алгоритмы переменного порядка и шага для решения умеренно жестких задач. Выбор метода осуществляется по точности и устойчивости соответствующей числен ной схемы. Эффективность алгоритмов достигается за счет примене ния на каждом участке наиболее эффективной численной формулы.

5.1. АЛГОРИТМ НА ОСНОВЕ ТРЕХСТАДИЙНОЙ СХЕМЫ Изучим вопросы повышения эффективности алгоритмов интегри рования за счет использования на различных участках численных схем с различными областями устойчивости. При этом будем предполагать, что все численные формулы основаны на одних и тех же стадиях, что позволяет достаточно просто организовать переключение с одной схе мы на другую. Итак, рассмотрим трехстадийную численную формулу типа Рунге–Кутты следующего вида:

yn +1 = yn + p1k1 + p2 k2 + p3k3, k1 = hf ( yn ), (5.1) 2 k2 = hf yn + k1, k3 = hf yn + k1 + k2.

3 Данный метод был подробно изложен в предыдущей главе.

Метод второго порядка Нетрудно видеть, что схема (5.1) будет иметь второй порядок точ ности, если коэффициенты pi, 1 i 3, вычислить по формулам p1 = 1 / 4, p2 = (3 18 g ) / 4, p3 = 9 g / 2, (5.2) Г л а в а 5. АЛГОРИТМЫ ИНТЕГРИРОВАНИЯ ПЕРЕМЕННОГО ПОРЯДКА И ШАГА где от параметра g зависят размер и форма области устойчивости.

Функция роста Q3,2 ( z ) формулы (5.1) имеет вид Q3,2 ( z ) = = 1 + z + z 2 / 2 + gz 3. В дальнейшем для определенности будем рас сматривать методы с максимальным интервалом устойчивости, что соответствует значению g, равному 1 / 16. На рис. 5.1–5.33 слева при ведены линии уровня | Q3,2 ( z ) |= s при s, равном 1;

0,7;

0,3;

0,1;

справа – объемные изображения областей устойчивости. Область устойчивости численной схемы (5.1), (5.2) при g = 1 / 16 показана на рис. 5.1.

Рис. 5.1. Область устойчивости метода (5.1), (5.2) при g = 1/ Введем обозначения:

|1 6g | |1 6g | An = || k2 k1 ||, An = || hf ( yn +1 ) k1 ||, 4 | ( k3 k2 )i | Vn = 3 max, (5.3) 1i N | ( k2 k1 )i | где || || – некоторая норма в R N. Тогда для контроля точности и при выборе величины шага интегрирования схемы (5.1) с коэффициентами (5.2) можно применять неравенства An, An, (5.4) а для контроля устойчивости – следующее:

Vn D. (5.5) 5.1. Алгоритм на основе трехстадийной схемы При g = 1 / 16 постоянную величину D можно выбрать равной 6, по тому что длина интервала устойчивости схемы (5.1) с коэффициентами (5.2) примерно равна 6. При практической реализации алгоритма ин тегрирования во второй формуле (5.3) обычно используется постоян ная | 1 6 g | /4 вместо | 1 6 g | /6 с целью уменьшения возможных по вторных вычислений решения.

Метод первого порядка Перейдем к построению формулы интегрирования более грубой в смысле точности, но с более широкой областью устойчивости. Для оп ределенности снова будем рассматривать метод с максимальным ин тервалом устойчивости. Для построения такой численной схемы при меним (5.1) для решения линейного тестового уравнения y = y, y (0) = y0, Re( ) 0. Будем иметь yn+1 = Q3,1 ( z ) yn, z = h, Q3,1 ( z ) = 1 + ( p1 + p2 + p3 ) z + 2( p2 + p3 ) z 2 / 3 + 0, 4 p3 z 3.

Область устойчивости (5.1) ограничена кривой | Q3,1 ( z ) | = 1 в ком плексной плоскости z = h. Из вида функции устойчивости Q3,1 ( z ) следует, что конфигурация области устойчивости схемы (5.1) зависит от коэффициентов pi, 1 i 3. Целью является получение такого на бора коэффициентов p1, p2 и p3, чтобы схема (5.1) имела макси мальный интервал устойчивости и первый порядок точности. Второе условие будет выполнено, если p1 + p2 + p3 = 1. В результате требует ся построить полином вида Q3,1 ( z ) = 1 + z + c2 z 2 + c3 z 3, (5.6) для которого неравенство | Q3,1 ( z ) | 1 выполняется на максимальном интервале c2 = 2( p2 + p3 ) / 3, c3 = 2 p3 / 5. (5.7) Г л а в а 5. АЛГОРИТМЫ ИНТЕГРИРОВАНИЯ ПЕРЕМЕННОГО ПОРЯДКА И ШАГА Рассмотрим одновременно многочлен Чебышева T3 ( x) третьей степени T3 ( x) = 4 x3 3x, (5.8) для которого выполнено неравенство | T3 ( x) | 1 при x [1,1]. Прове дем замену переменных в (5.8), полагая x = 1 2 z /. Получим T3 ( z ) = 1 18 z / + 48 z 2 / 2 32 z 3 / 3, (5.9) при этом отрезок [, 0] отображается на отрезок [ 1,1]. Нетрудно пока зать, что среди всех полиномов вида (5.6) для (5.9) при = 18 нера венство | T3 ( z ) | 1 выполняется на максимальном интервале [, 0].

Потребуем совпадения коэффициентов (5.6) и (5.9), т. е.

2( p2 + p3 ) / 3 = 48 / 2, 2 p3 / 5 = 32 / 3. Отсюда и из условия первого порядка точности имеем p1 = 18 / 72 / 2, p2 = 72 / 2 + 144 / 3, p3 = 144 / 3.

Подставляя сюда = 18, получим коэффициенты схемы (5.1), т. е.

p1 = 7 / 9, p2 = 16 / 81, p3 = 2 / 81, (5.10) при которых ее область устойчивости расширена до 18 по действи тельной оси. Функция роста Q3,1 ( z ) метода (5.1) с коэффициентами (5.10) имеет вид Q3,1 ( z ) = 1 + x + 4 z 2 / 27 + 4 z 3 / 729.

Область устойчивости метода (5.1), (5.10) показана на рис. 5.2. Так как длина интервала устойчивости формулы (5.1) с коэффициентами (5.2) при g = 1 / 16 равна примерно 6, то за счет замены (5.2) на (5.10) на участке установления получим сокращение вычислительных затрат приблизительно в 3 раза.

Теперь построим неравенства для контроля точности и устойчиво сти схемы (5.1) с коэффициентами (5.10). Для этого разложим стадии ki, 1 i 3, в ряд Тейлора и подставим в первую формулу (5.1).

С применением (3.32) имеем 5.1. Алгоритм на основе трехстадийной схемы yn+1 = yn + hf n + c2 h 2 f n f n + O (h3 ), (5.11) где c2 вычисляется по первой формуле (5.7). В случае использования набора коэффициентов (5.10) имеет место c2 = 4 / 27.

Рис. 5.2. Область устойчивости метода (5.1), (5.10) Пусть yn = y (tn ). Тогда, сравнивая (3.31) и (5.11) до членов с h включительно, получим, что локальная ошибка n,1 схемы первого порядка точности имеет вид n,1 = (1 2c2 ) h 2 f f / 2 + O(h3 ). Формулу первого порядка предполагается применять на участке установления, где ошибки за счет неточности (5.1) невелики. Поэтому для контроля точности вычислений можно применять оценку локальной ошибки.

Учитывая, что k2 k1 = 2h 2 f n f n / 3 + O(h3 ), hf ( yn +1 ) k1 = h 2 f n f n + O(h3 ), для контроля точности и при выборе величины шага интегрирования будем применять соответственно неравенства dAn, dAn, d = | 3 6c2 | / |1 6 g |, (5.12) где An и A вычисляются по формулам (5.3). В случае (5.1), (5.10) имеем d = 152 / 45. Так как длина интервала устойчивости схемы (5.1), (5.10) равна 18, то для ее контроля устойчивости можно применять неравенство (5.5) при D 18. Заметим, что в неравенстве (5.5) исполь зуются постоянные 6 и 18, ограничивающие длину интервала устой чивости, в то время как для комплексных z = h их нужно уменьшить.

Поэтому при настройке на задачи они могут быть изменены. Здесь для определенности оставим их такими.

Г л а в а 5. АЛГОРИТМЫ ИНТЕГРИРОВАНИЯ ПЕРЕМЕННОГО ПОРЯДКА И ШАГА Отметим, что применение оптимальных полиномов типа (5.9) в ка честве функции устойчивости приводит к тому, что внутри интервала устойчивости всегда найдется одна или несколько точек, в которых она равна +1 или 1. Это в свою очередь означает, что соответствую щая область устойчивости будет «почти» многосвязной и небольшие возмущения могут приводить к ее сокращению. Ниже будет описан алгоритм построения многочленов, для которых длина интервалов ус тойчивости близка к максимальной, в то же время область устойчиво сти в какой-то мере лишена описанного недостатка. При этом все при веденные выше рассуждения остаются в силе, а для того, чтобы воспользоваться другим многочленом устойчивости, нужно пересчи тать коэффициенты pi, 1 i 3, и соответствующим образом подпра вить константы в неравенствах (5.4), (5.5) и (5.12). Данное замечание относится к приведенным ниже алгоритмам на основе других формул.

Алгоритм переменного порядка и шага Теперь перейдем к формулировке алгоритма интегрирования на основе численной схемы (5.1) с коэффициентами (5.2) и (5.10). Так как в (5.1), (5.2) и (5.1), (5.10) используются одни и те же стадии ki, 1 i 3, то при расчете по любой из них имеется полная информация о другой численной схеме. Благодаря этому можно сформулировать не сколько вариантов алгоритмов интегрирования с различными спосо бами выбора эффективной численной схемы на шаге интегрирования.

В частности, решение о выборе той или иной схемы можно принимать после вычисления стадий ki, 1 i 3, на основании неравенства для контроля точности (5.4), (5.12) и устойчивости (5.5). Здесь приведен пример алгоритма интегрирования, когда схема выбирается до начала вычислений на каждом шаге. Как показывают расчеты, эффективности обоих типов алгоритмов различаются незначительно. Первый шаг все гда выполняется по схеме (5.1), (5.2) как более точной. С учетом того, что имеет место An = O (h 2 ), An = O (h 2 ) и Vn = O(h), алгоритм на ос нове схемы (5.1) с коэффициентами (5.2) формулируется следующим образом.

Шаг 1а. Вычисляются k1 и k2 по формуле (5.1).

Шаг 2а. Вычисляется оценка An по первой формуле (5.3).

5.1. Алгоритм на основе трехстадийной схемы Шаг 3а. Вычисляется значение параметра q1 по формуле q1 An =.

Шаг 4а. Если q1 1, то происходит повторное вычисление с шагом q1hn, т. е. осуществляется возврат на шаг 1а.

Шаг 5а. Вычисляется k3 и приближение к решению yn +1 по фор муле (5.1) с коэффициентами (5.2).

Шаг 6а. Вычисляется оценка ошибки An по второй формуле (5.3).

Шаг 7а. Вычисляется значение параметра q2 по формуле q2 An =.

Замечание. Для повышения надежности расчетов здесь при q2 можно сделать возврат на шаг 1а с шагом интегрирования q2 hn.

Шаг 8а. Вычисляется значение Vn по последней формуле (5.3).

Шаг 9а. Вычисляются значения параметров q1, q2, r1 и r2 соответ 2 ственно по формулам q1 An = / d, q2 An = / d, r1Vn = 18 и r2Vn = 6.

Шаг 10а. Вычисляются значения hnew1 и hnew2 по формулам hnew1 = max ( hn, min(q1, q2, r1 )hn ),. (5.13) hnew2 = max ( hn, min(q1, q2, r2 )hn ) Шаг 11а. Если hnew2 hnew1, то выполняется следующий шаг ин тегрирования при hn +1 = hnew2. В противном случае hn +1 = hnew1 и управление передается алгоритму на основе численной схемы (5.1), (5.10). Переход со схемы (5.1), (5.2) осуществляется только в том слу чае, когда прогнозируемый шаг увеличивается и из условия точности и устойчивости для схемы (5.1), (5.10) становится больше шага для схе мы (5.1), (5.2).

Сформулируем алгоритм интегрирования на основе схемы (5.1), (5.10).

Шаг 1б. Вычисляются k1 и k2 по формуле (5.1).

Шаг 2б. Вычисляется оценка ошибки An по первой формуле (5.3).

Шаг 3б. Вычисляется значение параметра q1 по формуле q1 An = / d.

Г л а в а 5. АЛГОРИТМЫ ИНТЕГРИРОВАНИЯ ПЕРЕМЕННОГО ПОРЯДКА И ШАГА Шаг 4б. Если q1 1, то происходит повторное вычисление с шагом q1hn, т. е. происходит возврат на шаг 1б.

Шаг 5б. Вычисляются k3 и приближение к решению по формуле (5.1) с коэффициентами (5.10).

Шаг 6б. Вычисляется оценка ошибки An по второй формуле (5.3).

Шаг 7б. Вычисляются значение параметра q2 по формуле q2 An = / d.

Замечание. Для повышения надежности здесь при q2 1 можно сделать возврат на шаг 1б с шагом интегрирования q2 hn. Кроме того, здесь же можно вычислить приближение к решению по схеме (5.1), (5.2) как более точной, заново вычислить hf ( yn +1 ) и оценку An. В случае выполнения второго неравенства (5.4) полученное приближе ние может быть принято.

Шаг 8б. Вычисляется значение Vn по последней формуле (5.3).

Шаг 9б. Вычисляются значения параметров q1, q2, r1 и r2 по фор мулам 2 q1 An =, q2 An =, r1Vn = 18, r2Vn = 6.

Шаг 10б. Вычисляются параметры hnew1 и hnew2 по формуле (5.13).

Шаг 11б. Если hnew1 hnew2, то hn +1 = hnew1, и выполняется сле дующий шаг. В противном случае hn +1 = hnew2 и управление передает ся на схему (5.1), (5.2).

Ниже данный алгоритм переменного порядка и шага будем назы вать DISPD. В описанном алгоритме в некоторых сложных ситуациях предпочтение отдается алгоритму на основе схемы (5.1), (5.2) как бо лее точному. Заметим также, что по отдельности каждый из алгорит мов на основе (5.1), (5.2) или (5.1), (5.10) будет менее эффективен DISPD по следующей причине. Если шаг ограничен лишь требованием точности, то метод (5.1), (5.2) эффективнее (5.1), (5.10). Однако как только начинается участок установления, где шаг интегрирования вы бирается лишь из условия устойчивости, алгоритм (5.1), (5.10) сразу становится эффективнее (5.1), (5.2) в 3 раза. В DISPD оба метода при 5.2. Алгоритм с применением стадий метода Рунге–Кутты–Мерсона меняются на тех участках, где они наиболее эффективны. Если сравни вать алгоритмы на основе (5.1), (5.2) и на основе (5.1), (5.10) между собой, то ни одному из них предпочтение отдать нельзя. Если шаг ог раничен требованием точности на большей части интервала интегри рования, то, безусловно, эффективнее (5.1), (5.2). Если шаг в основном ограничен устойчивостью, то предпочтительнее алгоритм на основе (5.1), (5.10).

В алгоритме DISPD реализуется ситуация, когда на всем проме жутке интегрирования используется одна из схем: (5.1), (5.2) или (5.1), (5.10). С помощью признака легко реализуются ситуации:

• расчеты по численным формулам (5.1), (5.2) или (5.1), (5.10) с контролем или без контроля устойчивости;

• расчеты по формулам (5.1), (5.2) и (5.1), (5.10) с автоматическим выбором численной схемы.

Затем с помощью эксперимента можно выбрать вариант алгоритма, наиболее эффективного на данных задачах.

5.2. АЛГОРИТМ С ПРИМЕНЕНИЕМ СТАДИЙ МЕТОДА РУНГЕ–КУТТЫ–МЕРСОНА Из приведенных выше рассуждений следует, что эффективность алгоритмов интегрирования, основанных на схемах высокого порядка точности, можно повысить за счет применения численных формул бо лее низкого порядка, но с более широкой областью устойчивости.

Расширение области устойчивости и тем самым увеличение величины шага интегрирования не приводят к увеличению вычислительных за трат исходного алгоритма, так как это связано с заменой набора коэф фициентов pi, 1 i m, и с незначительными затратами на оценку максимального собственного числа матрицы Якоби. Отметим, что эф фективность исходного метода не ухудшается, потому что при данном подходе он является составной частью более сложного алгоритма и по признаку может быть применен отдельно. Все это позволяет пересмот реть известные численные схемы с точки зрения построения на их ос нове алгоритмов переменного порядка и шага.

Для изучения возникающих при этом проблем рассмотрим пяти стадийный метод типа Рунге–Кутты, стадии которого совпадают с со ответствующими стадиями метода Мерсона, т. е.

Г л а в а 5. АЛГОРИТМЫ ИНТЕГРИРОВАНИЯ ПЕРЕМЕННОГО ПОРЯДКА И ШАГА yn+1 = yn + pi ki, k1 = hf ( yn ), k2 = hf yn + k1, i = 1 1 k3 = hf yn + k1 + k2, k4 = hf yn + k1 + k3, (5.14) 6 6 1 k5 = hf yn + k1 k3 + 2k4.

2 Метод четвертого порядка При значениях коэффициентов p1 = p5 = 1 / 6, p2 = p3 = 0, p4 = 2 / 3 (5.15) формула (5.14) совпадает с методом Мерсона, а ее многочлен устойчи вости Q5,4 ( z ) имеет вид Q5,4 ( z ) = 1 + z + z 2 / 2 + z 3 / 6 + z 4 / 24 + + z 5 / 144. Обозначим Cn =|| 2k1 9k3 + 8k4 k5 || /150, Vn = 6 max | ( k3 k2 )i | / | (k2 k1 )i |. (5.16) 1i N Тогда для контроля точности схемы (5.14) с коэффициентами (5.15) можно применять неравенство Cn 5/ 4, а для контроля устойчивости следующее:

Vn D, (5.17) где постоянную D можно выбрать равной 3,5, т. е. примерно равной длине интервала устойчивости (5.14), (5.15).

Метод первого порядка Теперь на основе стадий (5.14) построим формулу интегрирования более грубую в смысле точности, но с максимальным интервалом ус тойчивости. Для этого запишем матрицу B5 с элементами (3.21). Учи тывая значения коэффициентов ij, 2 i 5, 1 j i 1, получим 5.2. Алгоритм с применением стадий метода Рунге–Кутты–Мерсона 1 1 1 0 1/ 3 1/ 3 1/ 2 0 1/ 0 1 / 18 1 / 8 (5.18) 0 1 / 48 1 / 0 0 0 0 0 1 / Применим (5.14) для решения задачи y = y, y (0) = y0, Re() 0.

С использованием (3.23), (3.24) и (5.18) получим yn+1 = Q5,1 ( z ) yn, z = h, Q5,1 ( z ) = 1 + pi z + (2 p2 + 2 p3 + 3 p4 + 6 p5 ) z 2 / 6 + i = +(4 p3 + 9 p4 + 36 p5 ) z 3 / 72 + ( p4 + 8 p5 ) z 4 / 48 + p5 z 5 / 24. (5.19) Область устойчивости (5.14) ограничена кривой | Q5,1 ( z ) | = 1 в комплексной плоскости z = h. Выберем коэффициенты pi, 1 i 5, таким образом, чтобы схема (5.14) имела первый порядок точности и максимальный интервал устойчивости. Первое условие будет выпол pi = 1. В результате задача построения численной схемы нено, если i = с максимальным интервалом устойчивости сводится к построению по линома вида Q5,1 ( z ) = 1 + z + ci z i, (5.20) i = для которого неравенство | Q5,1 ( z ) | 1 выполняется на максимальном интервале [, 0]. Приравнивая затем коэффициенты при одинаковых степенях z в (5.19) и (5.20), получим линейную систему уравнений относительно pi, 1 i 5, т. е.

p1 + p2 + p3 + p4 + p5 = 1, p2 / 3 + p3 / 3 + p4 / 2 + p5 = c2, (5.21) p3 / 18 + p4 / 8 + p5 / 2 = c3, p4 / 48 + p5 / 6 = c4, p5 / 24 = c5.

Разрешая (5.21), получим коэффициенты метода первого порядка точ ности с максимальным интервалом устойчивости, т. е.

Г л а в а 5. АЛГОРИТМЫ ИНТЕГРИРОВАНИЯ ПЕРЕМЕННОГО ПОРЯДКА И ШАГА p1 = 1 3c2 + 24c4 48c5, p2 = 3c2 18c3 + 36c4, (5.22) p3 = 18c3 108c4 + 216c5, p4 = 48c4 192c5, p5 = 24c5.

Для получения коэффициентов (5.20) рассмотрим одновременно многочлен Чебышева T5 ( x) пятой степени T5 ( x) = 16 x5 20 x3 + 5 x, для которого выполнено неравенство | T5 ( x) | 1 при x [1,1]. Заме ним переменные, полагая x = 1 2 z /. Получим 50 400 2 1120 3 1280 4 512 T5 ( z ) = 1 z+ z z+ z z, (5.23) 2 3 4 при этом отрезок [, 0] отображается на отрезок [ 1,1]. Потребуем совпадение коэффициентов (5.20) и (5.23) при одинаковых степенях z.

Тогда с применением (5.19) получим, что коэффициенты pi, 1 i 5, схемы (5.14) и коэффициенты многочлена (5.23) связаны соотноше ниями (5.21), где = 50, c2 = 400 / 2, c3 = 1120 / 3, (5.24) 4 c4 = 1280 /, c5 = 512 /.

Подставляя (5.24) в (5.22), получим коэффициенты p1 = 5, 248365568 101, p2 = 3, 260928 101, p3 = 1,395154944 101, p4 = 9,5158272 103, (5.25) p5 = 3,93216 105, при которых схема (5.14) имеет первый порядок точности, а ее интер вал устойчивости расширен до 50 по действительной оси. Функция устойчивости схемы (5.14), (5.25) совпадает с (5.23) при = 50 и имеет вид Q5,1 ( z ) = 1 + z + 0,16 z 2 + 0,00896z 3 + 0,0002048z 4 + 0,0000016384z 5.

Область устойчивости метода первого порядка показана на рис. 5.3.

Теперь построим неравенство для контроля точности и устойчиво сти схемы (5.14) с коэффициентами (5.25). Для этого разложим стадии 5.2. Алгоритм с применением стадий метода Рунге–Кутты–Мерсона ki, 1 i 5, в ряды Тейлора по степеням h и подставим в первую формулу (5.14). С применением (3.32) получим yn+1 = yn + hf n + c2 h 2 f n f n + O (h3 ).

(5.26) Рис. 5.3. Область устойчивости метода (5.14), (5.25) В случае использования коэффициентов (5.25) имеем c2 = 0,16.

Пусть yn = y (tn ). Тогда, сравнивая (3.31) и (5.26) до членов с h включительно, получим, что локальная ошибка n,1 схемы (5.14) пер вого порядка точности имеет вид n,1 = (1 / 2 c2 )h 2 f f + O(h3 ). Фор мулу первого порядка предполагается применять на участке установ ления, где ошибки за счет неточности (5.14) невелики. Поэтому для контроля точности вычислений можно использовать оценку локальной ошибки.

Введем обозначения An = | (3 6c2 ) | || k2 k1 || /2, An = = | (1 2c2 ) | || hf ( yn+1 ) k1 || /2. (5.27) Учитывая, что имеет место k2 k1 = h 2 f n f n / 3 + O (h3 ), (5.28) hf ( yn +1 ) k1 = h 2 f n f n + O (h3 ), для контроля точности вычислений и при выборе величины шага ин тегрирования будем применять соответственно неравенства An, An. (5.29) Г л а в а 5. АЛГОРИТМЫ ИНТЕГРИРОВАНИЯ ПЕРЕМЕННОГО ПОРЯДКА И ШАГА При практической реализации алгоритма во второй формуле (5.27) используется постоянная 0,5 | 3 6c2 | вместо 0,5 | 1 2c2 |. Более жест кое неравенство при выборе величины шага, чем при контроле точно сти вычислений, позволяет избежать некоторых повторных вычисле ний решения (возвратов) вследствие нарушения первого неравенства (5.29). Так как длина интервала устойчивости схемы (5.14), (5.25) рав на 50, то для контроля устойчивости можно применять неравенство (5.17) при D 50.

В силу того, что интервал устойчивости (5.14), (5.25) примерно в 14 раз шире (5.14), (5.15), при переходе с (5.14), (5.15) на (5.14), (5.25) шаг интегрирования, вообще говоря, может быть увеличен в 14 раз. Как показывают расчеты, в ряде случаев такое резкое увели чение шага приводит к тому, что нарушается первое неравенство (5.29) и шаг может быть уменьшен до такой величины, при которой осущест вляется обратный переход на формулы (5.14), (5.15). Такие неоправ данные переходы с одной схемы на другую и обратно приводят к снижению эффективности алгоритма интегрирования.

Чтобы избежать таких переходов, включим в состав алгоритма ин тегрирования схему второго порядка, у которой интервал устойчиво сти шире, чем интервал устойчивости (5.14), (5.15), но же, чем интер вал устойчивости (5.14), (5.25). Переход с четвертого порядка на первый и обратно теперь осуществим посредством формулы второго порядка точности, что позволяет выбирать шаг интегрирования плав но. Заметим, что для схемы третьего порядка точности на основе ста дий ki, 1 i 5, определенных в (5.14), интервал устойчивости рас ширяется незначительно по отношению к (5.14), (5.15), и ее рассматривать не будем.

Метод второго порядка Перейдем к рассмотрению численной формулы вида (5.14) второго порядка точности. Используя (3.32), приближенное решение, опреде ляемое по схеме (5.14), можно записать в виде 5 1 1 yn+1 = yn + pi hf n + p2 + p3 + p4 + p5 h 2 f n f n + 3 3 i = 5.2. Алгоритм с применением стадий метода Рунге–Кутты–Мерсона 1 1 + h3 p2 + p3 + p4 + p5 f n f n2 + 18 18 1 + p3 + p4 + p5 f n 2 f n + O(h 4 ), (5.30) 18 2 где элементарные дифференциалы вычислены в точке yn. Пусть yn = y (tn ). Тогда, сравнивая (3.31) и (5.30), получим условия второго порядка точности схемы (5.14), т. е.

i =1 pi = 1, p2 / 3 + p3 / 3 + p4 / 2 + p5 = 1 / 2. (5.31) Для контроля точности применим неравенство вида (2.30). Сравни вая (3.31) и (5.30) до членов с h3 включительно при требовании yn = y (tn ), видим, что для этого необходимо выполнение соотношения p2 / 9 + p3 / 9 + p4 / 4 + p5 = 1 / 3, (5.32) обеспечивающего согласование локальных ошибок методов первого и второго порядков точности. Прежде чем выписать оценку ошибки схе мы (5.14) при условиях (5.31) и (5.32), рассмотрим функцию устойчи вости Q5,2 ( z ) при полученных равенствах на коэффициенты. Учиты вая в (5.19) соотношения (5.31) и (5.32), получим Q5,2 ( z ) = 1 + z + z 2 / 2 + z 3 / g1 + z 4 / 24 + z 5 / g 2, (5.33) где g1 и g 2 – произвольные постоянные, (1 2 p1 ) g 2 = 96. Коэффици енты pi, 1 i 5, как функции g1 и g 2, определяются из (5.31), (5.32) по формулам p1 = 0,5( g 2 96) / g 2, p2 = 3 18 g1, (5.34) p3 = (72 g1 18 p1 9) / 4, p4 = 4 p1, p5 = (1 2 p1 ) / 4, а локальная ошибка схемы (5.14) имеет вид n,2 = g1 1 ( g1 6)h f 2 f / 6 + O( h 4 ). Тогда для контроля точности применим оценку ошибки вида (2.30), т. е., учитывая (5.27), (5.28), для выбора величины Г л а в а 5. АЛГОРИТМЫ ИНТЕГРИРОВАНИЯ ПЕРЕМЕННОГО ПОРЯДКА И ШАГА шага интегрирования можно использовать соответственно следующие неравенства:

An / d, An / d, d = ( g1 6) / (3 g1 6 g1c2 ), (5.35) где значения An и An вычислены по формулам (5.27), а c2 есть коэф фициент при слагаемом z 2 в полиноме устойчивости метода первого порядка.

Исследуем устойчивость (5.14), (5.34). Как известно, максимальная длина интервала устойчивости формулы второго порядка с пятью вы числениями правой части примерно равна 19. Однако при вычислении ki, 1 i 5, по формулам (5.14) коэффициент при слагаемом z 4 в функции устойчивости (5.33) оказывается фиксированным, равным 1/24. Поэтому интервал устойчивости (5.14), (5.34) ограничен числом 8,6, которое получено на ЭВМ, при этом g1 = 4, 48, g 2 = 393,1. Функ ция устойчивости имеет вид Q5,2 ( z ) = 1 + z + 0,5 z 2 + z 3 / 4,58 + z 4 / 24 + z 5 / 393,1.

Область устойчивости метода (5.14), (5.36) показана на рис. 5.4.

Рис. 5.4. Область устойчивости метода (5.14), (5.36) Подставляя данные значения g1 и g 2 в (5.34), получим p1 = 3,77893665732 101, p2 = 9,30131004367 101, p3 = 2,03904914358 102, p4 = 1,51157466294, (5.36) p5 = 6,1053167133 102.

5.2. Алгоритм с применением стадий метода Рунге–Кутты–Мерсона Если параметр с2 определяется по формуле (5.24), т. е. метод пер вого порядка также используется с максимальным интервалом устой чивости, то в (5.35) имеет место d 0,17. Тогда для контроля устойчи вости схемы (5.14), (5.36) можно применять неравенство (5.17) при D 8,6.

Алгоритм переменного порядка и шага Сформулируем алгоритм переменного порядка и шага на основе построенных формул первого, второго и четвертого порядков точно сти. Дальше будем учитывать, что имеют место соотношения An = O (h 2 ), An = O (h 2 ) и Vn = O(h). Будем предполагать, что Cn (см. формулу (5.16)) при увеличении шага ведет себя как O(h5 ), а при уменьшении как O (h 4 ). Это позволяет избежать некоторых повтор ных вычислений решения вследствие невыполнения точности расче тов. В неравенстве (5.17) будем использовать постоянные D1, D2 и D4, связанные с размером области устойчивости соответственно ме тодов первого, второго и четвертого порядков точности.

Зададим два числа l1 и l2, определяющие количества шагов без увеличения шага интегрирования и без смены численной формулы со ответственно после возврата (повторного вычисления решения вслед ствие невыполнения точности вычислений) и после успешного шага.

Кроме того, определим две переменные l1 и l2, которые при первом обращении к алгоритму положим равными нулю.

Режим работы описанного ниже алгоритма DISPM(IJK) будем оп ределять тремя числами I, J и K. Параметр I может принимать значения 0, 1, 2, 4. При I = 0 расчеты проводятся с переменным по рядком, при I 0 по фиксированной схеме, порядок точности p ко торой совпадает с I, т. е. p = I. При J = 0 устойчивость контролиру ется, при J 0 не контролируется. Параметр K может принимать значения 0, 1 и 2. При K = 0 оценка максимального собственного чис ла n max матрицы Якоби вычисляется с помощью подпрограммы, со ставленной вычислителем;

при K, равном 1, оценка c max вычисля n ется степенным методом;

при K = 2 осредняется по формуле Г л а в а 5. АЛГОРИТМЫ ИНТЕГРИРОВАНИЯ ПЕРЕМЕННОГО ПОРЯДКА И ШАГА n n max = i =1 n max / n, где n – номер текущей точки интегрирования.

cp c В случае J 0 параметр K может быть произвольным. В данных обо значениях рассмотренные выше алгоритмы на основе схемы (5.14) с параметрами (5.15) являются частными случаями алгоритма перемен ного порядка и шага. В частности, алгоритм на основе схемы (3.65) с неравенством для контроля точности вычислений соответствует DISPM(410), а этот же алгоритм с контролем устойчивости с помощью неравенства (4.16) соответствует DISPM (40K), где K может прини мать одно из приведенных выше значений.

Далее будем предполагать, что при первом обращении к алгоритму вычислена стадия k1 в точке (t0, y0 ) по формуле (5.14). Алгоритм ин тегрирования на основе схемы (5.14), (5.25) формулируется так.

Шаг 1а. Вычисляется стадия k2 по формуле (5.14).

Шаг 2а. Вычисляется оценка ошибки An по первой формуле (5.27).

q1 An =.

Шаг 3а. Вычисляется параметр q1 по формуле Шаг 4а. Если q1 1, то l1 полагается равным l1 и происходит по вторное вычисление с шагом q1hn, т. е. происходит переход на шаг 1а.

Шаг 5а. Вычисляются стадии ki, 3 i 5, по формулам (5.14).

Шаг 6а. Вычисляется приближение к решению yn+1 по формулам (5.14), (5.25).

Шаг 7а. Вычисляется стадия kn+11 = hn f ( yn +1 ).

, Шаг 8а. Вычисляется оценка ошибки An по второй формуле (5.27).

q2 An =.

Шаг 9а. Вычисляется параметр q2 по формуле Шаг 10а. Если q2 1, то l1 полагается равным l1 и происходит по вторное вычисление с шагом q2 hn, т. е. происходит переход на шаг 1а.

Шаг 11а. Параметрам l1 и l2 присваиваются значения (l1 1) и (l2 1).

Шаг 12а. Если l1 0 или l2 0, то hn +1 полагается равным hn и выполняется следующий шаг интегрирования по схеме (5.14), (5.25), т. е. происходит переход на шаг 1а.

Шаг 13а. Параметр l2 полагается равным l2.

5.2. Алгоритм с применением стадий метода Рунге–Кутты–Мерсона Шаг 14а. В зависимости от признака K вычисляется оценка Vn максимального собственного числа матрицы Якоби.

Шаг 15а. Вычисляется параметр r по формуле rVn = D1, причем если I 0 и J 0, то r = max( q1, q2 ).

hn +1 = Шаг 16а. Вычисляется новый шаг интегрирования = max ( hn, min(q1, q2, r )hn ).

Шаг 17а. Если I = 1, то выполняется следующий шаг по схеме (5.14), (5.25).

Шаг 18а. Если выполняется неравенство min(q1, q2 )Vn D2, то выполняется следующий шаг интегрирования по схеме первого поряд ка. В противном случае осуществляется переход на формулу второго порядка точности.

Алгоритм интегрирования на основе схемы (5.14), (5.36) формули руется следующим образом.

Шаг 1б. Вычисляется стадия k2 по формуле (5.14).

Шаг 2б. Вычисляется оценка ошибки An по первой формуле (5.27).

Шаг 3б. Вычисляется параметр q1 по формуле q1 An = / d, где значение d определено в (5.35).

Шаг 4б. Если q1 1, то l1 полагается равным l1 и происходит по вторное вычисление решения с шагом q1hn, т. е. осуществляется пере ход на шаг 1б.

Шаг 5б. Вычисляются стадии ki, 3 i 5, по формуле (5.14).

Шаг 6б. Вычисляется решение по формулам (5.14), (5.36).

Шаг 7б. Вычисляется стадия kn+11 = hf ( yn +1 ).

, Шаг 8б. Вычисляется оценка ошибки An по второй формуле (5.27).

q2 An = / d.

Шаг 9б. Вычисляется параметр q2 по формуле Шаг 10б. Если q2 1, то l1 полагается равным l1 и происходит по вторное вычисление с шагом q2 hn, т. е. осуществляется переход на шаг 1б.

l2 присваиваются величины (l1 1) и Шаг 11б. Параметрам l1 и (l2 1).

Г л а в а 5. АЛГОРИТМЫ ИНТЕГРИРОВАНИЯ ПЕРЕМЕННОГО ПОРЯДКА И ШАГА Шаг 12б. Если l1 0 или l2 0, то hn +1 полагается равным hn и выполняется следующий шаг интегрирования по схеме (5.14), (5.36), т. е. осуществляется переход на шаг 1б.

Шаг 13б. Параметр l2 полагается равным l2.

Шаг 14б. Вычисляются оценка ошибки Cn по формуле (5.16) и, в зависимости от параметра K, оценка Vn максимального собственного числа матрицы Якоби.

Шаг 15б. Вычисляется параметр r по формуле rVn = D2, причем если I 0 и J 0, то r = max(q1, q2 ).

Шаг 16б. Вычисляется новый шаг интегрирования hn +1 = = max ( hn, min(q1, q2, r )hn ).

Шаг 17б. Если I = 2, то выполняется следующий шаг по схеме (5.14), (5.36).

Шаг 18б. Вычисляются значения параметров v1 и v2 по формулам v1 = min(q1, q2 ), v2 = min(q1, q2, r ).

Шаг 19б. Если v1Vn D2, v2 An, то управление передается на схему первого порядка. Если v1Vn D4, v1 Cn 5/ 4, то осуществляет ся переход на метод четвертого порядка. В остальных случаях выпол няется следующий шаг интегрирования по формулам (5.14), (5.36).

Алгоритм интегрирования на основе схемы (5.14), (5.15) формули руется следующим образом.

Шаг 1в. Вычисляются стадии ki, 1 i 5, по формулам (5.14).

Шаг 2в. Вычисляется оценка ошибки Cn по первой формуле (5.16).

Шаг 3в. Вычисляется значение параметра q по формуле q 4 Cn = 5 / 4.

Шаг 4в. Если q 1, то l1 полагается равным l1 и происходит по вторное вычисление с шагом qhn, т. е. происходит переход на шаг 1в.

Шаг 5в. Перевычисляется значение параметра q по формуле q 5Cn = 5 / 4.

Шаг 6в. Вычисляется приближение к решению yn +1 по схеме (5.14), (5.15).

5.3. Алгоритм с применением стадий метода Рунге–Кутты–Фельберга пятого порядка Шаг 7в. Вычисляется стадия kn+11 = hf ( yn +1 ).

, Шаг 8в. Параметрам l1 и l2 присваиваются величины (l1 1) и (l2 1).

Шаг 9в. Если l1 0 или l2 0, то hn +1 полагается равным hn и выполняется следующий шаг интегрирования по схеме (5.14), (5.15), т. е. происходит переход на шаг 1в.

Шаг 10в. Параметр l2 полагается равным l2.

Шаг 11в. Вычисляются оценка ошибки An по формуле (5.27) и оценка Vn максимального собственного числа по формуле (5.16).

Шаг 12в. Вычисляется значение параметра r по формуле rVn = D4, причем если I 0 и J 0, то r = q.

Шаг 13в. Вычисляется новый шаг интегрирования hn +1 = = max ( hn, min( r, q )hn ).

Шаг 14в. Если I = 5, то выполняется следующий шаг по схеме (5.14), (5.15).

Шаг 15в. Если qVn D4, q 2 dAn, то управление передается на метод второго порядка. В противном случае выполняется следующий шаг интегрирования по схеме четвертого порядка точности.

Для конкретных расчетов может быть рекомендован следующий набор параметров D1 = 50, D2 = 8,6, D4 = 3,5.

5.3. АЛГОРИТМ С ПРИМЕНЕНИЕМ СТАДИЙ МЕТОДА РУНГЕ–КУТТЫ–ФЕЛЬБЕРГА ПЯТОГО ПОРЯДКА Ниже построено неравенство для контроля устойчивости шести стадийного метода Фельберга пятого порядка точности. На основе стадий этого метода построена численная схема первого порядка с расширенной областью устойчивости. Сформулирован алгоритм ин тегрирования переменного порядка и шага.

Метод Фельберга пятого порядка Рассмотрим задачу Коши для систем обыкновенных дифференци альных уравнений умеренной жесткости Г л а в а 5. АЛГОРИТМЫ ИНТЕГРИРОВАНИЯ ПЕРЕМЕННОГО ПОРЯДКА И ШАГА y = f (t, y ), y (t0 ) = y0, t0 t tk, где y и f вещественные N -мерные вектор-функции, t независи мая переменная. Для простоты ниже часть выкладок проведем для ав тономной системы y = f ( y ). Для решения данной задачи будем при менять метод типа Рунге–Кутты следующего вида:

yn+1 = yn + pi ki, (5.37) i = где k1 = hf (tn, yn ), k2 = hf tn + h, yn + k1, 4 3 3 k3 = hf tn + h, yn + k1 + k2, 8 32 12 1932 k4 = hf tn + h, yn + k1 k2 + k3, (5.38) 13 2197 2197 439 3680 k5 = hf tn + h, yn + k1 8k2 + k3 k4, 216 513 1 8 3544 k6 = hf tn + h, yn k1 + 2k2 k3 + k4 k5.

2 27 2565 4104 При значениях коэффициентов 16 p51 =, p52 = 0, p53 =, 135 (5.39) 28561 9 p54 =, p55 =, p56 = 56430 50 формула (5.37) совпадает с методом Фельберга и имеет пятый порядок точности, а ее функция устойчивости Q6,5 ( z ) имеет вид Q6,5 ( z ) = 1 + z + z 2 / 2 + z 3 / 6 + z 4 / 24 + z 5 / 120 + z 6 / 2080.

5.3. Алгоритм с применением стадий метода Рунге–Кутты–Фельберга пятого порядка Область устойчивости метода (5.37) показана на рис. 5.5.

Рис. 5.5. Область устойчивости метода Фельберга (5.38), (5.39) Известен другой набор коэффициентов:

25 1408 2197 p41 =, p42 = 0, p43 =, p44 =, p45 =, p46 = 0, 216 2565 4104 при которых схема (5.37) имеет четвертый порядок точности. Функция устойчивости Q6,4 ( z ) данного метода имеет вид Q6,4 ( z ) = 1 + z + z 2 / 2 + z 3 / 6 + z 4 / 24 + z 5 / 104, а область устойчивости показана на рис. 5.6.

Рис. 5.6. Область устойчивости метода четвертого порядка Так как в каждой точке имеем два приближения к решению по ме тодам четвертого и пятого порядков, то для контроля точности можно ( p5i p4i )ki n,5 = использовать оценку ошибки вида 24, i = где ki, 1 i 6, определены формулами (5.38);

|| || некоторая норма Г л а в а 5. АЛГОРИТМЫ ИНТЕГРИРОВАНИЯ ПЕРЕМЕННОГО ПОРЯДКА И ШАГА в R N. В результате для контроля точности вычислений можно приме нять неравенство n,5, где требуемая точность расчетов. Из (5.38) видно, что приращение k5 вычисляется в точке tn +1, и поэтому нет необходимости в дополнительном неравенстве для выбора величи ны шага интегрирования.

Теперь перейдем к построению неравенства для контроля устойчи вости. Из анализа результатов расчетов алгоритмами с контролем ус тойчивости следует, что при получении оценки максимального собст венного числа матрицы Якоби можно ограничиться рассуждениями для линейной задачи с постоянными коэффициентами вида y = Ay + b, y (0) = y0. Запишем ki, 1 i 3, применительно к данной задаче. Будем иметь k1 = h( Ayn + b), k2 = k1 + (hA)k1 / 4, k3 = k1 + 3(hA)k1 / 8 + 9(hA) 2 k1 / 128.

k2 k1 = (hA)k1 / 4, (32k3 48k2 + 16k1 ) / 9 = Отсюда получим Vn = max1i N {| (32k3 48k2 + = (hA) k1 / 4. Введем обозначение +16k1 )i | / | (k2 k1 )i |} 9. В результате, с использованием степенного метода оценки максимального собственного числа, для контроля ус тойчивости (5.37) с коэффициентами (5.39) можно применять неравен ство Vn 3,6. В данном неравенстве используется число 3,6, которое примерно равно длине интервала устойчивости схемы (5.37), (5.39).

Отметим, что по мнимой оси область устойчивости также ограни чена числом 3,6. Заметим, так как в случае нелинейной задачи в пред ставлениях ki, 1 i 3, в виде рядов Тейлора имеются дополнитель ные слагаемые, т. е.

k2 k1 = h 2 f n f n / 4 + h3 f n f n2 / 32 + O(h 4 ), (32k3 48k2 + 16k1 ) / 9 = h3 f n 2 f n / 4 + 4h3 f n f n2 / 3 + O(h 4 ), то они будут вносить искажения в оценку максимального собственного числа. Поэтому при выборе шага интегрирования будем учитывать грубость оценки.

5.3. Алгоритм с применением стадий метода Рунге–Кутты–Фельберга пятого порядка Зададимся некоторым числом q 1. Учитывая, что n,5 = O(h5 ) и Vn = O(h), определим два числа q иr по формулам q5 n,5 =, rVn = 3,6. Тогда шаг по точности и устойчивости можно выбрать по формуле hn +1 = min( q, r )hn. Однако в силу грубости оценки максимального собственного числа данное соотношение приводит к неоправданным колебаниям величины шага интегрирования. Результа ты расчетов показывают, что для прогноза величины шага более на дежна следующая формула:


hn +1 = max(hn, min(q, r )hn ). (5.40) Из анализа результатов расчетов жестких задач алгоритмом с кон тролем точности и устойчивости следует значительное повышение эф фективности за счет дополнительного контроля устойчивости. Однако на участках установления точность вычислений значительно лучше задаваемой точности. В такой ситуации выгоднее считать методом бо лее низкого порядка, но с более широкой областью устойчивости.

Метод первого порядка На основе стадий ki, 1 i 6, определяемых формулами (5.38), по строим формулу интегрирования более грубую в смысле точности, но с максимальным интервалом устойчивости. Учитывая значения пара метров ij, 2 i 6, 1 j i 1, запишем матрицу B6 вида (3.21), т. е.

1 1 1 1 1 1 3 12 0 4 8 9 72 1 0 128 169 B6 = 513 5 0 0 2197 0 0 0 0 0 0 0 Г л а в а 5. АЛГОРИТМЫ ИНТЕГРИРОВАНИЯ ПЕРЕМЕННОГО ПОРЯДКА И ШАГА Теперь коэффициенты P6 = ( p1,..., p6 )T метода (5.37) и коэффици енты C6 = (c1,..., c6 )T многочлена устойчивости Q6,1 ( z ) связаны соот ношением типа (3.26), т. е. B6 P6 = C6. Применим (5.37) для решения скалярного тестового уравнения y = y, Re() 0. В результате по лучим yn+1 = Q6,1 ( z ) yn, z = h, 6 1 3 Q6,1 ( z ) = 1 + pi z + p2 + 3 p3 + p4 + p5 + p6 z 2 + 4 8 13 i = 9 1 513 72 1 5 p4 + p5 + p6 z 3 + p6 z 4 + + p3 + p4 + p 128 169 2 8 2197 12 5 11 p6 z 5 + p6 z 6.

+ p5 + (5.41) 104 1248 Область устойчивости (5.37) ограничена кривой | Q6,1 ( z ) | = 1 в комплексной плоскости z = h. Выберем коэффициенты pi, 1 i 6, таким образом, чтобы схема (5.37) имела первый порядок точности и максимальный интервал устойчивости. Первое условие будет выпол pi = 1. В результате задача построения численной схемы нено, если i = с максимальным интервалом устойчивости сводится к построению по линома вида Q6,1 ( x) = 1 + z + ci z i, (5.42) i = для которого неравенство | Q6,1 ( z ) | 1 выполняется на максимальном интервале [, 0]. Приравнивая затем коэффициенты при одинаковых степенях z в (5.41) и (5.42), получим линейную систему алгебраиче ских уравнений относительно pi, 1 i 6. Разрешая данную линей ную систему с верхней треугольной матрицей, получим коэффициенты метода первого порядка точности с максимальным интервалом устой чивости. Здесь они не приведены в явном виде в силу громоздкости соответствующих выражений.

5.3. Алгоритм с применением стадий метода Рунге–Кутты–Фельберга пятого порядка Для получения коэффициентов полинома (5.42) рассмотрим одно временно многочлен Чебышева T6 ( x) шестой степени T6 ( x) = 32 x 6 48 x 4 + 18 x 2 1, (5.43) для которого выполняется неравенство | T6 ( x) | 1 при x [1,1]. Про ведем замену переменных в (5.43), полагая x = 1 2 z /. Получим 72 840 2 3584 T6 ( z ) = 1 z+ z z+ 2 6912 4 6144 5 2048 + z z+ z, (5.44) 4 5 при этом отрезок [, 0] отобразится на отрезок [ 1,1]. Потребуем сов падения коэффициентов многочленов (5.42) и (5.44) при одинаковых степенях z. Тогда получим, что параметры pi, 1 i 6, схемы (5.37) и коэффициенты (5.42) связаны соотношением B6 P6 = C6, где 840 3584 6912 6144 = 72, c2 =, c3 =, c4 =, c5 =, c6 =.

2 3 4 После этого, разрешая систему B6 P6 = C6, получим коэффициенты p1 = +0, 41975960186956;

p2 = 0, 44944365216575;

p3 = +0,1296419611922;

p4 = 0,12199235635231 102 ;

(5.45) p5 = 0,66250690732054 104 ;

p6 = 0,11118997045939 105, при которых схема (5.37) имеет первый порядок точности, а ее интер вал устойчивости расширен до 72 по действительной оси. Функция устойчивости схемы (5.37), (5.45) совпадает с (5.44) при = 72 и имеет вид 35 2 14 Q6,1 ( z ) = 1 + z + z+ z+ 216 14 1 z5 + z6.

+ z+ 3888 314928 Г л а в а 5. АЛГОРИТМЫ ИНТЕГРИРОВАНИЯ ПЕРЕМЕННОГО ПОРЯДКА И ШАГА Область устойчивости данного метода показана на рис. 5.7.

Рис. 5.7. Область устойчивости метода (5.37), (5.45) Построим неравенство для контроля точности схемы (5.37) с коэф фициентами (5.45). Для этого разложим ki, 1 i 6, в ряды Тейлора и подставим в первую формулу (5.37). В результате получим yn+1 = yn + hf n + c2 h 2 f n f n + O(h3 ).

(5.46) В случае (5.45) имеет место c2 = 35 / 216. Пусть yn = y (tn ). Сравнивая (5.46) с разложением точного решения в ряд Тейлора до членов с h включительно, получим, что локальная ошибка n,1 схемы (5.37) пер вого порядка имеет вид n,1 = (1 2c2 ) h 2 f / 2 + O (h3 ).

f Метод первого порядка предполагается использовать на участке установления, где ошибки за счет неточности (5.37) невелики. Поэтому для контроля точности вычислений можно использовать оценку ло кальной ошибки.

Введем обозначения,1 =| 2 4c2 ||| k2 k1 ||,,1 =| 1 2c2 | n n k2 k1 = h 2 f n f n / 4 + O(h3 ) || hf ( yn+1 ) k1 ||. Учитывая, что и hf ( yn +1 ) k1 = h 2 f n f n + O(h3 ), для контроля точности вычислений бу дем применять неравенство,1, а при выборе величины шага до n полнительно будем проверять,1. При реализации алгоритма ин n тегрирования в формуле для вычисления,1 используется постоянная n | 2 4c2 | вместо | 1 2c2 |. Более жесткое неравенство при выборе 5.3. Алгоритм с применением стадий метода Рунге–Кутты–Фельберга пятого порядка шага, чем при контроле точности, позволяет избежать некоторых по вторных вычислений решения (возвратов) вследствие нарушения нера венства,1. Так как длина интервала устойчивости схемы (5.37), n (5.45) равна 72, то для контроля устойчивости можно применять нера венство Vn 72.

Метод второго порядка В силу того, что интервал устойчивости (5.37), (5.45) примерно в 20 раз больше (5.37), (5.39), при переходе с (5.37), (5.39) на (5.37), (5.45) шаг интегрирования, вообще говоря, может быть увеличен в 20 раз. Как показывают расчеты, в ряде случаев такое резкое увеличе ние шага приводит к тому, что нарушается неравенство,1 и шаг n может быть уменьшен до такой величины, что осуществляется обрат ный переход на формулы (5.37), (5.39). Чтобы избежать таких нео правданных переходов, включим в состав алгоритма еще схему второ го порядка, у которой интервал устойчивости шире, чем у схемы (5.37), (5.39), но уже интервала устойчивости (5.37), (5.45). Переход с пятого порядка на первый и обратно будем осуществлять посредством формулы второго порядка точности. Это позволяет выбирать шаг ин тегрирования более плавно. Заметим, что для схемы третьего порядка точности на основе величины ki, 1 i 6, определяемой формулами (5.38), интервал устойчивости расширяется незначительно по отноше нию к (5.37), (5.39), и поэтому ее рассматривать не будем.

Нетрудно видеть, что метод (5.37) имеет второй порядок точности, если его параметры связаны соотношением B6 P6 = C6, где c1 = 1 и c2 = 0,5. В случае набора коэффициентов c1 = 1;

c2 = 0,5;

c3 = 0,8799401907 101;

c4 = 0,6616916777 102 ;

c5 = 0, 2217607053 103 ;

c6 = 0, 2731155893 105, при которых длина интервала устойчивости шестистадийной схемы типа Рунге–Кутты второго порядка точности максимальная и равна примерно 28,5. Рассуждая по аналогии с построением метода первого порядка точности, получим коэффициенты схемы (5.37) второго по рядка, которые имеют вид Г л а в а 5. АЛГОРИТМЫ ИНТЕГРИРОВАНИЯ ПЕРЕМЕННОГО ПОРЯДКА И ШАГА p1 = 0,38402730263584;

p2 = 0, 28983444250334;

p3 = 1,0619634532176;

p4 = 0,36673326913333 101;

(5.47) p5 = 0, 46504946987375 102 ;

p6 = 0, 20657470028730 103.

Рис. 5.8. Область устойчивости метода (5.37), (5.47) Локальная ошибка n,2 схемы (5.37) с коэффициентами (5.47) име ет вид n,2 = (1 6c3 )h3 f 2 f / 6 + O (h 4 ), где c3 = 0,8799401907 102. Тогда для контроля точности вычислений и при выборе величины шага интегрирования можно использовать не равенства,1 / d и,1 / d, где d =| (1 6c3 ) / (3 6c2 ) | ;

n n c2 = 0,16203703703704. При использовании методов первого и второго порядков с максимальными интервалами устойчивости имеем d 0, 24. Для контроля устойчивости схемы (5.37), (5.47) можно ис пользовать неравенство Vn 28,5 (рис. 5.8).

5.4. АЛГОРИТМ С ПРИМЕНЕНИЕМ СТАДИЙ МЕТОДА РУНГЕ–КУТТЫ–ФЕЛЬБЕРГА СЕДЬМОГО ПОРЯДКА Ниже построено неравенство для контроля устойчивости тринадца тистадийного метода Фельберга седьмого порядка точности. На основе первых семи стадий построен метод первого порядка с расширенной областью устойчивости. Сформулирован алгоритм интегрирования переменного порядка и шага.

5.4. Алгоритм с применением стадий метода Рунге–Кутты–Фельберга седьмого порядка Метод Фельберга седьмого порядка Рассмотрим задачу Коши для системы дифференциальных уравне ний y = f (t, y ), y (t0 ) = y0, t0 t tk, (5.48) где y и f – N -мерные вещественные вектор-функции, t – независи мая переменная. Для решения (5.48) будем использовать явные форму лы типа Рунге–Кутты следующего вида:

i yn+1 = yn + pmi ki, ki = hf tn + i h, yn + ij k j, (5.49) i =1 j = где 1 = 0, 2 = 2 / 27, 3 = 1 / 9, 4 = 1 / 6, 5 = 5 / 12, 6 = 1 / 2, 7 = 5 / 6, 8 = 1 / 6, 9 = 2 / 3, 10 = 1 / 3, 11 = 1, 12 = 0, 13 = 1, 21 = 2 / 27, 31 = 1 / 36, 32 = 1 / 12, 41 = 1 / 24, 42 = 0, 43 = 1 / 8, 51 = 5 / 12, 52 = 0, 53 = 25 / 16, 54 = 25 / 16, 61 = 1 / 20, 62 = 63 = 0, 64 = 1 / 4, 65 = 1 / 5, 71 = 25 / 108, 72 = 73 = 0, 74 = 125 / 108, 75 = 65 / 27, 76 = 125 / 54, 81 = 31 / 300, (5.50) 82 = 83 = 84 = 0, 85 = 61 / 225, 86 = 2 / 9, 87 = 13 / 900, 91 = 2, 92 = 93 = 0, 94 = 23 / 108, 95 = 704 / 45, 96 = 107 / 9, 97 = 67 / 90, 98 = 3, 10,1 = 91 / 108, 10,2 = 10,3 = 0, 10,4 = 23 / 108, 10,5 = 976 / 135, 10,6 = 311 / 54, 10,7 = 19 / 60, 10,8 = 17 / 6, 10,9 = 1 / 12, 111 = 2383 / 4100, 11,2 = 11,3 = 0, 11,4 = 341 / 164,, 11,5 = 4496 / 1025, 11,6 = 301 / 82, 11,7 = 2133 / 4100, 11,8 = 45 / 82, 11,9 = 45 / 164, 1110 = 18 / 41, 12,1 = 3 / 205,, Г л а в а 5. АЛГОРИТМЫ ИНТЕГРИРОВАНИЯ ПЕРЕМЕННОГО ПОРЯДКА И ШАГА 12,2 = 12,3 = 12,4 = 12,5 = 0, 12,6 = 6 / 41, 12,7 = 3 / 205, 12,8 = 3 / 41, 12,9 = 3 / 41, 12,10 = 6 / 41, 12,11 = 0, 13,1 = 1777 / 4100, 13,2 = 13,3 = 0, 13,4 = 341 / 164, 13,5 = 4496 / 1025, 13,6 = 289 / 82, 13,7 = 2193 / 4100, 13,8 = 51 / 82, 13,9 = 33 / 164, 13,10 = 12 / 41, 13,11 = 0, 13,12 = 1.

При значениях коэффициентов p71 = 41 / 840, p72 = p73 = p74 = p75 = 0, p76 = 34 / 105, p77 = p78 = 9 / 35, p79 = p7,10 = 9 / 280, (5.51) p7,11 = 41 / 840, p7,12 = p7,13 = схема (5.49), (5.50) имеет седьмой порядок точности.


Контроль точности и устойчивости Известно, что численные формулы (5.49), (5.50) с коэффициентами p81 = p82 = p83 = p84 = p85 = 0, p86 = 34 / 105, p87 = p88 = 9 / 35, (5.52) p89 = p8,10 = 9 / 280, p8,11 = 0, p8,12 = p8,13 = 41 / 840, имеют восьмой порядок точности. Теперь локальную ошибку n ме тода седьмого порядка можно оценить по формуле n = ( p8i p7i ) ki. В результате для контроля точности вычислений i = применяется неравенство || n ||, (5.53) где || || – некоторая норма в R N, – требуемая точность расчетов.

Учитывая, что n = O(h8 ), шаг h ac по точности выбирается по фор муле h ac = qh, где q находится из уравнения q8 || n || =. Если q 1, 5.4. Алгоритм с применением стадий метода Рунге–Кутты–Фельберга седьмого порядка то происходит повторное вычисление решения (возврат) с шагом h, равным qh. В противном случае получаем приближенное решение, а прогнозируемый шаг вычисляется по формуле h ac = qh. Неравенство (5.53) хорошо зарекомендовало себя при решении многих практиче ских задач, здесь мы будем его использовать.

Теперь построим неравенство для контроля устойчивости. Для это го применим (5.49) для решения линейной задачи y = Ay с постоян ной матрицей A. Первые три стадии k1, k2 и k3 схемы (5.49) приме нительно к данной задаче имеют вид k1 = Xyn, k2 = ( X + 2 X 2 / 27) yn, k3 = ( X + X 2 / 9 + X 3 / 162) yn, где X = hA. Нетрудно видеть, что име ют место соотношения 12k3 18k2 + 6k1 = 2 X 3 yn / 27, k2 k1 = 2 X 2 yn / 27.

Теперь оценку максимального собственного числа матрицы Якоби можно вычислить степенным методом. Введем обозначение vn = max {| (12k3 18k2 + 6k1 )i | / | (k2 k1 )i |}. (5.54) 1i N Для контроля устойчивости метода Фельберга можно применять неравенство vn D, где постоянная D ограничивает интервал устой чивости. Устойчивость методов типа Рунге–Кутты исследуется на ска лярном тестовом уравнении y = y, Re() 0. (5.55) Применяя (5.49) с коэффициентами (5.51) для решения (5.55) получим, что функция устойчивости Q13,7 ( x) метода седьмого порядка точности имеет вид Q13,7 ( x) = 1 + c7, i xi, x = h, где i = c7,1 = 1;

c7,2 = 0,5;

c7,3 = 0,16666666666667;

c7,4 = 0, 41666666666667 101;

c7,5 = 0,83333333333333 102 ;

c7,6 = 0,13888888888889 102 ;

c7,7 = 0,19841269841270 103 ;

Г л а в а 5. АЛГОРИТМЫ ИНТЕГРИРОВАНИЯ ПЕРЕМЕННОГО ПОРЯДКА И ШАГА c7,8 = 0, 23165371472663 104 ;

c7,9 = 0, 23671439526314 105 ;

c7,10 = 0,51829448771964 107 ;

c7,11 = 0, 43191207309970 107.

Область устойчивости данного метода показана на рис. 5.9, откуда видно, что интервал устойчивости метода седьмого порядка приблизи тельно равен 5.

Рис. 5.9. Область устойчивости метода (5.49), (5.51) Функция устойчивости Q13,8 ( x) метода восьмого порядка точности (5.49), (5.52) имеет вид Q13,8 ( x) = 1 + c8, i xi, x = h, где i = c8,1 = 1;

c8,2 = 0,5;

c8,3 = 0,16666666666667;

c8,4 = 0, 41666666666667 101;

c8,5 = 0,83333333333333 102 ;

c8,6 = 0,13888888888889 102 ;

c8,7 = 0,19841269841270 103 ;

c8,8 = 0, 24801587301587 104 ;

c8,9 = 0, 23490700935724 105 ;

c8,10 = 0, 23620053064283 106 ;

c8,11 = 0, 25914724385982 107 ;

c8,12 = 0.14397069103323 107.

Область устойчивости метода (5.49), (5.52) показана на рис. 5.10, откуда видно, что интервал устойчивости метода восьмого порядка 5.4. Алгоритм с применением стадий метода Рунге–Кутты–Фельберга седьмого порядка тоже приблизительно равен 5. Поэтому положим в неравенстве для контроля устойчивости D = 5.

Рис. 5.10. Область устойчивости метода (5.49), (5.52) Учитывая, что vn = O (h), шаг h st по устойчивости можно выби рать по формуле h st = rh, где r вычисляется из равенства rvn = D.

Оценка (5.54) является грубой, поэтому контроль устойчивости ис пользуется как ограничитель на размер шага интегрирования. В ре зультате прогнозируемый шаг вычисляется по формуле h : = max h, min( h ac, h st ). (5.56) Данная формула позволяет стабилизировать поведение шага на участ ке установления решения, где определяющую роль играет устойчи вость.

Метод первого порядка На основе стадий численной схемы (5.49) построим метод первого порядка точности с более широкой областью устойчивости. Для этого перепишем (5.49) в виде i m yn+1 = yn + pi ki, yn,i = yn + ij k j, (5.57) i =1 j = 1 i 13, ki = hf (tn + i h, yn,i ), Г л а в а 5. АЛГОРИТМЫ ИНТЕГРИРОВАНИЯ ПЕРЕМЕННОГО ПОРЯДКА И ШАГА где yn,i называют промежуточными, или внутренними, схемами мето да (5.57), yn,1 = yn. Введем в рассмотрение матрицу Bm с элементами bij вида (см. формулу (3.21)) b1i = 1, 1 i m, bki = 0, 2 k m, 1 i k 1, i bki = ij bk 1, j, 2 k m, k i m, j = k где ij, 2 i m, 1 j i 1, определены соотношениями (5.50). При меняя (5.57) для решения линейного уравнения (5.55), получим yn+1 = Qm ( x) yn, x = h, где функция устойчивости Qm ( x) имеет вид m m Qm ( z ) = 1 + ci xi, ci = bij p j. (5.58) i =1 j =i Введем обозначения Cm = (c1,..., cm )T, Pm = ( p1, …, pm )T. Тогда из (5.58) следует, что коэффициенты ci, 1 i m, многочлена устойчиво сти и коэффициенты pi, 1 i m, метода (5.57) связаны соотношением Bm Pm = Cm. (5.59) При заданных ij, 2 i m, 1 j i 1, и ci, 1 i m, коэффициенты pi можно определить из (5.59). Требование первого порядка точности m pi = 1. Поэтому положим c1 = 1. Остальные схемы (5.57) означает i = коэффициенты ci, 2 i m, используем для задания размера и формы области устойчивости.

Известно, что при c1 = 1 для любого m можно так выбрать коэф фициенты многочлена устойчивости, что соответствующая область устойчивости будет расширена до величины 2m 2 по вещественной оси. Однако в этом случае в экстремальных точках xi, 1 i m 1, имеет место | Qm ( xi ) | = 1. Это приводит к тому, что область устойчи вости становится «почти» многосвязной. Поэтому во многих случаях коэффициенты ci, 2 i m, выбираются так, чтобы в экстремальных 5.4. Алгоритм с применением стадий метода Рунге–Кутты–Фельберга седьмого порядка точках выполнялось соотношение Qm ( xi ) = F, | F | 1. За счет подхо дящего выбора F можно влиять на величину области по мнимой оси.

При значении m = 13 максимальный интервал устойчивости равен 338. Подставляя соответствующие ci, 1 i 13, в (5.59), получим ко эффициенты pi, 1 i 13, метода (5.57) первого порядка точности с максимальным интервалом устойчивости. С учетом того, что размер области устойчивости по вещественной оси метода (5.49), (5.51) равен 5, для задач, в которых шаг в основном ограничен по устойчивости, эффективность может быть повышена примерно в 67 раз. Однако из результатов расчетов даже простейших задач вида (5.55) следует, что шаг интегрирования существенно меньше максимально возможного.

Это связано с тем, что области устойчивости внутренних схем (5.57) малы, и при больших шагах получаем большие погрешности в проме жуточных вычислениях.

Из результатов численных экспериментов следует, что наиболее эффективным является метод первого порядка (5.57) с коэффициента ми ij (5.50) при m = 7. Коэффициенты ci, 2 i 7, многочлена ус тойчивости Q7 ( x) с максимальным интервалом устойчивости, равным 98,0, имеют вид c1 = 1;

c2 = 0,16326530612245;

c3 = 0,99958350687216 102 ;

c4 = 0, 29142376293650 103 ;

c5 = 0, 43614440711587 105 ;

c6 = 0,32366931882441 107 ;

c7 = 0,94364232893419 1010, при этом из (5.59) получим p1 = 0,36275196357533;

p2 = +0,19206362739513;

p3 = +0,84515155604427;

p4 = +0,32207503966499;

(5.60) p5 = +0,33295308199673 102 ;

p6 = +0,13204058892796 103 ;

p7 = +0,16906205162854 106.

Область устойчивости данного метода показана на рис. 5.11.

Г л а в а 5. АЛГОРИТМЫ ИНТЕГРИРОВАНИЯ ПЕРЕМЕННОГО ПОРЯДКА И ШАГА Рис. 5.11. Область устойчивости метода (5.57), (5.60) Приведем еще один набор коэффициентов c1 = 1;

c2 = 0,17242757067512;

c3 = 0,11240089243768 101;

c4 = 0,34986371648830 103 ;

c5 = 0,55971678758885 105 ;

c6 = 0, 44431955012931 107 ;

c7 = 0,13862209013567 109, при которых в экстремальных точках выполняется условие | Q7 ( x) | = 0,9.

Подставляя данные значения ci, 2 i 7, в (5.59), получим коэффи циенты p1 = 0,38057403389775;

p2 = +0,96333718256431 101;

p3 = +0,89767190586740;

p4 = +0,38213980246925;

(5.61) p5 = +0, 42473101431975 102 ;

p6 = +0,18104880747501 103 ;

p7 = +0, 24835399856329 метода (5.57) первого порядка точности почти с максимальным интер валом устойчивости, равным приблизительно 90.

Область устойчивости приведена на рис. 5.12. Область устойчиво сти построенного метода первого порядка точности по вещественной оси примерно в 18 раз шире области устойчивости численной схемы (5.49), (5.51). Кроме того, метод первого порядка по числу вычислений правой части задачи (5.48) почти в 2 раза дешевле (5.49), (5.51). По этому для задач, в которых шаг ограничен в основном по устойчиво сти, предполагается теоретическое повышение эффективности в 36 раз.

Построим неравенства для контроля точности и устойчивости (5.57), (5.61). В неравенстве для контроля точности будем применять 5.4. Алгоритм с применением стадий метода Рунге–Кутты–Фельберга седьмого порядка оценку локальной ошибки n,1. С применением разложений точного и приближенного решений в ряды Тейлора по степеням h включительно получим n,1 вида n,1 = (1 2c2 )h 2 ( ft + f y f ) / 2 + O(h3 ), где дифференциалы ft и f y f вычислены на точном решении y (tn ).

Рис. 5.12. Область устойчивости метода (5.57), (5.61) Нетрудно видеть также, что имеет место соотношение k2 k1 = 2h 2 ( f nt + f ny f n ) / 27 + O(h3 ), где дифференциалы f nt и f ny f n вычислены на приближенном реше нии yn. Тогда для контроля точности (5.57), (5.61) можно применять неравенство || k2 k1 || / | d (1 2c2 ) |, (5.62) где d = 27 / 4, || || – некоторая норма в R N, – требуемая точность расчетов.

В построенном неравенстве стадия k1 вычисляется в точке tn, а стадия k2 – в точке (tn + 2h / 27). Так как ни одна стадия не вычисля ется в точке tn+1, то при быстром изменении решения это может при вести к потере точности вычислений. Поэтому в алгоритме интегриро вания контроль неравенства (5.62) используется как предварительный.

С учетом соотношения hf ( yn +1 ) k1 = h 2 ( f nt + f ny f n ) / 2 + O(h3 ) Г л а в а 5. АЛГОРИТМЫ ИНТЕГРИРОВАНИЯ ПЕРЕМЕННОГО ПОРЯДКА И ШАГА окончательное решение по точности принимается проверкой нера венства | 1 2c2 | || hf ( yn +1 ) k1 ||. (5.63) Использование двух неравенств для контроля точности вычисле ний позволяет существенно сократить число повторных вычислений решения вследствие нарушения точности расчетов. Дополнительное сокращение возвратов достигается выбором d = 1 в формуле (5.62).

Далее, так как интервал устойчивости численной схемы (5.57), (5.61) ограничен числом 90, для ее контроля устойчивости можно применять неравенство vn 90, где vn определяется по формуле (5.54).

Алгоритм переменного порядка и шага Методы первого порядка с расширенными областями устойчивости эффективны на участках установления, где шаг ограничен по устойчи вости. Очевидно, что если шаг ограничен по точности, что характерно для переходных участков (погранслоев), то явный метод Эйлера будет в 7 раз эффективнее построенного метода. Естественно, что при высо кой точности расчетов на переходных участках эффективнее будет ме тод (5.49), (5.51). Существенного повышения эффективности можно достигнуть за счет применения каждого метода на том участке, где он наиболее эффективен. В качестве критерия переключения с метода на метод можно использовать неравенство для контроля устойчивости.

При расчетах по методу (5.49), (5.51) переход на численную схему (5.57), (5.61) осуществляется при нарушении неравенства vn 5. При расчетах методом первого порядка обратный переход происходит в случае выполнения vn 5. Вычисления методом первого порядка со провождаются дополнительным (наряду с точностью) контролем нера венства vn 90, а шаг выбирается по формуле вида (5.56).

5.5. АЛГОРИТМ С ПРИМЕНЕНИЕМ СТАДИЙ МЕТОДА ДОРМАНДА–ПРИНСА ВОСЬМОГО ПОРЯДКА Ниже построено неравенство для контроля устойчивости тринадца тистадийного метода Дорманда–Принса восьмого порядка точности.

На основе первых семи стадий построен метод первого порядка с рас 5.5. Алгоритм с применением стадий метода Дорманда–Принса восьмого порядка ширенной областью устойчивости. Сформулирован алгоритм интегри рования переменного порядка и шага.

Метод Дорманда–Принса Рассмотрим задачу Коши для системы дифференциальных уравне ний (5.48). Для решения (5.48) будем использовать явные формулы типа Рунге–Кутты следующего вида:

i yn+1 = yn + pmi ki, ki = hf tn + i h, yn + ij k j, (5.64) i =1 j = где 1 = 0, 2 = 1 / 18, 3 = 1 / 12, 4 = 1 / 8, 93 5 = 5 / 16, 6 = 3 / 8, 7 = 59 / 400, 8 =, 9 =, 200 13 10 =, 11 =, 12 = 13 = 1, 20 42 = 52 = 62 = 63 = 72 = 73 = 82 = 83 = 92 = 93 = 0, 10,2 = 10,3 = 11,2 = 11,3 = 12,2 = 12,3 = 13,2 = 13,3 = 13,12 = 0, 1 1 1 1 3 21 =, 31 =, 32 =, 41 =, 43 =, 51 =, 18 48 16 32 32 75 75 3 3 53 =, 54 =, 61 =, 64 =, 65 =, 64 64 80 16 29443841 71 =, 74 =, (5.65) 614563906 28693883 23124283 75 =, 76 =, 81 =, 1125000000 1800000000 61564180 22789713 84 =, 85 =, 86 =, 158732637 633445777 Г л а в а 5. АЛГОРИТМЫ ИНТЕГРИРОВАНИЯ ПЕРЕМЕННОГО ПОРЯДКА И ШАГА 180193667 39632708 87 =, 91 =, 94 =, 1043307555 573591083 421739975 100302831 95 =, 96 =, 97 =, 2616292301 723423059 800635310 246121993 98 =, 10,1 =, 10,4 =, 3783071287 1340847787 309121744 12992083 10,5 =, 10,6 =, 10,7 =, 1061227803 490766935 393006217 123872331 10,8 =, 10,9 =, 111 =,, 1396673457 1001029789 8478235783 1311729495 11,4 =, 11,5 =, 11,6 =, 508512852 1432422823 48777925059 15336726248 11,7 =, 11,8 =, 11,9 =, 3047939560 1032824649 3065993473 185892177 1110 =, 12,1 =, 12,4 =,, 597172653 718116043 477755414 703635378 12,5 =, 12,6 =, 12,7 =, 1098053517 230739211 5232866602 4093664535 12,8 =, 12,9 =, 12,10 =, 850066563 808688257 65686358 403863854 12,11 =, 13,1 =, 13,4 =, 487910083 491063109 411421997 652783627 13,5 =, 13,6 =, 13,7 =, 543043805 914296604 13158990841 3936647629 13,8 =, 13,9 =, 13,10 =, 6184727034 1978049680 13,11 =.

5.5. Алгоритм с применением стадий метода Дорманда–Принса восьмого порядка При значениях коэффициентов 13451932 p81 =, p82 = p83 = p84 = p85 = 0, p86 =, 455176623 1757004468 656045339 p87 =, p88 =, p89 =, (5.66) 5645159321 265891186 465885868 53011238 p8,10 =, p8,11 =, p8,12 =, p8,13 = 322736535 667516719 схема (5.64), (5.65) имеет восьмой порядок точности.

Контроль точности и устойчивости Известно, что численные формулы (5.64), (5.65) с коэффициентами 14005451 p71 =, p72 = p73 = p74 = p75 = 0, p76 =, 335480064 181606767 561292985 p77 =, p78 =, p79 =, (5.67) 758867731 797845732 760417239 p7,10 =, p7,11 =, 1151165299 528747749 p7,12 =, p7,13 = 2220607170 имеют седьмой порядок. Тогда для контроля точности схемы восьмого n порядка можно использовать оценку ошибки вида n = i =1 ( p8i p7i )ki. В результате для контроля точности вычисле ний применяется следующее неравенство:

|| n ||, (5.68) где || || – некоторая норма в R N, – требуемая точность расчетов.

Учитывая соотношение n = O(h8 ), шаг h ac по точности выбирается по формуле h ac = qh, (5.69) Г л а в а 5. АЛГОРИТМЫ ИНТЕГРИРОВАНИЯ ПЕРЕМЕННОГО ПОРЯДКА И ШАГА где значение q находится из уравнения q8 || n || =. Если q 1, то происходит повторное вычисление решения (возврат) с шагом h, рав ным qh. В противном случае вычисляется приближенное решение, а прогнозируемый шаг вычисляется по формуле (5.69). Неравенство (5.68) хорошо зарекомендовало себя при решении многих практиче ских задач и будет использоваться ниже.

Построим неравенство для контроля устойчивости. Для этого при меним (5.64) для решения линейной задачи y = Ay, y (0) = y0 с посто янной матрицей A. Первые три стадии k1, k2 и k3 схемы (5.64) при менительно к данной задаче имеют вид k1 = Xyn, k2 = ( X + X 2 / 18) yn, k3 = ( X + X 2 / 12 + X 3 / 288) yn, где X = hA. Нетрудно видеть, что име ет место 8(2k3 3k2 + k1 ) = X 3 yn / 18 и k2 k1 = X 2 yn / 18. Теперь оценку максимального собственного числа матрицы Якоби можно вы числить степенным методом. Введем обозначение vn = 8 max {| (2k3 3k2 + k1 )i | / | ( k2 k1 )i |}. (5.70) 1i N Тогда для контроля устойчивости метода Дорманда–Принса можно применять неравенство vn D, (5.71) где постоянная D ограничивает интервал устойчивости.

Применяя (5.64) с коэффициентами (5.66) для решения (5.55), по лучим, что функция устойчивости Q8 ( x) метода восьмого порядка точности имеет вид Q8 ( x) = 1 + c8, i xi, x = h, где i = c8,1 = 1;

c8,2 = 0,5;

c8,3 = 0,16666666666667;

c8,4 = 0, 41666666666667 101;

c8,5 = 0,83333333333333 102 ;

c8,6 = 0,13888888888889 102 ;

c8,7 = 0,19841269841270 103 ;

c8,8 = 0, 24801587301587 104 ;

5.5. Алгоритм с применением стадий метода Дорманда–Принса восьмого порядка c8,9 = 0, 27521279901047 105 ;

c8,10 = 0, 24231996586959 106 ;

c8,11 = 0, 24389718205443 107 ;

c8,12 = 0, 20346152896858 109.

Область устойчивости метода восьмого порядка показана на рис. 5.13.

Рис. 5.13. Область устойчивости метода (5.64), (5.66) восьмого порядка точности Из рис. 5.9 видно, что интервал устойчивости метода седьмого по рядка приблизительно равен 5. Функция устойчивости Q7 ( x) метода седьмого порядка точности имеет вид Q7 ( x) = 1 + c7, i xi, x = h, где i = c7,1 = 1;

c7,2 = 0,5;

c7,3 = 0,16666666666667;

c7,4 = 0, 41666666666667 101;

c7,5 = 0,83333333333333 102 ;

c7,6 = 0,13888888888889 102 ;

c7,7 = 0,19841269841270 103 ;

c7,8 = 0, 25044253893357 104 ;

c7,9 = 0, 25810595673311 105 ;

c7,10 = 0, 27974369233314 106 ;

c7,11 = 0,10950482491122 107 ;

c7,12 = 0,10214473550959 10 9.

Область устойчивости метода седьмого порядка показана на рис. 5.14.

Интервал устойчивости метода седьмого порядка тоже приблизитель но равен 5. Поэтому положим в неравенстве (5.71) D = 5.

Г л а в а 5. АЛГОРИТМЫ ИНТЕГРИРОВАНИЯ ПЕРЕМЕННОГО ПОРЯДКА И ШАГА Рис. 5.14. Область устойчивости метода (5.64), (5.67) седьмого порядка точности Учитывая, что vn = O (h), шаг h st по устойчивости можно выби рать по формуле h st = rh, где r вычисляется из равенства rvn = D.

Оценка (5.70) является грубой, поэтому здесь контроль устойчиво сти используется как ограничитель на размер шага интегрирования. В результате прогнозируемый шаг вычисляется по формуле (5.56). Отме тим, что формула (5.56) применяется для прогноза величины шага ин тегрирования hn +1 после успешного вычисления решения c предыду щим шагом hn и поэтому фактически не приводит к увеличению вычислительных затрат.

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

Отметим, что в зависимости от интервала интегрирования и размерно сти системы дифференциальных уравнений одна и та же задача может быть отнесена к нежесткой, умеренной жесткости или жесткой. Если интервал интегрирования таков, что решение представляет собой толь ко переходный участок, соответствующий максимальному собствен ному числу матрицы Якоби, а участок установления отсутствует, то такая задача к жестким не относится. Ее решение может быть вычис лено явным методом с небольшой областью устойчивости. Однако ес ли участок установления велик, то могут возникнуть проблемы с при менением явных методов, даже с контролем устойчивости.



Pages:     | 1 | 2 || 4 | 5 |   ...   | 9 |
 





 
© 2013 www.libed.ru - «Бесплатная библиотека научно-практических конференций»

Материалы этого сайта размещены для ознакомления, все права принадлежат их авторам.
Если Вы не согласны с тем, что Ваш материал размещён на этом сайте, пожалуйста, напишите нам, мы в течении 1-2 рабочих дней удалим его.