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

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

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


Pages:     | 1 |   ...   | 2 | 3 || 5 | 6 |   ...   | 9 |

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

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

5.5. Алгоритм с применением стадий метода Дорманда–Принса восьмого порядка Метод первого порядка На основе стадий численной схемы (5.64) построим метод первого порядка точности с более широкой областью устойчивости. Для этого перепишем (5.64) в виде i m yn+1 = yn + pi ki, yn,i = yn + ij k j, (5.72) i =1 j = ki = hf (tn + i h, yn,i ), где yn,i называют промежуточными или внутренними схемами метода (5.72).

Применяя (5.72) для решения линейного скалярного уравнения (5.55), получим yn+1 = Qm ( x) yn, x = h, где функция устойчивости Qm ( x) имеет вид m m Qm ( z ) = 1 + ci xi, ci = bij p j. (5.73) i =1 j =i Обозначим Cm = (c1,..., cm )T и Pm = ( p1, …, pm )T. Тогда из (5.73) следует, что коэффициенты ci, 1 i m, многочлена устойчивости и коэффициенты pi, 1 i m, метода (5.72) связаны соотношением Bm Pm = Cm, (5.74) где элементы bij, 1 i, j m, матрицы Bm определены выше. Теперь при заданных ij, 2 i m, 1 j i 1, и ci, 1 i m, коэффициенты pi, 1 i m, можно определить из (5.74). Требование первого поряд m pi = 1. Поэтому положим c1 = 1.

ка точности схемы (5.72) означает i = Остальные коэффициенты ci, 2 i m, используем для задания раз мера и формы области устойчивости. Далее будем считать их задан ными.

При m = 13 максимальный интервал устойчивости равен 338. Под ставляя соответствующие ci, 1 i 13, в (5.74), получим коэффициен Г л а в а 5. АЛГОРИТМЫ ИНТЕГРИРОВАНИЯ ПЕРЕМЕННОГО ПОРЯДКА И ШАГА ты pi, 1 i 13, метода (5.72) первого порядка точности с максималь ным интервалом устойчивости. С учетом того, что размер области ус тойчивости по вещественной оси метода (5.64), (5.66) равен 5, для за дач, в которых шаг в основном ограничен по устойчивости, эффективность может быть повышена примерно в 67 раз. Однако из результатов расчетов даже простейших задач, как и в тринадцатиста дийном методе Фельберга, следует, что шаг интегрирования сущест венно меньше максимально возможного. Это связано с тем, что облас ти устойчивости внутренних схем (5.72) малы, и при больших шагах получаем большие погрешности в промежуточных вычислениях.

Из результатов численных экспериментов следует, что наиболее эффективным является метод первого порядка (5.72) с коэффициента ми ij (5.65) при 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.74) получим p1 = 0, 43635713190292;

p2 = 0,39757930691747;

p3 = 1,1027283617527;

p4 = +0,72030701125358;

p5 = +0,10208963607634 101;

p6 = +0,56373316433595 103 ;

p7 = +0,12836904213518 103.

Приведем еще один набор коэффициентов:

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, 5.5. Алгоритм с применением стадий метода Дорманда–Принса восьмого порядка при которых в экстремальных точках xi, 1 i 6, многочлена устой чивости выполняется условие | Q7 ( x) | = 0,9. Подставляя данные зна чения ci, 2 i 7, в (5.74), получим коэффициенты p1 = 0, 41342955189830;

p2 = 0,57548324135785;

p3 = +1,1243725642680;

p4 = +0,85058623738482;

p5 = +0,12991731772814 101;

(5.75) p6 = +0,77368430693719 103 ;

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

Область устойчивости построенного метода первого порядка точ ности по вещественной оси примерно в 18 раз шире области устойчи вости численной схемы (5.64), (5.66). Кроме того, метод первого по рядка по числу вычислений правой части задачи (5.48) почти в 2 раза дешевле данной численной схемы. Поэтому для задач, в которых шаг ограничен в основном по устойчивости, предполагаемое теоретическое повышение эффективности в 36 раз, как и в тринадцатистадийном ме тоде Фельберга. Заметим, что коэффициенты методов первого порядка и коэффициенты многочленов устойчивости в алгоритмах на основе стадий численных формул Рунге–Кутты–Фельберга и Дорманда– Принса совпадают, что вполне естественно.

Построим неравенства для контроля точности и устойчивости ме тода (5.72). В неравенстве для контроля точности будем применять оценку локальной ошибки n,1. С применением разложений точного и приближенного решений в ряды Тейлора по степеням h до членов с h2 n,1 n,1 = включительно получим, что имеет вид = (1 2c2 )h 2 ( ft + f y f ) / 2 + O(h3 ), где элементарные дифференциалы ft и f y f вычислены на точном решении y (tn ). Нетрудно видеть так же, что имеет место соотношение k2 k1 = h 2 ( f nt + f ny f n ) / 18 + O (h3 ), где элементарные дифференциалы f nt и f ny f n вычислены на прибли Г л а в а 5. АЛГОРИТМЫ ИНТЕГРИРОВАНИЯ ПЕРЕМЕННОГО ПОРЯДКА И ШАГА женном решении yn. Тогда для контроля точности (5.72) можно при менять неравенство || k2 k1 || / | d (1 2c2 ) |, (5.76) где d = 9, || || – некоторая норма в R N, – точность расчетов. В по строенном неравенстве стадия k1 вычисляется в точке tn, а стадия k2 – в точке (tn + h / 18). Так как ни одна стадия не вычисляется в точке tn+1, то при быстром изменении решения это может приводить к поте ре точности вычислений. Поэтому в алгоритме интегрирования кон троль (5.76) используется как предварительный. С учетом соотноше ния hf ( yn +1 ) k1 = h 2 ( f nt + f ny f n ) / 2 + O(h3 ) окончательное решение по точности принимается проверкой следующего неравенства:

|| hf ( yn +1 ) k1 || / | 1 2c2 |.

Использование двух неравенств для контроля точности вычисле ний позволяет существенно сократить число повторных вычислений решения вследствие нарушения точности расчетов. Дополнительное сокращение возвратов достигается выбором d = 1 в формуле (5.76).

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

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

При расчетах по методу (5.64), (5.66) переход на численную схему (5.72), (5.75) осуществляется при нарушении неравенства vn 5. При 5.5. Алгоритм с применением стадий метода Дорманда–Принса восьмого порядка расчетах методом первого порядка обратный переход происходит в случае выполнения vn 5. Вычисления методом первого порядка со провождаются дополнительным (наряду с точностью) контролем нера венства vn 90, а шаг выбирается по формуле типа (5.56).

Метод Дорманда–Принса получил широкое распространение. Его программная реализация входит, пожалуй, во все известные библиоте ки программ. Численное исследование данного алгоритма проводилось многими авторами. Нам не хотелось бы лишний раз подтверждать преимущества данного метода сравнением с другими алгоритмами с автоматическим выбором шага интегрирования при высокоточных расчетах, ибо это слишком известный факт. В частности, в монографии [177] отмечено, что метод Дорманда–Принса дает «прекрасные чис ленные результаты...», если он применяется для решения нежестких задач. Поэтому здесь ограничились применением данного метода для решения жестких задач.

Тестовые расчеты показали повышение эффективности за счет до полнительного контроля устойчивости, а также за счет расчетов с пе ременным порядком. Наиболее показательна задача моделирования реакции Белоусова–Жаботинского (орегонатор), при решении которой алгоритмом без контроля устойчивости происходит почти 800 000 по вторных вычислений решения за счет неустойчивости численной схе мы. Это означает увеличение вычислительных затрат примерно на 9 000 000 дополнительных вычислений правой части исходной задачи, ибо повторное вычисление решения сопровождается двенадцатью до полнительными вычислениями функции f – хотя метод тринадцати стадийный, стадия k1 зависит от размера шага линейно, и поэтому в случае возврата ее вычислять не нужно. Контроль устойчивости по зволяет существенно сократить неоправданные возвраты и тем самым повысить эффективность алгоритма интегрирования. По сути дела, контроль устойчивости расширяет возможности метода Дорманда– Принса применительно к задачам умеренной жесткости.

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

Применение на участке установления методов низкого порядка точности с расширенными областями устойчивости позволяет значи тельно увеличить размер шага интегрирования без увеличения вычис лительных затрат. На переходных участках, где определяющую роль играет точность вычислений, эффективны методы более высокого по рядка точности, но с небольшой областью устойчивости. Комбиниро вание методов низкого и высокого порядков с помощью неравенства для контроля устойчивости позволяет значительно повысить эффек тивность расчетов.

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

5.6. АЛГОРИТМ НА ОСНОВЕ ДВУХСТАДИЙНОЙ СХЕМЫ В приведенных выше алгоритмах интегрирования при построении неравенства для контроля устойчивости требовалось не менее трех стадий. Целью данного раздела является построение неравенства для контроля устойчивости двухстадийного метода. Техника построения неравенства для контроля устойчивости несколько отличается от при мененной выше.

5.6. Алгоритм на основе двухстадийной схемы Метод второго порядка точности Для численного решения задачи Коши для системы обыкновенных дифференциальных уравнений (2.2) рассмотрим явный двухстадийный метод типа Рунге–Кутты вида yn+1 = yn + p1k1 + p2 k2, k1 = hf ( yn ), k2 = hf ( yn + k1 ). (5.77) В случае неавтономной системы y = f (t, y ) схема (5.79) записыва ется в виде yn +1 = yn + p1k1 + p2 k2, k1 = hf (tn, yn ), k2 = hf (tn + h, yn + k1 ).

Ниже для сокращения выкладок будем рассматривать задачу (2.2). Од нако построенные методы можно применять для решения неавтоном ных задач.

Получим соотношения на коэффициенты метода (5.77) второго по рядка точности. Для этого разложим стадии k1 и k2 в ряды Тейлора по степеням h до членов с h3 включительно и подставим в первую фор мулу (5.77). В результате получим yn+1 = yn + ( p1 + p2 )hf n + p2 h 2 f n f n + +2 p2 h3 f n f n2 / 2 + O (h 4 ), (5.78) где элементарные дифференциалы вычислены на приближенном ре шении yn. Точное решение y (tn +1 ) в окрестности точки tn имеет вид y (tn +1 ) = y (tn ) + hf + h 2 f / 2 + h3 ( f 2 + f 2 f ) / 6 + O(h 4 ), f f где элементарные дифференциалы вычислены на точном решении y (tn ). Пусть yn = y (tn ). Сравнивая полученные ряды для точного и приближенного (5.78) решений до членов с h 2, запишем условия вто рого порядка точности схемы (5.77), которые имеют вид p1 + p2 = 1 и p2 = 1 / 2. При данных соотношениях локальная ошибка имеет вид n = (2 3)h3 f 2 / 12 + h3 f 2 f / 6 + O(h 4 ). Построим неравенство для f контроля точности вычислений метода второго порядка. Для этого рассмотрим вспомогательную схему yn+11 = yn + k1 первого порядка, Г л а в а 5. АЛГОРИТМЫ ИНТЕГРИРОВАНИЯ ПЕРЕМЕННОГО ПОРЯДКА И ШАГА точности. С помощью идеи вложенных методов ошибку n,2 метода второго порядка можно оценить по формуле n,2 = yn +1 yn +11 =, = p2 (k2 k1 ). Для повышения надежности данной оценки выберем = 1. Тогда стадия k1 вычисляется в точке tn, а k2 в точке tn +1.

Как показывают расчеты, использование производных решения в крайних точках шага приводит к более надежному неравенству для контроля точности вычислений. При = 1 коэффициенты метода вто рого порядка определяются однозначно = 1 и p1 = p2 = 1 / 2, а локальная ошибка n и неравенство для контроля точности вычисле ний имеют соответственно вид n = h3 f 2 f / 6 h3 f / 12 + O(h 4 ) и f || k2 k1 || /2, где || || – некоторая норма в R N, – точность интег рирования.

Теперь построим неравенство для контроля устойчивости числен ной формулы (5.77). Для этого рассмотрим вспомогательную стадию k3 = hf ( yn+1 ). Заметим, что k3 совпадает со стадией k1, которая при меняется на следующем шаге интегрирования, и поэтому ее использо вание не приводит к дополнительным вычислениям правой части систе мы (2.2). Запишем стадии k1, k2 и k3 применительно к задаче y = Ay, y (0) = y0, где A есть матрица с постоянными коэффициентами. В ре зультате получим: k1 = Xyn, k2 = ( X + X 2 ) yn, k3 = ( X + X 2 + X 3 / 2) yn, X = hA. Легко видеть, что имеют место равенства где k2 k1 = X 2 yn, 2(k3 k2 ) = X 3 yn. Тогда оценку максимального собст венного числа vn,2 = h n max матрицы Якоби системы (2.2) можно вы числить по формуле vn,2 = 2 max {| (k3 k2 )i | / | (k2 k1 )i |}. Интервал 1i N устойчивости численной схемы (5.77) приблизительно равен 2. Поэто му для ее контроля устойчивости можно применять неравенство vn,2 2.

Полученная оценка является грубой, поэтому ниже контроль ус тойчивости используется как ограничитель на размер шага интегриро вания. В результате прогнозируемый шаг hn +1 вычисляется следую щим образом. Новый шаг h ac по точности определим по формуле 5.6. Алгоритм на основе двухстадийной схемы h ac = q1hn, где hn есть последний успешный шаг интегрирования, а q1, учитывая соотношение k2 k1 = O(hn ), задается уравнением h st q1 || k2 k1 || =. Шаг по устойчивости зададим формулой st h = q2 hn, где q2 с учетом равенства vn,2 = O (hn ) определяется соот ношением q2 vh,2 = 2. В результате прогнозируемый шаг hn +1 вычис ляется по следующей формуле: hn +1 = max[hn, min(h ac, h st )]. Данная формула позволяет стабилизировать поведение шага на участке уста новления решения, где определяющую роль играет устойчивость.

Метод первого порядка точности Для численного решения задачи (2.2) рассмотрим схему вида yn+1 = yn + r1k1 + r2 k2, k1 = hf ( yn ), k2 = hf ( yn + k1 ). (5.79) Заметим, что при r1 = r2 = 1 / 2 численная формула (5.79) имеет второй порядок точности. Построим менее точную схему с максимальным ин тервалом устойчивости. Для этого применим (5.79) для решения ска лярного тестового уравнения y = y, y (0) = y0, Re( ) 0. Получим yn +1 = Q( x) yn, где функция устойчивости Q ( x) имеет вид Q( z ) = = 1 + (r1 + r2 ) z + r2 z 2, z = h. Требование первого порядка точности приводит к соотношению r1 + r2 = 1, которое ниже будем считать вы полненным. Теперь выберем r2 таким образом, чтобы метод (5.79) имел максимальный интервал устойчивости. Для этого рассмотрим многочлен Чебышева T2 ( x) = (2 x 2 1) на промежутке [1,1]. Проведем замену переменных, полагая x = 1 2 z /. Получим T2 ( z ) = 1 8 z / + +8 z 2 / 2, при этом отрезок [, 0] отображается на отрезок [ 1,1]. Не трудно показать, что среди всех многочленов вида P2 ( z ) = 1 + z + c2 z для T2 ( z ) неравенство | T2 ( z ) | 1 выполняется на максимальном ин тервале [, 0], = 8. Потребуем совпадения коэффициентов Q( z ) и T2 ( z ) при = 8. Это приводит к соотношениям r1 + r2 = 1, r2 = 1 / 8.

В результате имеем коэффициенты r1 = 7 / 8 и r2 = 1 / 8 метода порядка Г л а в а 5. АЛГОРИТМЫ ИНТЕГРИРОВАНИЯ ПЕРЕМЕННОГО ПОРЯДКА И ШАГА точности с максимальным интервалом устойчивости, локальная ошиб ка n которого имеет вид n = 3h 2 f / 8 + O (h3 ). Для контроля точно f сти численной формулы первого порядка используем оценку локаль ной ошибки. С учетом того, что k2 k1 = h 2 f n f n + O (h3 ), а также вида локальной ошибки получим неравенство для контроля точности сле дующего вида: || k2 k1 || 8 / 3, где || || – некоторая норма в R N, – требуемая точность расчетов.

Функция устойчивости Q2,1 ( z ) метода первого порядка имеет вид Q2,1 ( z ) = 1 + z + z 2 / 8, а область устойчивости показана на рис. 5.15.

Рис. 5.15. Область устойчивости метода первого порядка точности Построим неравенство для контроля устойчивости метода первого порядка. Для этого снова рассмотрим вспомогательную стадию k3 = hf ( yn+1 ). Запишем стадии k1, k2 и k3 применительно к задаче y = Ay, где A есть матрица с постоянными коэффициентами. В ре зультате получим k1 = Xyn, k2 = ( X + X 2 ) yn, k3 = ( X + X 2 + X 3 / 8) yn, где X = Ay. Легко видеть, что k2 k1 = X 2 yn, 8(k3 k2 ) = X 3 yn. Тогда оценку максимального собственного числа vn,1 = h n max матрицы Якоби системы (1) можно вычислить по формуле vn,1 = 8 max {| (k3 k2 )i | / | ( k2 k1 )i |}.

1i N Интервал устойчивости численной схемы (5.79) равен 8. Поэтому для ее контроля устойчивости можно применять неравенство vn,1 8.

5.7. Многочлены устойчивости Алгоритм интегрирования переменного порядка и шага На основе построенных методов первого и второго порядков точ ности легко сформулировать алгоритм переменного порядка и шага.

Расчеты всегда начинаются методом второго порядка как более точ ным. Переход на схему первого порядка осуществляется при наруше нии неравенства vn,2 2. Обратный переход на метод второго порядка происходит в случае выполнения неравенства vn,1 2. При расчетах по методу первого порядка наряду с точностью контролируется устойчи вость, а выбор прогнозируемого шага производится по аналогии с ме тодом второго порядка точности.

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

Основные соотношения Для упрощения выкладок рассмотрим задачу Коши для автоном ной системы обыкновенных дифференциальных уравнений y = f ( y ), y (t0 ) = y0, t0 t tk, Г л а в а 5. АЛГОРИТМЫ ИНТЕГРИРОВАНИЯ ПЕРЕМЕННОГО ПОРЯДКА И ШАГА для решения которой применим явные методы типа вида i m i +1, j k j, 1 i m 1, yn +1 = yn + pmi ki, yn,i = yn + (5.80) j =1 i = ki = hf ( yn,i 1 ), 1 i m, yn,0 = yn.

Все полученные ниже результаты можно использовать для неавто номных задач. Потребуется матрица Bm с элементами bij вида (см.

главу 3) b1i = 1, 1 i m, bki = 0, 2 k m, 1 i k 1, i 1 (5.81) bki = ij bk 1, j, 2 k m, k i m, j = k где ij – коэффициенты схемы (5.80).

Устойчивость одношаговых методов исследуется на линейном уравнении y = y, y (0) = y0, t 0, (5.82) с комплексным. Применяя вторую формулу (5.80) к (5.82), получим m yn+1 = Qm ( z ) yn, z = h, Qm ( z ) = 1 + cmi z i, i = (5.83) m bij pmj, 1 i m.

cmi = j =i Обозначая Cm = (cm1,..., cmm )T и Pm = ( pm1, …, pmm )T, последнее со отношение (5.83) можно переписать в виде Bm Pm = Cm, (5.84) где элементы матрицы Bm определены соотношениями (5.81). Для промежуточных численных схем (5.80) имеем k yn,k = Qk ( z ) yn, Qk ( z ) = 1 + cki z i, i = (5.85) k bij k +1, j, 1 k m 1.

cki = j =i 5.7. Многочлены устойчивости Используя обозначения k = (k +11,…,k +1,k )T и ck = (ck1,…, ckk )T,, получим, что коэффициенты ij промежуточных схем (5.80) и коэф фициенты соответствующих многочленов устойчивости связаны соот ношениями Bk k = ck, 1 k m 1. (5.86) Из сравнения (5.81) и (5.85) следует, что bki = ci 1,k 1, т. е. элементы (k + 1) -го столбца матрицы Bm совпадают с коэффициентами много члена устойчивости Qk ( z ). Отсюда следует, что если заданы коэффи циенты многочленов устойчивости основной и промежуточной чис ленных схем, то коэффициенты методов (5.80) однозначно определяются из линейных систем (5.84) и (5.86) с верхними треуголь ными матрицами Bi, 1 i m.

Разлагая точное и приближенное решения в ряды Тейлора, запишем y (tn+1 ) = y (tn ) + hf + h 2 f / 2 + h3 f 2 + f 2 f 6 + O (h 4 ), f f m m yn+1 = yn + b1 j pmj hf n + b2 j pmj h 2 f n f n + j =1 j =2 m m + b3 j pmj h3 f n 2 f n + 0,5 b2 j pmj h3 f n f n2 + O(h 4 ), (5.87) j =3 j =2 где элементарные дифференциалы вычислены на точном y (tn ) и при ближенном yn решениях соответственно. Из сравнения соотношений (5.87) при условии yn = y (tn ) видно, что численная формула (5.80) бу m b1 j pmj = 1.

дет иметь первый порядок точности, если Требование j = второго порядка точности (5.80) означает выполнение дополнительно m b2 j pmj = 1 / 2. Наконец схема (5.80) будет иметь третий го условия j = порядок, если дополнительно выполняются еще два соотношения:

m 1 m2 b3 j pmj =, b2 j pmj =. (5.88) 6 j =2 j = Г л а в а 5. АЛГОРИТМЫ ИНТЕГРИРОВАНИЯ ПЕРЕМЕННОГО ПОРЯДКА И ШАГА Отсюда следует, что для построения m -стадийных методов перво го и второго порядков точности в линейной системе (5.84) следует по ложить cm1 = 1 и cm1 = 1, cm 2 = 0,5 соответственно. Задача о построе нии m -стадийных методов третьего порядка сводится к совместному решению линейной системы (5.84) и второго уравнения (5.88) при ус cm1 = 1, cm 2 = 1 / 2, cm3 = 1 / 6. Учитывая, что b2 j = j, ловии 2 j m, нетрудно видеть, что совместность (5.84) и (5.88) эквива лентна совместности соотношений m m 1 j pmj = 2, 2j pmj = 3, j =2 j = где величина j, 2 j m, определена формулами (5.81).

Конструирование областей устойчивости Для того чтобы воспользоваться соотношениями (5.84) и (5.86), требуются коэффициенты многочленов устойчивости. Функция устой чивости явного m -стадийного метода типа Рунге–Кутты представляет собой полином степени m. Через заданные коэффициенты схемы (5.80) можно вычислить коэффициенты этого многочлена по формуле (5.84). Здесь рассмотрим задачу получения таких коэффициентов мно гочлена устойчивости, чтобы область устойчивости схемы имела за данную, естественно «разумную», форму и размер.

Итак, пусть заданы два числа k и m, k m. Рассмотрим много члены вида k m Qm,k ( x) = 1 + ci xi + ci xi, (5.89) i =1 i = k + где коэффициенты ci, 1 i k, заданы, а ci, k + 1 i m, свободные.

Обычно ci, 1 i k, определяются из требований аппроксимации. По этому для определенности будем предполагать, что ci = 1 / i !, 1 i k.

Обозначим экстремальные точки (5.89) через x1,…, xm1, причем x1 x2... xm 1. Неизвестные коэффициенты ci, k + 1 i m, опре 5.7. Многочлены устойчивости делим из условия, что многочлен (5.89) в экстремальных точках xi, k i m 1, принимает заданные значения, т. е.

Qm,k ( xi ) = Fi, k i m 1, (5.90) где F ( x) есть некоторая заданная функция, Fi = F ( xi ). Для этого на xi, k i m 1, и ci, k + 1 j m, рассмотрим следующую алгебраи ческую систему уравнений Qm,k ( xi ) = Fi, Qm,k ( xi ) = 0, k i m 1, (5.91) m где Q m,k ( x) = ici xi 1. Перепишем (5.91) в виде, удобном для расче i = тов на ЭВМ. Для этого обозначим через y, z, g и r векторы с коорди натами k yi = xk +i 1, zi = ck +i, gi = Fk +i 1 1 c j yij, j = (5.92) k ri = jc j yij 1, 1 i m k, j = через E1, …, E5 диагональные матрицы с элементами на диагонали вида mk k e1 = k + i, e2 =, e3 = jc j yij 1 + (k + j ) z j yik + j 1, ii ii ii yi j =1 j = mk k j ( j 1)c j yij 2 + (k + j )(k + j 1) z j yik + j 2, ii e4 = (5.93) j =2 j = e5 = ( 1) k +i 1, 1 i m k, ii а через A – матрицу с элементами aij = yik + j, 1 i, j m k. Заметим, что компоненты векторов (5.92), матриц (5.93) и A зависят от чисел m и k, причем g = g ( y ), r = r ( y ), E2 = E2 ( y ), E3 = E3 ( y, z ), E4 = E4 ( y, z ), A = A( y ).

Г л а в а 5. АЛГОРИТМЫ ИНТЕГРИРОВАНИЯ ПЕРЕМЕННОГО ПОРЯДКА И ШАГА Индексы и аргументы опущены для упрощения записи. С использова нием введенных обозначений задачу (5.91) можно записать в виде Az g = 0, E2 AE1 z r = 0. (5.94) Система (5.94) является плохо обусловленной, что приводит к оп ределенным трудностям при использовании для ее решения метода простой итерации. Для сходимости метода Ньютона требуются хоро шие начальные условия, что в данном случае проблемно. Если в (5.91) положить Fi = ( 1)i, k i m 1, т. е. поставить задачу нахождения полинома с максимальным размером интервала устойчивости, то во прос о вычислении начального условия y 0 решается достаточно про сто с использованием значений экстремальных точек многочлена Че бышева, рассматриваемого на отрезке 2m 2, 0, где m есть степень (5.89). Начальные условия можно вычислить по формуле i yi = m 2 cos 1, 1 i m 1. (5.95) m Подставляя (5.95) в первую формулу (5.94), получим коэффициенты полинома Чебышева, для которого | Qm,1 ( x) | 1 при x [2m 2, 0]. При любом k в качестве начальных условий можно взять (5.95) и, как по казывают расчеты, будет хорошая сходимость. Если же F (1)i, k i m 1, то выбор начальных условий является, вообще говоря, искусством.

Опишем один способ решения (5.91) или (5.94), который не нужда ется в хороших начальных условиях. Для численного решения (5.94) используем метод установления, который заключается в том, что для стационарной задачи строится некоторый нестационарный процесс, решение которого с течением времени устанавливается к решению ис ходной задачи. Итак, рассмотрим задачу ( ) dy = E5 E2 AE1 A1 g r, y (0) = y0, (5.96) dt где элементы матрицы E5 определены в (5.93). Ясно, что после опре деления стационарной точки (5.96) коэффициенты многочлена устой 5.7. Многочлены устойчивости чивости можно вычислить из первой системы (5.94). Заметим, что при использовании матрицы E5 все собственные числа матрицы Якоби задачи (5.96) имеют отрицательные вещественные части, т. е. задача (5.96) устойчивая. Из результатов расчетов следует, что (5.96) является жесткой задачей. Методы решения таких задач предполагают вычис ление матрицы Якоби, что в случае (5.96) связано с определенными трудностями. Поэтому для ее решения используем метод второго по рядка точности с численным расчетом и замораживанием матрицы Якоби, который применительно к системе y = f ( y ) имеет вид yn+1 = yn + ak1 + (1 a) k2, Dn k1 = hn f ( yn ), Dn k2 = k1. (5.97) Здесь a = 1 0,5 2;

k1 и k2 стадии метода;

Dn = E ahn An ;

E единичная матрица, hn шаг интегрирования;

An некоторая матри () 2 f n = f ( yn ) / y ца, представимая в виде An = f n + hn Rn + O hn ;

матрица Якоби системы (5.96);

Rn не зависящая от шага интегриро вания произвольная матрица. Если матрица Якоби f n k вычислена k An = f ( yn k ) / y = f n khn Rn + O (hn ), шагов назад, имеем где 0 k I h, I h максимальное число шагов с замороженной матрицей, Rn = f n f n. Если матрица Якоби рассчитывается численно с шагом s j = d j hn, где d j постоянные, то получим, что элементы rn,ij мат rn,ij = 0,5d j 2 fi ( yn ) / y 2, рицы Rn определяются по формулам j 1 i, j N. Так как при записи схемы (5.97) матрица Rn произвольная, то из приведенных выше рассуждений следует, что вопрос о замора живании и численной аппроксимации матрицы Якоби можно рассмат ривать одновременно.

Для контроля точности (5.97) можно применять оценку ошибки ви да n = (a 1 / 3)h 2 f n f n + O (hn ). Так как k2 k1 = ahn f n f n + O( hn ), то 3 2 n с точностью до членов O(hn ) можно оценить по приближенной формуле ( jn ) = a 1 (a 1 / 3) D1 jn [k2 k1 ], 1 jn 2.

n Г л а в а 5. АЛГОРИТМЫ ИНТЕГРИРОВАНИЯ ПЕРЕМЕННОГО ПОРЯДКА И ШАГА В результате для контроля точности схемы (5.97) можно применять неравенство || ( jn ) || 3a / | 3a 1 |, 1 jn 2, (5.98) где – требуемая точность интегрирования;

|| || – некоторая норма в R N ;

целочисленная переменная jn выбирается наименьшей, при которой выполняется неравенство (5.98). Шаг численного диффе sj, 1 j N, ренцирования выбирается по формуле ( ) s j = max 1014, 107 | y j |. Решение жестких задач осуществляется с двойной точностью в силу плохой обусловленности матрицы Якоби системы (5.96). Постоянная 107 введена в формулу с целью выдви j жения ее на середину разрядной сетки. Теперь j -й столбец an матри цы An вычисляется по формуле j an = [ f ( y1,…, y j + s j,…, y N ) f ( y1,…, y j,…, y N )] / s j, 1 j N, т. е. для задания An требуется N вычислений правой части задачи (5.96).

В задаче Коши (5.96) начальные условия для всех m и k выбира ются по формуле (5.95). Заметим, что A1 не вычисляется, а решается линейная система алгебраических уравнений Aw = g. Теперь перейдем к численному исследованию влияния функции F (см. формулу (5.90)) на размер и форму области устойчивости. Сначала докажем следую щую теорему.

Теорема 5.1. Пусть в (5.91) имеет место Fi = (1)i, k i m 1, и пусть | Qm,k ( xi ) | 1 при 1 i k 1. Тогда среди всех полиномов вида (5.89) для построенного многочлена Qm,k ( x) выполняется неравенст во | Qm,k ( x) | 1 на максимальном интервале [ m, 0], где m xm1, Q ( m ) = (1) m.

5.7. Многочлены устойчивости Для доказательства рассмотрим одновременно многочлен Pm,k ( x) вида k m Pm,k ( x) = 1 + ci xi + bi xi, x 0, i =1 i = k + где bi, k + 1 i m, есть некоторый набор коэффициентов. Пусть | Pm,k ( x) | 1 на промежутке [ m, 0], m 0. Предположим, что | Pm,k ( xi ) | 1 при k i m 1 и в то же время m m, где 0 1, a m есть функция от, т. е. m = m (). Рассмотрим Qm,k ( x) и Pm,k ( x) на промежутке [ m (), 0]. Многочлен m Qm,k ( x) Pm,k ( x) = x k +1 (ci bi ) xi k (5.99) i = k + имеет степень меньше m либо равную m. В то же время нуль является его (k + 1) -кратным корнем, и так как имеют место соотношения sign[Qm,k ( xi ) Pm,k ( xi )] = (1)i, k i m 1, sign Qm,k ( m () ) Pm,k ( m () ) = ( 1) m, то многочлен меняет знак между любыми двумя точками xk, …, xm1, m (). Следовательно, на промежутке [ m (), xk ] много член Qm,k ( x) Pm,k ( x) имеет дополнительно (m k ) корней. В ре зультате получим, что из условия m () m следует, что полином (5.99) степени меньше m либо равной m имеет (m + 1) корень. Из по лученного противоречия следует, что m () m. Окончательное до казательство получим, если выберем m как точную верхнюю грань m () по всем. Теорема доказана.

Таким образом, задача о построении метода с максимальным ин тервалом устойчивости сводится к решению (5.91) при условии Г л а в а 5. АЛГОРИТМЫ ИНТЕГРИРОВАНИЯ ПЕРЕМЕННОГО ПОРЯДКА И ШАГА Fi = ( 1)i, k i m 1. Из результатов расчетов следует, что коэффи циент cm многочлена (5.89) убывает с ростом m и, в частности, при m = 13 и k = 2 величина cm имеет порядок 1020. Поэтому примене ние (5.89) при m 13 представляется проблематичным на современ ных ЭВМ.

Результаты численных расчетов Кратко опишем результаты численного эксперимента по определе нию влияния функции F (см. (5.91)) на размер и форму области ус тойчивости. Ниже приведены: m – степень многочлена устойчивости или число стадий метода (5.80), k – число заданных коэффициентов многочлена устойчивости или порядок точности метода, – длина интервала устойчивости, ci, k + 1 i m, – коэффициенты многочлена устойчивости, причем ci = 1 / i! при 1 i k. В расчетах выбрано отно сительно небольшое значение параметра m, равное 4. Сначала для за дания F рассмотрим следующую формулу:

Fi = (1)i u, k i m 1, (5.100) где u – некоторое положительное число. На рис. 5.16–5.28 слева при ведены изолинии | Q( x) | = s при s, равном 1;

0,7;

0,3;

0,1. Справа даны объемные изображения, которые в совокупности с общепринятыми рисунками на плоскости позволяют более глубоко понять структуру области устойчивости и более правильно оценить эффективность чис ленной схемы с точки зрения устойчивости.

Если рассматривать задачу построения многочленов с максималь ным интервалом устойчивости, то нужно положить u = 1. Соответст вующие области устойчивости показаны на рис. 5.16–5.18. На рис. 5. приведена область устойчивости общеизвестного метода типа Рунге– Кутты четвертого порядка точности. Из рис. 5.16–5.18 видно, что об ласти устойчивости методов с максимальным интервалом устойчиво сти «почти» многосвязные. Поэтому, если требуется, чтобы интервал устойчивости был достаточно велик, можно выбрать u = 0,9. Для это го случая области устойчивости приведены на рис. 5.20–5.22.

5.7. Многочлены устойчивости Рис. 5.16. Область устойчивости при значениях параметров:

m = 4;

k = 1;

= 32;

c2 = 0,15625;

c3 = 0,78125 102 ;

c4 = 0,1220703125 Рис. 5.17. Область устойчивости при значениях параметров:

m = 4;

k = 2;

= 12,05;

c3 = 0,78084483452775 101;

c4 = 0,36084539218728 Рис. 5.18. Область устойчивости при значениях параметров:

m = 4;

k = 3;

= 6,03;

c4 = 0,18455702268873 Г л а в а 5. АЛГОРИТМЫ ИНТЕГРИРОВАНИЯ ПЕРЕМЕННОГО ПОРЯДКА И ШАГА Рис. 5.19. Область устойчивости при значениях параметров:

m = 4;

k = 4;

= 2, Рис. 5.20. Область устойчивости при значениях параметров:

m = 4 ;

k = 1 ;

= 29,99 ;

c2 = 0,16491828888770 ;

c3 = 0,87735109261266 102 ;

c4 = 0,14625153854464 Рис. 5.21. Область устойчивости при значениях параметров:

m = 4;

k = 2;

= 11,65;

c3 = 0,80023047470068 10 1;

c4 = 0,38169926858491 5.7. Многочлены устойчивости Рис. 5.22. Область устойчивости при значениях параметров:

m = 4;

k = 3;

= 5,91;

c4 = 0,18738403407411 При u = 0,5 области устойчивости более существенно расширены по мнимой оси, их изображения приведены на рис. 5.23–5.25.

Рис. 5.23. Область устойчивости при значениях параметров:

m = 4 ;

k = 1 ;

= 21,80 ;

c2 = 0, 21254252143474 ;

c3 = 0,15291951589356 101;

c4 = 0,35076567511682 Рис. 5.24. Область устойчивости при значениях параметров:

m = 4;

k = 2;

= 9,75;

c3 = 0,91036842279128 101;

c4 = 0,51553618898240 Г л а в а 5. АЛГОРИТМЫ ИНТЕГРИРОВАНИЯ ПЕРЕМЕННОГО ПОРЯДКА И ШАГА Рис. 5.25. Область устойчивости при значениях параметров:

m = 4;

k = 3;

= 5,33;

c4 = 0, 20266812077634 При u = 0,3 области устойчивости еще больше расширены по мни мой оси, но интервал устойчивости сокращается (см. рис. 5.26–5.28).

Рис. 5.26. Область устойчивости при значениях параметров:

m = 4 ;

k = 1 ;

= 17, 46 ;

c2 = 0, 24956372436363 ;

c3 = 0, 22021550790745 101;

c4 = 0,63043330899815 Рис. 5.27. Область устойчивости при значениях параметров:

m = 4;

k = 2;

= 8,35;

c3 = 0,10176908319354;

c4 = 0,67355258124605 5.7. Многочлены устойчивости Рис. 5.28. Область устойчивости при значениях параметров:

m = 4;

k = 3;

= 4,96;

c4 = 0, 21481671634505 10 При значениях F, равных 0,1;

0,3;

0,1, корни многочлена устойчи вости комплексные. Область устойчивости показана на рис. 5.29.

Рис. 5.29. Область устойчивости при значениях параметров:

m = 4 ;

k = 1 ;

= 10,59 ;

c2 = 0,37216982743909 ;

c3 = 0,52440027420933 101;

c4 = 0, 24749608242196 Отметим, что при уменьшении значения u область устойчивости по вещественной оси сокращается, а по мнимой оси – «раздувается».

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

Из анализа графических изображений областей устойчивости сле дует, что форма, размер и структура области устойчивости зависят от расположения корней многочлена устойчивости (5.89) в комплексной плоскости {h}. На расположение корней можно влиять выбором функции F. При решении задач вида (5.80), собственные числа мат рицы Якоби которых имеют большие мнимые части и решения кото рых носят осциллирующий характер, часто не требуется значительное Г л а в а 5. АЛГОРИТМЫ ИНТЕГРИРОВАНИЯ ПЕРЕМЕННОГО ПОРЯДКА И ШАГА расширение интервала устойчивости. В этом случае шаг из условия точности выбирается достаточно малым и поэтому расширение облас ти устойчивости требуется в основном по мнимой оси. При наличии чисто мнимых собственных чисел нужно, чтобы на некотором участке мнимой оси выполнялось условие | Qm,k ( x) | = 1. При повышении по рядка точности, т. е. с ростом k, это условие выполняется само собой.

Для методов низкого порядка расширения области устойчивости по мнимой оси можно добиться двумя путями, которые можно комбини ровать.

Во-первых, за счет сближения действительных корней полинома устойчивости. При вычислении F по формуле (5.100) это достигается выбором достаточно малого значения параметра u. Однако это приво дит к значительному сокращению интервала устойчивости и, что более существенно, для методов первого порядка точности не приводит к включению мнимой оси в область устойчивости.

Во-вторых, за счет выбора коэффициентов c j, k + 1 j m, таким образом, чтобы корни многочлена устойчивости были комплексными.

При четном m для этого достаточно выбрать функцию F так, чтобы выполнялось неравенство 0 F 1. В случае нечетного m один ко рень обязательно действительный и, в зависимости от того, как он рас положен относительно комплексных корней, изменяются размер, фор ма и структура области устойчивости. В частности, при m = 5, если действительный корень расположен между комплексными корнями, область устойчивости «почти прямоугольная». При значениях F, рав ных 0,1;

0,3;

0,3;

0,1, область устойчивости показана на рис. 5.30.

Рис. 5.30. Область устойчивости при значениях параметров:

m = 5 ;

k = 1 ;

= 17,06 ;

c2 = 0,38097306606596 ;

c3 = 0,59001328331104 101;

c4 = 0,38792470657322 102 ;

c5 = 0,90974734200777 5.7. Многочлены устойчивости При значениях F, равных 0.3;

0,1;

0,3, область устойчивости приведена на рис. 5.31. В этом случае комплексно-сопряженная пара корней находится между вещественными корнями.

Рис. 5.31. Область устойчивости при значениях параметров:

m = 4 ;

k = 1 ;

= 14, 48 ;

c2 = 0, 26137220423406 ;

c3 = 0, 26563274252502 101;

c4 = 0,91728980912922 Отметим одну важную особенность многочленов устойчивости.

При увеличении порядка точности при m = k область устойчивости становится многосвязной. При k = 6 происходит отделение пары ком плексно-сопряженных корней, которые в правой полуплоскости обра зуют «островки устойчивости», вторая пара корней отрывается при k = 11 (рис. 5.32 и 5.33). Представляется, что этот процесс носит регу лярный характер. Становится понятным, почему при повышении по рядка точности область устойчивости не увеличивается достаточно быстро.

Рис. 5.32. Область устойчивости при значениях параметров:

m = 6;

k = 6;

= 3, Г л а в а 5. АЛГОРИТМЫ ИНТЕГРИРОВАНИЯ ПЕРЕМЕННОГО ПОРЯДКА И ШАГА Рис. 5.33. Область устойчивости при значениях параметров:

m = 11;

k = 11;

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

Глава МЕТОДЫ ТИПА РОЗЕНБРОКА П ри решении жестких задач достаточно большой размерности L -устойчивыми методами основные затраты приходятся на обращение матрицы Якоби системы дифференциальных уравнений.

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

Данные методы получили широкое распространение благодаря доста точно хорошим свойствам точности и устойчивости и простоте реали зации на ЭВМ. Идея любого метода типа Розенброка заключается во введении матрицы Якоби непосредственно в численную формулу. Это приводит к необходимости решения на каждом шаге линейной систе мы алгебраических уравнений в отличие от неявных и полуявных ме тодов типа Рунге–Кутты, в которых применяется итерационный про цесс Ньютона для решения нелинейных систем. Введение матрицы Якоби в численную схему существенно упростило программную реа лизацию алгоритмов интегрирования. Однако это привело и к негатив ным последствиям: возникли проблемы с замораживанием матрицы Якоби, т. е. с использованием одной матрицы Якоби на нескольких шагах интегрирования.

6.1. ЧИСЛЕННЫЕ СХЕМЫ Ниже в основном будут рассматриваться численные методы реше ния задачи Коши для автономной системы дифференциальных уравне ний вида y = f ( y ), y (t0 ) = y0, t0 t tk, (6.1) где y и f – вещественные N-мерные вектор-функции;

t – независимая переменная, изменяющаяся на заданном интервале [t0, tk ]. Для реше ния задачи (6.1) будем применять методы типа Розенброка вида i m yn+1 = yn + pi ki, Di ki = hf yn + ij k j, (6.2) i =1 i = Г л а в а 6. МЕТОДЫ ТИПА РОЗЕНБРОКА i где Di = E ai hf yn + ij k j ;

E – единичная матрица;

f = f / y – j = матрица Якоби системы (6.1);

ki, 1 i m, – стадии метода;

pi, ai, ij, ij – коэффициенты, определяющие свойства точности и устойчи вости (6.2). В настоящее время методы типа Розенброка трактуются более широко. Под ними понимаются все численные схемы, в которых матрица Якоби или ее аппроксимация вводится непосредственно в формулу интегрирования. Методы, отличные от (6.2), подробно будут рассматриваться в главе 7.

Численные схемы вида (6.2) можно получить из класса полуявных методов типа Рунге–Кутты, если в них при вычислении каждой стадии ограничиться одной итерацией метода Ньютона. Однако в (6.2) при вычислении стадий необходимо решать только линейные системы ал гебраических уравнений, в то время как в неявных или полуявных ме тодах типа Рунге–Кутты требуется использовать итерационный про цесс типа ньютоновского.

Наиболее эффективные реализации методов типа Розенброка воз никают при выборе некоторых коэффициентов следующим образом:

a1 = a2 = = am = a, ij = 0, 2 i m, 1 j i 1.

Тогда на каждом временном шаге требуется декомпозиция одной мат рицы размерности N, где N есть размерность исходной системы дифференциальных уравнений (6.1). Ниже будем считать это требова ние выполненным.

Дальнейшее сокращение вычислительных затрат достигается за счет применения одной и той же матрицы Якоби на нескольких шагах интегрирования, т. е. за счет использования замороженной матрицы Якоби. Учитывая сделанные замечания, численные схемы (6.2) можно записать в следующем виде:

m yn+1 = yn + pi ki, Dn = E ahAn, a 0, i = (6.3) i Dn ki = hf yn + ij k j, j = где матрица An представима в виде 6.1. Численные схемы An = f n + hBn + O(h 2 ).

(6.4) Здесь f n = f ( yn ) / y – матрица Якоби;

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

В случае замораживания матрицы Якоби имеем An = f ( yn k ) / y = f n khf n f n + O(h 2 ), (6.5) где с учетом (6.4) имеет место Bn = kf n f n, 0 k I f, I f есть мак симальное число шагов с замороженной матрицей Якоби 2 f ( yn ) N fn fn = fi ( yn ).

2 yi i = Если матрица Якоби вычисляется численно с шагом r j = c j h, где c j, 1 j N, – постоянные числа, то получим An = f n + hGn + O(h 2 ), ij Bn = Gn, где элементы g n матрицы Gn определяются по формулам g n = 0,5c j 2 fi ( yn ) / 2 y j, 1 i, j N.

ij В случае замораживания численной матрицы Якоби можно записать An = f n + h(Gn kf n f n ) + O(h 2 ), Bn = Gn kf n f n.

Так как при записи схемы (6.3) матрица Bn является произвольной, то из приведенных выше формул следует, что вопрос о замораживании и численной аппроксимации матрицы Якоби можно рассматривать од новременно.

Ниже потребуется представление приближенного решения yn+1, полученного с помощью схемы (6.3), в виде ряда Тейлора по степеням h до членов с h3 включительно. Введем обозначения i ij, 1 = 0, i = 2 i m, j = (6.6) i ij j, 1 = = 0, = 3 i m.

2 i j = Г л а в а 6. МЕТОДЫ ТИПА РОЗЕНБРОКА С использованием (6.6) стадии ki, 1 i m, схемы (6.3) можно за писать в виде ( ) ki = hf n + (a + i )h 2 f n f n + a 2 + 2ai + h3 f n 2 f n + i + i h3 f n f n2 / 2 + ah3 Bn f n + O(h 4 ), 2 (6.7) где элементарные дифференциалы вычисляются на приближенном ре шении yn. Правильность формулы (6.7) тривиально показывается ме тодом математической индукции. В силу громоздкости выкладок дока зательство не приводится. Подставляя (6.7) в первую формулу (6.3), получим представление приближенного решения в виде ряда Тейлора, т. е.

m m yn +1 = yn + pi hf n + (a + i ) pi h 2 f n f n + i =1 i =1 m m ( ) + a 2 + 2ai + i pi h3 f n2 f n + i pi h3 f n f n2 / 2 + i =1 i =1 m + a pi h3 Bn f n + O(h 4 ). (6.8) i = Далее формула (6.8) будет использоваться при построении конкретных методов интегрирования. Отметим, что если в (6.3) вместо An приме няется матрица Якоби f n, то в (6.8) слагаемое вида h3 Bn f n отсутству ет.

Теорема 6.1. Пусть в численной схеме (6.3) используется матрица An, представимая в виде (6.4), где Bn есть произвольная матрица, не зависящая от величины шага интегрирования. Тогда при любом значе нии параметра m нельзя построить m -стадийную схему (6.3) выше второго порядка точности.

Разложение точного решения y (tn +1 ) в окрестности точки tn в ряд Тейлора по степеням h до членов с h3 включительно имеет вид y (tn+1 ) = y (tn ) + hf + h 2 f / 2 + h3[ f 2 f + f 2 ] / 6 + O(h 4 ), f f (6.9) где элементарные дифференциалы вычислены на точном решении y (tn ). Для доказательства теоремы сравним (6.8) и (6.9) до членов с 6.1. Численные схемы h3 включительно при условии yn = y (tn ). Тогда схема (6.3) будет ап проксимировать задачу (6.1) с третьим порядком, если m m m pi = 1, a pi + i pi =, i =1 i =1 i = m m m a 2 pi + 2a i pi + pi =, (6.10) i i =1 i =2 i = m m i2 pi = 3, a pi = 0.

i =2 i = Так как a 0 по определению, то доказательство теоремы следует из противоречия первого и последнего равенств (6.10). Теорема доказана.

Теорема 6.2. Пусть схема (6.3) используется с замораживанием матрицы Якоби, т. е. матрица An представима в виде (6.4), где Bn имеет вид (6.5). Тогда при любом значении m нельзя построить m стадийную схему (6.3) с постоянными коэффициентами pi, 1 i m, ij, 2 i m, 1 j m 1, выше второго порядка точности.

Для доказательства снова рассмотрим соотношения (6.10). С исполь зованием (6.5) получим, что в этом случае в (6.10) вместо двух послед них уравнений имеем следующее соотношение на коэффициенты:

m 0,5 i pi = 1 / 6 + ak, (6.11) i = где параметр k определяет число шагов с замороженной матрицей Якоби. При получении (6.11) использовалось первое равенство (6.10).

Учитывая (6.6), из (6.11) получим, что либо величина ij, 2 i m, 1 j i 1, либо величина pi, 1 i m, зависит от числа k. Теорема доказана.

Исходя из приведенных выше рассуждений, будем рассматривать методы (6.3) не выше второго порядка точности.

Г л а в а 6. МЕТОДЫ ТИПА РОЗЕНБРОКА 6.2. АЛГОРИТМЫ ИНТЕГРИРОВАНИЯ НА ОСНОВЕ ДВУХСТАДИЙНОЙ ЧИСЛЕННОЙ СХЕМЫ Для решения задачи Коши (6.1) рассмотрим двухстадийную схему вида yn+1 = yn + p1k1 + p2 k2, Dn = E ahf n, (6.12) Dn k1 = hf ( yn ), Dn k2 = hf ( yn + 21k1 ).

Приближенное решение, полученное с помощью формулы (6.12), можно записать с помощью (6.8), т. е.

yn +1 = yn + ( p1 + p2 )hf n + (ap1 + ap2 + 21 p2 )h 2 f n f n + ( ) + h3 2 p2 f n f n2 + 2 a 2 p1 + a 2 p2 + 2a21 p2 f n 2 f n 2 + O(h 4 ). (6.13) 21 Cравнивая (6.9) и (6.13) при условии yn = y (tn ) до членов с h2 включи тельно, получим условия второго порядка точности схемы (6.12), т. е.


p1 + p2 = 1, 21 p2 = 1 / 2 a. (6.14) При записи второго равенства (6.14) использовалось первое.

Исследуем устойчивость (6.12). Применяя ее для решения задачи y = y, y (0) = y0, (6.15) с комплексным, Re( ) 0, получим yn+1 = Q ( x) yn, x = h, где функция устойчивости Q( x) имеет вид Q( x) = {1 + (1 2a ) x + (a 2 2a + 0,5) x 2 } (1 ax) 2. (6.16) Отсюда получим условие L -устойчивости схемы (6.12), т. е.

a 2 2 a + 1 / 2 = 0. (6.17) 6.2. Алгоритмы интегрирования на основе двухстадийной численной схемы Запишем локальную ошибку n,2 схемы (6.12). Сравнивая (6.9) и (6.13) при условии yn = y (tn ) до членов с h3 и учитывая (6.14) и (6.17), получим n,2 = h3 a f 2 f + [21 (2a 1) / 4 + 1 / 6] f 2 + O(h 4 ), f (6.18) где элементарные дифференциалы вычислены на точном решении y (tn ).

Контроль точности (6.12) будем осуществлять методом первого порядка zn +1 = yn + k1, (6.19) где k1 определен в (6.12). Локальная ошибка n,1 схемы (6.19) имеет вид n,1 = (1 / 2 a)h 2 f f + O(h3 ).

(6.20) Из сравнения (6.18) и (6.20) следует, что условия согласования (2.22) локальных ошибок методов первого и второго порядков будут выполнены, если положить равным нулю соотношение на коэффици енты при слагаемом h3 f 2 в формуле (6.18). В результате получим f набор параметров 21 = 2 / (3 6a ), p1 = 1 3(1 2a )2 / 4, p2 = 3(1 2a) 2 / 4, (6.21) при которых (6.12) имеет второй порядок, а ее локальная ошибка имеет вид n,2 = ( a 1 / 3)h3 f 2 f + O (h 4 ). (6.22) Уравнение (6.17) имеет два корня: a1 = 1 0,5 2 и a2 = 1 + 0,5 2. Вы берем a = a1, так как в этом случае меньше коэффициент в главном члене локальной ошибки (6.22). Подставляя a1 в (6.21), получим ( ) a = 1 2 / 2, 21 = 2 2 + 1 / 3, (6.23) ( ) ( ) 2 p1 = 1 3 2 1 / 4, p2 = 3 2 1 4.

Г л а в а 6. МЕТОДЫ ТИПА РОЗЕНБРОКА Для контроля точности схемы (6.12) с коэффициентами (6.23) можно использовать оценки ошибки вида (2.29) или (2.30). Для определения данных оценок рассмотрим разность ( yn +1 zn +1 ), т. е.

yn+1 zn +1 p2 ( k2 k1 ) = (1 / 2 a )h 2 f n f n + O (h3 ).

(6.24) Тогда оценку (2.29) можно вычислить с использованием формул (2.16) и (2.21), где имеет место c1 = 1 / 2 a и c2 = a 1 / 3. Оценка (2.30) мо жет быть вычислена непосредственно по формуле (6.24), т. е. n,2 име ет вид n,2 = 0,5(3a 1)(1 2a )(k2 k1 ). (6.25) Исследуем асимптотическое поведение (6.25). Применяя (6.12) для решения задачи (6.15), получим (6.16). В силу L -устойчивости схемы (6.12) имеет место | Q( x) | 0 при x, т. е. при x прибли женное решение стремится к точному. Так как при решении задачи (6.15) имеет место k2 k1 = 21 x 2 yn / (1 ax) 2, то относительная ошиб ка | n,2 / yn | 2 + 2 при x. Данный факт может привести к неоправданным повторным вычислениям решения (возвратам) вслед ствие невыполнения точности расчетов, хотя фактическая точность может выполняться. Поэтому для контроля точности (6.12) с коэффи циентами (6.23) будем использовать неравенство вида || D1 jn (k2 k1 ) || C, 1 jn 2, (6.26) где параметр jn выбирается наименьшим, при котором выполняется неравенство (6.26), а число C определяется по формуле C = 2 [ (3a 1)(1 2a) ] = 14 2 + 20 39.

Из (6.23) видно, что 21 1. С одной стороны, такой выбор значе ния параметра 21 может привести к улучшению алгоритма контроля точности вычисления решения, потому что, как видно из формул (6.25), он может раньше «почувствовать» изменение производной 6.3. Замораживание матрицы Якоби в методах типа Розенброка решения и вовремя уменьшить величину шага интегрирования. Одна ко, с другой стороны, такой выбор величины параметра 21 приводит к тому, что производная решения вычисляется вне интервала [tn, tn +1 ], а это может ухудшать свойства численной схемы. Итоговая численная схема L -устойчивая, а промежуточная yn+11 = yn + 21k1 (6.27), этим свойством не обладает. Это может приводить к дополнительной ошибке и, как следствие, к понижению эффективности алгоритма интегрирования. Кроме того, свойство внутренней устойчивости – необходимое условие B -сходимости. Потребуем L -устойчивости схе мы (6.27). Применяя ее для решения (6.15), получим yn+11 = [1 + (21 a) x ] yn (1 ax). Отсюда следует, что (6.27) будет, L -устойчивой, если 21 = a. Выберем снова a = 1 0,5 2. Тогда из (6.14) и равенства 21 = a получим набор коэффициентов p1 = 21 = a = 1 2 / 2, p2 = 2 / 2, (6.28) при которых (6.12) имеет второй порядок, а ее локальная ошибка n, имеет вид n,2 = h3[( a 1 / 3) f 2 f + (9a 1) f 2 / 12] + O(h 4 ).

f (6.29) В результате для контроля точности схемы (6.12) с параметрами (6.28) будем применять неравенство (6.26), где параметр C вычисляется по формуле C = 0,5 2 | 3 6a | / | 3a 1 | = 3 2 + 3 7.

6.3. ЗАМОРАЖИВАНИЕ МАТРИЦЫ ЯКОБИ В МЕТОДАХ ТИПА РОЗЕНБРОКА Как уже отмечалось выше, создать метод (6.3) с матрицей An, удовлетворяющей условиям (6.4) и (6.5), с постоянными коэффициен тами выше второго порядка точности нельзя. Поэтому построим алго Г л а в а 6. МЕТОДЫ ТИПА РОЗЕНБРОКА ритм интегрирования с использованием одной и той же матрицы Яко би на нескольких шагах на основе численной схемы второго порядка.

Для этого рассмотрим следующую численную формулу:

yn +1 = yn + p1k1 + p2 k2, Dn i = E ahn i f ( yn i ), (6.30) Dni k1 = hn i f ( yn ), Dn i k2 = hn i f ( yn + 21k1 ), 0 i I f, с коэффициентами p1 = 21 = a = 1 2 / 2, p2 = 2 / 2. (6.31) Здесь hn i – шаг интегрирования;

f = f / y – матрица Якоби исход ной системы дифференциальных уравнений;

I f – целое число, опре деляющее максимальное число шагов с замороженной матрицей;

i – целочисленная переменная, которая может последовательно прини мать значения 0, 1, 2, …, I f.

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

Запишем представление приближенного решения в окрестности точки yn в виде ряда Тейлора до членов с h3 включительно. С ис пользованием (6.5), (6.6) и (6.8) будем иметь yn+1 = yn + ( p1 + p2 )hf n + (ap1 + ap2 + 21 p3 )h 2 f n f n + ( ) + h3 a 2 p1 + a 2 p2 + 2a21 p2 f n 2 f n + ( ) + 2 p2 / 2 iap1 iap2 f n f n2 + O(h 4 ).

(6.32) Поступая по аналогии с построением метода (6.12), (6.28), получим коэффициенты (6.31) L -устойчивой схемы (6.30) второго порядка точ ности, для которой функция внутренней устойчивости является L устойчивой.

6.3. Замораживание матрицы Якоби в методах типа Розенброка Пусть yn = y (tn ). Тогда из сравнения (6.32) с представлением точ ного решения в окрестности точки tn в виде ряда Тейлора y (tn +1 ) = y (tn ) + hf + h 2 f / 2 + ( f 2 f + f 2 ) / 6 + O(h 4 ) f f до членов с h3 включительно следует, что локальную ошибку n+ схемы (6.30) с коэффициентами (6.31) можно представить в виде n +1 = n +1 + n +1, где 1 3 2 n +1 = a h f n f n + O (h ), 1 1 n+1 = + ia + a 21 h f n f n. (6.33) 6 2 4 При записи (6.33) использовались соотношения (6.14) и (6.17).

С применением метода первого порядка (6.19) и соотношения (6.24) получим, что оценка ошибки n имеет вид ( ) n = 2 1 (k2 k1 ) / 3. (6.34) Ниже для схемы (6.30) будет построена оценка ошибки типа (6.34), причем первые члены разложений в ряды Тейлора построенной оценки и (6.34) будут совпадать. Чтобы воспользоваться оценкой типа (6.34), нужно добиться выполнения равенства n+1 = 0, что можно сделать двумя способами.

1. Выбрать 21 из условия равенства нулю соотношения на коэф фициенты в выражении для n+1.

2. Используя уже сделанные вычисления – чтобы не увеличивать затраты на шаг интегрирования – перестроить схему (6.30).

В первом случае получим 21 = (2 + 12ai ) / (3 6a) и 21 при i, что, вообще говоря, приводит к вычислению функции f за пределами шага интегрирования. Поэтому остановимся на прежнем вы боре 21. Для изучения второго случая обозначим k3 = hn i f ( yn + ak1 ) и разложим k3 в ряд Тейлора Г л а в а 6. МЕТОДЫ ТИПА РОЗЕНБРОКА k3 = hf n + ah 2 f n f n + a 2 h3 f n 2 f n + a 2 h3 f n f n2 / 2 + O(h3 ).

(6.35) С использованием (6.14) из (6.7) и (6.35) легко видеть, что k3 k1 = ( a 2 / 2 + ia)h3 f n f n2 + O(h 4 ), и схему (6.30) можно преобразовать к виду yn+1 = yn + ak1 + (1 a )k2 + p (i )( k3 k1 ), (6.36) где p (i ) = 1 + 6ia + 3a 3 2a 2 ( 3a 2 + 6ia ). Для этой схемы очевидно n+1 = 0, она обладает вторым порядком аппроксимации и, так как на линейной задаче внесенная добавка равна нулю, формально является L -устойчивой.

Фактически же схему (6.36) L -устойчивой считать нельзя, потому что при больших шагах интегрирования происходит настолько ощути мое усиление ошибок округлений, что (6.36) становится непригодной для практических расчетов жестких задач. Причина этого – в примене нии неитерированного значения функции f. Такой вывод подтвер ждают расчеты даже простейших линейных уравнений. Указанный недостаток можно устранить, рассмотрев схему k yn+1 = yn + ak1 + (1 a) k2 + p(i ) Dn i (k3 k1 ), (6.37) причем по отношению к ошибкам округлений при k = 1 и k = 2 схема (6.37) является соответственно A -устойчивой и L -устойчивой. Однако это приводит к одному или двум дополнительным обратным ходам в методе Гаусса. Чтобы в некоторой мере компенсировать затраты, бу дем использовать схему (6.30), а при выборе величины шага интегри рования одновременно с контролем ошибки типа (6.34) будем требо вать выполнения неравенства 1 j || Dn i n (k3 k1 ) || C f, 1 jn J f, (6.38) где || || – некоторая норма в R N, J f – целое число, jn – целочислен ная переменная, – требуемая точность расчетов, C f – постоянная.

6.3. Замораживание матрицы Якоби в методах типа Розенброка Очевидно, что при jn = 3 вычислительные затраты на шаг интег рирования совпадают с (6.37) при k = 2. Нетрудно понять также, что одновременное использование (6.34) и (6.38) эквивалентно контролю оценки локальной ошибки n+1 и контролю аналога глобальной ошибки, построенной по n+1. Можно оценить C f, например, поло жить C f равным C f = 0,5/ | p(i ) |, хотя, по нашему мнению, эту кон станту следует выбирать при настройке алгоритма на конкретный тип задач из условия минимизации повторных вычислений решения (воз вратов) с помощью численного эксперимента. Аналогичное замечание справедливо и относительно выбора J f.


Теперь перейдем к получению неравенства для контроля точности вычислений. Прежде всего отметим, что, так как 21 0,3, то в фор муле для оценки ошибки (6.34) не используется информация в точке t = tn +1, что для задач с большой жесткостью может приводить к поте ре точности вычислений. В данном случае ситуация усугубляется при менением «испорченной» матрицы Якоби. Чтобы избавиться от этого недостатка, используем еще одно вычисление функции f. Обозначим k4 = hni Dni f ( yn+1 ) и разложим k4 в ряд Тейлора в окрестности точ ки yn до членов с h 2 включительно, т. е.

k4 = hf n + (1 + a )h 2 f n f n + O( h3 ).

(6.39) Теперь вместо (6.34) можно применять следующую оценку ошибки:

n = ( 3 2 4 ) ( k4 k1 ) 6, (6.40) причем с использованием (6.7) и (6.39) нетрудно убедиться, что первые члены разложений в ряды Тейлора ошибок (6.34) и (6.40) совпадают.

Отметим, что оценка (6.34) может использоваться как предварительная с целью уменьшения возможных возвратов. В конечном итоге для кон троля точности и при выборе величины шага схемы (6.30) с коэффици ентами (6.31) будем применять неравенство (6.38) и неравенство вида 1l || Dn in (k4 k1 ) || 3(3 2 + 4), 1 ln 2. (6.41) Г л а в а 6. МЕТОДЫ ТИПА РОЗЕНБРОКА Заметим, что в практических расчетах неравенство (6.38), как пра вило, выполняется при jn = 1, а jn 1 в основном имеет место либо при прохождении экстремальных точек, либо при резком возрастании величины шага интегрирования. Далее, так как мы используем более сильное ограничение на выбор величины шага – одновременный кон троль неравенств (6.38) и (6.41), чем при контроле точности – контроль только оценки (6.41), то размер шага может уменьшаться даже в том случае, если требуемая точность выполняется. В результате подходя щим выбором C f и J f можно фактически избежать повторных вы числений решения. Кроме того, дополнительное вычисление функции f ( yn+1 ), применяемое для оценки ошибки (6.40), в случае успешного шага к увеличению вычислительных затрат не приводит, потому что оно используется для выполнения следующего шага. Более того, в слу чае расчетов с замороженной матрицей k4 совпадает со стадией k1, которая необходима для выполнения (n + 1) -го шага. Отметим также, что в случае линейной системы дифференциальных уравнений с по стоянными коэффициентами левая часть неравенства (6.38) равна ну лю, и значит, данное неравенство можно использовать в алгоритме ин тегрирования для определения типа задачи – линейная или нелинейная.

Теперь сформулируем алгоритм интегрирования на основе схемы (6.30) с коэффициентами (6.31). Для упрощения описания предвари тельная оценка не используется, хотя в программной реализации алго ритма осуществляется контроль дополнительного неравенства с при влечением (6.34). Введем обозначения:

1l Vn (ln ) = Dn in ( k4 k1 ), 1 ln 2, 1 j Wn ( jn ) = Dn i n (k3 k1 ), 1 jn J f.

Если задача линейная, то J f следует положить равным нулю, и кон троль неравенства (6.38) осуществляться не будет.

Итак, с учетом соотношения Vn (ln ) = O(h 2 ), Wn ( jn ) = O(h3 ) алго ритм интегрирования формулируется следующим образом.

6.3. Замораживание матрицы Якоби в методах типа Розенброка Шаг 1. Решение задачи (6.1) из точки tn в точку tn+1 вычисляется по формуле (6.30) с коэффициентами (6.31).

Шаг 2. Если J f = 0, то Wn полагается равным нулю и управление передается на шаг 4.

Шаг 3. Среди чисел 1,..., J f выбирается такое наименьшее jn, что выполняется неравенство || Wn ( jn ) || C f. Если данное неравенство не выполняется ни при каком значении jn, то jn полагается равным J f.

Шаг 4. Из чисел 1 и 2 выбирается такое наименьшее число ln, что выполняется неравенство || Vn (ln ) || (12 + 9 2 ). Если данное нера венство не выполняется ни при каком значении ln, то ln полагается равным 2.

Шаг 5. Вычисляется значение параметра r по формуле r 3 || Wn ( jn ) || = C f.

s 2 || Vn ( mn ) || = Шаг 6. Вычисляется параметр s по формуле = (12 + 9 2).

Шаг 7. Вычисляется значение параметра hnew по формуле hnew = min( s, r )hn.

Шаг 8. Если s 0, то параметр i полагается равным нулю (матри ца Якоби размораживается) и происходит повторное вычисление ре шения (возврат) с шагом hn = hnew.

Шаг 9. Значение i полагается равным i + 1.

Шаг 10. Если hn hnew H f hn, i I f, где H f – постоянная, H f 1, то вычисляется решение в следующей точке при hn +1 = hn с замороженной матрицей Якоби. Иначе i полагается равным нулю (матрица Якоби размораживается) и решение в следующей точке вы числяется при hn +1 = hnew.

Эффективность построенного алгоритма зависит от набора пара метров C f, J f, I f и H f, выбором которых можно настроить его на конкретный тип задач. Вопрос о выборе C f и J f обсуждался выше.

Г л а в а 6. МЕТОДЫ ТИПА РОЗЕНБРОКА Заметим здесь только, что если исходная задача мало изучена, то J f следует положить равным 3. Параметры I f и H f выбираются в зави симости от отношения стоимости вычисления функции f к стоимости вычисления и обращения матрицы Якоби. При увеличении I f и H f количество вычислений правой части возрастает, а число вычислений и обращений матрицы Якоби убывает. При I f = 0 или H f = 0 матри ца Якоби вычисляется на каждом шаге интегрирования. В конкретных расчетах алгоритм успешно использовался со значениями параметров C f = 1, 4, J f = 3, I f = 20 и H f = 1,6.

Глава КЛАСС (m, k)-МЕТОДОВ Т рудности с замораживанием матрицы Якоби в методах типа Розенброка привели к поискам их модификаций. Основное отличие рассматриваемых здесь методов заключается в том, что в них стадия метода не связывается с обязательным вычислением правой части исходной задачи. За счет этого проблема использования одной матрицы Якоби на нескольких шагах упрощается. В данной главе ис следуются (m, k ) -методы с одним (k = 1) и двумя (k = 2) вычисления ми правой части исходной задачи. Показано, что при k, равном 1 и 2, можно построить (m, k ) -методы соответственно второго и четвертого порядков точности, и получены соответствующие коэффициенты чис ленных схем. При k = 1 подробно исследован случай, когда (m,1) методы применяются для решения линейных задач с постоянными ко эффициентами.

7.1. ЧИСЛЕННЫЕ СХЕМЫ Пусть Z есть множество целых чисел и заданы m, k Z, k m.

Обозначим через M m множество чисел i Z таких, что 1 i m, а через M k и J i – подмножества из M m вида M k = {mi M m | m1 = 1, mi 1 mi, 2 i k, mk m}, (7.1) J i = {m j 1 M m | j 1, m j M k, m j i}, 1 i m.

Для численного решения задачи Коши y = f ( y ), y (t0 ) = y0, t0 t tk, (7.2) рассмотрим численные схемы следующего вида:

m yn+1 = yn + pi ki, Dn = E ahf n, i = Г л а в а 7. КЛАСС (m, k)-МЕТОДОВ i 1 i Dn ki = hf yn + ij k j + ij k j + hf n cij k j, i M k, (7.3) jJ j =1 j = i i ij k j + hf n cij k j, i M m Dn ki = ki 1 + Mk, jJ i j = где a, pi, ij, ij и cij – постоянные коэффициенты, h – шаг интег рирования, E – единичная матрица, f = f / y – матрица Якоби, k – количество вычислений функции f, m – количество обратных ходов в методе Гаусса. На каждом шаге интегрирования осуществляются од но вычисление матрицы Якоби и одна декомпозиция матрицы Dn. До пускается аппроксимация матрицы Якоби матрицей An, удовлетво ряющей условию (6.4). Так как k и m определяют затраты на шаг, а числа m1,..., mk из M k только распределяют их внутри шага, то ме тоды типа (7.3) называются (m, k ) -методами. Отметим, что при k = m и ij = cij = 0 (m, k ) -методы совпадают со схемами типа Розенброка (6.2), а при k = m и ij = 0 – с ROW-методами. Для (m, k ) -методов определены затраты на шаг интегрирования и правильно описана об ласть определения коэффициентов численных формул, что упрощает их исследование.

Заметим также, что при рассмотрении методов такого типа, как правило, все авторы изучали случай k = m, т. е. когда число стадий и количество вычислений правой части совпадают. В этом случае k стадийную схему (7.3) можно поставить в соответствие k -стадийной полуявной формуле типа Рунге–Кутты, при реализации которой на ка ждом шаге используется одна матрица размерности N. Относительно таких численных формул известно, что нельзя построить k -стадийную схему выше (k + 1) -го порядка точности. Очевидно, этот факт распро страняется и на (m, k ) -методы при m = k. В то же время методы типа Розенброка и полуявные методы типа Рунге–Кутты максимального порядка при k, равном 1, 2 и 3, построены. По-видимому, все это по служило причиной того, что методы (7.3) при k m не изучались. Если рассматривать схемы (7.3) при m k, то при k, равном 1 и 2, можно 7.2. Ряды Тейлора для стадий методов построить такие численные формулы, что их порядок точности равен 2k. Таким образом, среди методов (7.3) имеются численные схемы, которые в смысле точности не хуже неявных методов типа Рунге– Кутты и в то же время они требуют существенно меньших вычисли тельных затрат при реализации.

7.2. РЯДЫ ТЕЙЛОРА ДЛЯ СТАДИЙ МЕТОДОВ При исследовании (m, k ) -методов (7.3) потребуются представления в виде рядов Тейлора в окрестности точки yn до членов с h 4 включи тельно стадий вида Dn k1 = hf ( yn ), Dn k2 = k1, Dn k3 = hf ( yn + 31k1 + 32 k2 ) + 32 k2, Dn k4 = k3 + 42 k2 (7.4) и следующего:

Dn k1 = hf ( yn ), Dn k2 = hf ( yn + 21k1 ) + 21k1, Dn k3 = k2 + 31k1, Dn k4 = k3 + 41k1. (7.5) Для случая (7.4) имеем k1 = hf n + ah 2 f n f n + a 2 h3 f n2 f n + a3h 4 f n3 f n + O( h5 ), k2 = hf n + 2ah 2 f n f n + 3a 2 h3 f n 2 f n + 4a3h 4 f n3 f n + O( h5 ), k3 = (1 + 21 )hf n + (a + 31 + 32 + 3a32 )h 2 f n f n + ( ) + h3 a 2 + 2a31 + 3a32 + 6a 2 32 f n 2 f n + (31 + 32 )2 f n f n2 / 2 + ( ) + h 4 a3 + 3a 231 + 6a 232 + 10a 332 f n3 f n + + a (31 + 32 )(31 + 232 ) f n f n f n2 + (7.6) +10a332 f n3 f n + a (31 + 32 )(31 + 232 ) f n f n f n2 + + a(31 + 32 ) 2 f n f n f n2 / 2 + (31 + 32 )3 f n f n / 6 + O(h5 ), Г л а в а 7. КЛАСС (m, k)-МЕТОДОВ k4 = (1 + 32 + 42 )hf n + (2a + 31 + 32 + 4a32 + 3a 42 )h 2 f n f n + ( ) + h3 3a 2 + 3a31 + 4a32 + 10a 2 32 + 6a 2 42 f n2 f n + +0,5 ( 31 + 32 ) f n f n2 + ( ) + h 4 4a3 + 6a 231 + 10a332 + 10a3 42 f n3 f n + a(31 + 32 )(31 + 232 ) f n f n f n2 + + a (31 + 32 ) 2 f n f n f n2 + (31 + 32 )3 f n f n / 6 + O (h5 ).

Для стадий ki, 1 i 4, вычисляемых по формулам (7.5), получим k1 = hf n + ah 2 f n f n + a 2 h3 f n2 f n + a3h 4 f n3 f n + O( h5 ), k2 = (1 + 21 )hf n + (a + 21 + 2a 21 )h 2 f n f n + ( ) + h3 a 2 + 2a21 + 3a 2 21 f n 2 f n + 0,521 f n f n2 + ( ) h 4 a3 + 3a 221 + 4a3 21 f n3 f n + a2 f n f n f n2 + 0,5a2 f n f n f n2 + 21 + 3 f n f n / 6 + O (h5 ), k3 = (1 + 21 + 31 )hf n + ( +(2a + 21 + 3a 21 + 2a31 )h 2 f n f n + h3 3a 2 + 2a21 + 6a 2 21 + ) ( +3a 2 31 f 2 f n + 0,521 f n f n2 + h 4 4a3 + 6a 221 + 10a3 21 + (7.7) n ) +4a331 f n3 f n + a21 f n f n f n2 + a21 f n f n f n2 + 3 f n f n / 6 + O (h5 ), 2 2 21 k4 = (1 + 21 + 31 + 41 )hf n + (3a + 21 + 4a 21 + ( +3a31 + 2a 41 )h 2 f n f n + h3 6a 2 + 4a21 + 10a 2 21 + ) +6a 231 + 3a 2 41 f n 2 f n + 0,521 f n f n2 + 7.3. Численные схемы с одним вычислением правой части ( ) + h 4 10a3 + 10a 221 + 20a 2 21 + 10a 231 + 4a 2 41 f n3 f n + +1,5a2 f n f n f n2 + a21 f n f n f n2 + 3 f n f n2 / 6 + O( h5 ).

21 21 Ниже соотношения (7.4)–(7.7) будут использоваться при построе нии конкретных методов. Исследование (m, k ) -методов проведем по следующей схеме.

1. Выбираем числа k, m и p, где p есть порядок точности мето да. Задаемся множеством M k и записываем J i, которое теперь опре деляется однозначно.

2. Выбираем нужные представления ki из (7.6) или (7.7) и подстав ляем в (7.3).

3. При yn = y (tn ) сравниваем приближенное решение с рядом Тей лора (2.31) для точного решения до членов с h p включительно. Запи сываем соотношения на коэффициенты, обеспечивающие заданный порядок точности.

4. Применяем (7.3) к тестовому уравнению y = y, Re() 0, и за писываем соотношение на коэффициенты, обеспечивающее свойство L -устойчивости.

5. Исследуем совместность полученной нелинейной системы ал гебраических уравнений. Если система несовместна, то либо при дан ных m и k нельзя построить L -устойчивый метод p -го порядка, ли бо неудачно выбрано M k.

7.3. ЧИСЛЕННЫЕ СХЕМЫ С ОДНИМ ВЫЧИСЛЕНИЕМ ПРАВОЙ ЧАСТИ Выберем множества (7.1) вида M1 = {1}, J i =, 2 i m, т. е. рас смотрим методы m yn+1 = yn + pi ki, Dn k1 = hf ( yn ), Dn ki = ki 1, 2 i m, (7.8) i = где матрица Dn определяется второй формулой (7.3).

Г л а в а 7. КЛАСС (m, k)-МЕТОДОВ 7.4. ОБЩИЕ ПОЛОЖЕНИЯ Для изучения методов (7.8) введем в рассмотрение матрицу B b1i = bi1 = 1, i 1, (7.9) bij = bi 1, j + bi, j 1, i, j 2. (7.10) Лемма 7.1. Элементы матрицы B представимы в виде j i bi 1,k, bij = bk, j 1, bij = i, j 2. (7.11) k =1 k = Для доказательства достаточно расписать рекуррентное соотноше ние (7.10) и использовать условие (7.9).

Ниже через Bs,k будем обозначать матрицу с элементами (7.9) и (7.10), составленную из первых s строк и k столбцов матрицы B.

Лемма 7.2. Матрица Bm m невырожденная.

Для доказательства введем в рассмотрение матрицы Rk, 2 k m, у которых все элементы равны нулю, за исключением следующих:

rk = 1, 1 i m, rk,i 1 = 1, k i m, ii i (7.12) и матрицу R вида R = Rm Rm 1... R2. Очевидно, что матрица R невы рожденная. Покажем, что RBm,m есть верхняя треугольная матрица с единицами на главной диагонали, что доказывает лемму 7.2. Сначала умножим R2 на Bm,m. Учитывая (7.12), нужно из второй строки мат рицы Bm,m вычесть первую, из третьей – вторую и т. д., а результат следует поместить на место соответственно второй строкой, третьей и т. д. Тогда с использованием (7.9) получим, что в первом столбце мат рицы R2 Bm,m на первом месте стоит единица, а на остальных – нули.

Далее, из (7.10) имеем, что bi, j 1 = bij bi 1, j, 2 i m, 3 j m, т. е.

если в матрице R2 Bm,m вычеркнуть первую строку и столбец, то она совпадает с Bm,m, у которой вычеркнуты первая строка и последний столбец. Умножая последовательно матрицу R2 Bm,m на R3,..., Rm и проводя аналогичные рассуждения, получим доказательство леммы 7.2.

7.4. Общие положения Заметим, что так как матрицы Bm,m и Rk, 2 k m, целочислен ные, то R и RBm,m также целочисленные.

Лемма 7.3. Стадии метода (7.8) представимы в виде k j = ai 1bij hi f n(i 1) f n.

(7.13) i = Доказательство проведем методом математической индукции. Для этого запишем Dn 1 в виде ряда Тейлора по степеням h, т. е.

Dn 1 = ai hi f ni.

(7.14) i = Представление (7.13) при j, равном единице, очевидно. Пусть (7.13) выполняется при некотором j. Тогда для определения kn, j +1 умно жим (7.14) на (7.13) и соберем слагаемые при одинаковых степенях h.

В результате будем иметь i kn, j +1 = ai 1 bkj hi f n(i 1) f n.

k = i = Окончательно получим доказательство, если воспользуемся вторым соотношением (7.11). Лемма доказана.

С применением (7.13) приближенное решение, полученное по схе ме (7.8), можно представить в виде m yn+1 = yn + ai 1 bij p j hi f n(i 1) f n.

(7.15) j =1 i =1 Теорема 7.1. При любом значении m нельзя построить m стадийную схему (7.8) выше второго порядка точности.

Для доказательства положим yn = y (tn ). Тогда из сравнения ряда Тейлора для точного решения y (tn +1 ) и (7.15) следует, что в разложе нии точного решения в ряд Тейлора есть слагаемое вида h3 f 2, а в f представлении (7.15) оно отсутствует. Теорема доказана.

Г л а в а 7. КЛАСС (m, k)-МЕТОДОВ 7.5. МЕТОДЫ РЕШЕНИЯ ЛИНЕЙНЫХ ЗАДАЧ Для линейной задачи с постоянными коэффициентами y = Ay + b, y (t0 ) = y0, t0 t tk, (7.16) можно построить методы вида (7.8) более высокого порядка точности.

Точное решение y (tn +1 ) задачи (7.16) представимо в виде y (tn+1 ) = y (tn ) + hi f (i 1) f / i !, (7.17) i = где f = A – матрица Якоби задачи (7.16). Введем обозначения:

Pm = ( p1,..., pm )T, Vm (a) = (1, a 1 / 2!,..., a1m / m!)T. (7.18) Пусть yn = y (tn ). Тогда, сравнивая (7.15) и (7.17) до членов с h p вклю чительно, получим условия p -го порядка точности схемы (7.8), т. е.

m ai 1 bij p j = 1 / i !, 1 i p, (7.19) j = при этом ее локальная ошибка n, p имеет вид n, p = p h p +1 f p f + O(h p+ 2 ), (7.20) m p = 1 / ( p + 1)! a p b p +1, j p j. (7.21) j = Теорема 7.2. Пусть метод (7.8) применяется для решения линей ной задачи (7.16). Тогда при любом значении m можно построить m стадийную схему p -го порядка точности, p = m + 1.

Для доказательства рассмотрим соотношения (7.19) при p = m + 1.

С использованием (7.18) нетрудно видеть, что (7.19) можно переписать в виде m bm+1, j p j = 1 / [a m (m + 1)!].

Bm m Pm = Vm (a ), (7.22) j = Система (7.22) состоит из (m + 1) -го уравнения относительно (m + 1) неизвестных параметров a и pi, 1 i m. Учитывая, что мат рица Bm m невырожденная, определим из первого соотношения (7.22) 7.5. Методы решения линейных задач вектор Pm и подставим во второе равенство (7.22). В результате полу чим уравнение m -й степени относительно a. При фиксированном m данное уравнение может иметь m различных корней. Определяя один из них, значения pi, 1 i m, получим из линейной системы (7.22), что и завершает доказательство.

В (7.12) все операции над матрицей Bm,m можно проводить в це лой арифметике, что позволяет определить pi, 1 i m, с высокой точностью.

A-устойчивые методы (m + 1)-го порядка При четном m коэффициенты m -стадийной схемы (m + 1) -го по рядка точности могут получаться комплексными. Из анализа результа тов расчетов следует, что уравнение относительно параметра a имеет действительные корни. Ниже для каждого m, 1 m 10, приведено по одному значению параметра a, a = a (m), полученного методом деле ния отрезка пополам:

a(1) = 0,5;

a(2) = 0,7886151346;

a(3) = 1,0685790213;

a(4) = 1,3453664198;

a(5) = 0, 4732683913;

a(6) = 0,5566999421;

a(7) = 0,6395554058;

a(8) = 0,7220408535;

(7.23) a(9) = 0,8042736126;

a(10) = 0, 4173849297, при которых (7.8) с коэффициентами (7.22) имеет (m + 1) -й порядок.

Для решения линейных задач (7.16) можно найти коэффициенты явных m -стадийных методов типа Рунге–Кутты m -го порядка точно сти, т. е. там повышение порядка напрямую связано с ростом числа вычислений правой части исходной задачи. При применении (7.8) по вышение порядка точности приводит к увеличению числа обратных ходов в методе Гаусса. В случае достаточно большой размерности за дачи основные затраты в методах (7.8) приходятся на выполнение декомпозиции матрицы Dn. В то же время при решении жестких задач (7.16) схемы (7.8) имеют преимущество по сравнению с явными мето дами, так как они могут обладать свойством A -устойчивости.

Г л а в а 7. КЛАСС (m, k)-МЕТОДОВ Исследуем устойчивость численной формулы (7.8). Применяя схему (7.8) для решения линейного скалярного уравнения y = y, y (0) = y0, Re( ) 0, t 0, получим yn+1 = Q( x) yn, x = h, где m Q ( x) = 1 + pi x / (1 ax)i. (7.24) i = Переходя к пределу при x, имеем | Q( x) | | p1 a | / | a |. Отсю да при a 0 необходимое условие A -устойчивости схемы (7.8) имеет вид | p1 a | / a 1 или 0 p1 2a. (7.25) Непосредственно проверкой нетрудно убедиться, что при использова нии (7.23) коэффициенты схемы (7.8) удовлетворяют условию (7.25).



Pages:     | 1 |   ...   | 2 | 3 || 5 | 6 |   ...   | 9 |
 





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

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