Разработка и тестирование методов решения жестких обыкновенных дифференциальных уравнений Текст научной статьи по специальности «Математика»

УДК 519.63

Разработка и тестирование методов решения жестких

обыкновенных дифференциальных уравнений

12 2 © М.П. Галанин ' , С.Р. Ходжаева

1 Институт прикладной математики им. М.В. Келдыша РАН, Москва, 125047, Россия

2 Московский государственный технический университет им. Н.Э. Баумана, Москва, 105005, Россия

Приведены исследования (m^-метода, одностадийной комплексной схемы Розен-брока, метода конечных суперэлементов и явного четырехстадийного метода Рунге — Кутты применительно к решению жестких систем обыкновенных дифференциальных уравнений. Анализ тестовых расчетов показал, что лучшим выбором для систем с большим числом жесткости является одностадийная комплексная схема Розенброка (CROS). Метод конечных суперэлементов (МКСЭ) является «точным» для решения линейных систем обыкновенных дифференциальных уравнений, лучшим вспомогательным методом для его реализации является (4,2)-метод. Построен и протестирован вариант метода конечных суперэлементов для решения нелинейных задач, оказавшийся непригодным для задач большой жесткости.

Ключевые слова: обыкновенные дифференциальные уравнения, численное решение, жесткие системы, метод конечных суперэлементов, (4,2)-метод, CROS.

Введение. Постановка задачи. Рассмотрим задачу Коши для автономной системы обыкновенных дифференциальных уравнений (ОДУ):

u' = f(u), u(0) = U0, 0 < t < t0, (1)

T

где u — вектор-столбец \u1(t),..., un(t)} , причем известно, что данная система является жесткой.

Далее будут встречаться и неавтономные системы. Будем считать линейную систему ОДУ u' = Au (A — постоянная матрица n х n ) жесткой, если выполняются два следующих условия [1]:

1) все собственные числа Xi матрицы A имеют отрицательную действительную часть, т. е. ReXi < 0, i = 1,2,..., n;

2) число

maxlRe X J

S = 1<k <n' k|

min IRe X J

1<k <n k

велико. Число S называют жесткостью задачи.

Если A = A(t) и Xk =Xk (t), то вводят понятие жесткости на временном интервале. В этом случае число sup S(t) должно быть

/6(0,/q)

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

Не все задачи, называемые в литературе жесткими, соответствуют данному определению. Часто жесткими называют все задачи, решение которых содержит компоненты с резко различными характерными временами изменения [2].

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

В настоящей работе проведено сравнение численных методов решения жестких систем ОДУ: (4,2)-метода, МКСЭ, метода Рунге — Кутты четвертого порядка точности и одностадийной комплексной схемы Розенброка. Сравнение выполнено на основе решения тестовых задач. Для оценки работы методов выбран набор тестов, обладающих разной степенью жесткости [3]. Предложен и протестирован вариант МКСЭ для нелинейных задач.

Методы решения жестких систем ОДУ.

1. (ш,к)-метод [4]. Пусть заданы целые положительные числа т и к, к < т . Обозначим через Мт множество целых чисел г, удовлетворяющих условию 1 < г < т, а через Мк и Ji — подмножества из Мт следующего вида [4]:

Мк ={тг е Мт I1 = т1 < т2 < .. < тк < т);

Ji =\т^_1 е Мт | ] > 1, т] е Мк, mj < /}, 2 < г < т.

Тогда семейство (т, к )-методов записывается в виде

т

Уи+1 = У п + Е Ргкг; = Е _ аТп; г=1

Г г _1 1

Ъпкг = тГ Уп к3 + Еаз к j, г е Мк;

V 3=1 У

°пкг = кг_1 + Е а3к3, г е Мт \ Мк,

где т — шаг интегрирования решаемой задачи; у п — приближенное решение при / = I п; Е — единичная матрица размерностью п; {'п = д(у п)/ ду — матрица Якоби векторной функции {(у) а, рг, а3 и Ргз — вещественные коэффициенты, определяющие свойства точности и устойчивости (т, к) -метода.

Рассмотрим реализованный в данной работе Ь -устойчивый (4,2)-метод четвертого порядка точности [4]. Он выглядит следующим образом:

4

у и+1 = у пр»к», =Е - атЪп;

г=1

°пк1 = тГ (у п); °пк 2 = к1; Опк з = тг (у п+Рэ1к1 +Рэ2к 2)+аэ2к 2;

Опк 4 = к з +а,42к 2,

а коэффициенты а, р1, а^ 2, Р3к, 1 < I < 4, 3 < у < 4, 1 < к < 2 имеют следующие значения [4]: а = 0,57281606248213; р = 1,27836939012447; р2 =-1,00738680980438; р3 = 0,92655391093950; р4 =-0,33396131834691; р31 = 1,00900469029922; р32 =-0,25900469029921; а32 =-0,49552206416578; а42 =-1,28777648233922.

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

2. Одностадийная комплексная схема Розенброка (СЯ08). Общий вид многостадийных схем семейства Розенброка [5] имеет следующий вид:

Е -™тт*ы

т-1

У п + ТЕ Сткwк, гп + ТСт

к=1

5 ^ Ьтw т;

т=1

Г

w = Ъ т У

V

т-1

атк™к , + тат

к=1

где Ъы — матрица Якоби решаемой задачи; ^ — число стадий метода; атк, ат, Ьт, стк, ст — заданные комплексные (вообще говоря) числа.

Формулы перехода на новый временной слой для однопарамет-рического семейства одностадийных схем Розенброка имеют следующий вид:

Уп+1 = Уп + тКек; (Е-рт^(уп,г))к = Ъ(уп,г + та1), где Яе к — действительная часть вектора к с комплексными координатами; Е — единичная матрица; т — шаг метода; Р и а1 — числовые параметры, определяющие свойства схемы.

Данная схема обладает точностью 0(т2) при Яеа1 = ЯеР =

1

1

поэтому выберем а1 = и приведем исследование свойств системы в

зависимости от значения коэффициента р.

Выясним, при каких значениях параметра Р данная схема является ^-устойчивой (приведенные далее рассуждения основаны на [6-10]).

Функция устойчивости одностадийной комплексной схемы Ро-зенброка имеет вид [5]:

т=1±(ЬМ, 1 -ре,

а по определению схема является А-устойчивой, если |Я(£)| < 1 при Яе(^) < 0. Пусть Р = а + 7Ь, ^ = + 7^2. Запишем условие ^-устойчивости и получим допустимые значения р :

1 + (1 - а - ¡Ь)(^1 + 7^)

< 1.

1 - (а + /Ь)(^1 + 7^) После преобразований получим следующее неравенство:

( +^2 )-2а (^2 +^2 )<-2^1. Поскольку < 0, окончательно приходим к условию ^-устойчивости для СКОБ: Яе(р) > 2.

Напомним, что схема является /р-устойчивой, если она А-

устойчива и выполняется следующее условие: Я= р), р > 0

прии Яе0 [9].

Комплексная схема Розенброка является / 1-устойчивой, если параметр р лежит на правой половине окружности с центром в точке

Г10 1 - ГР 1+7

I — ,0 1 и радиусом —, а на концах этой полуокружности I р =-,

1 - 7 1

Р = I обладает свойством /2-устойчивости [5].

Запишем итоговый вид СЯОБ, для определенности выбрав Р =1 +7

2

У п+1 = У п Яе к; I Е - ^ х!и (уи, 0 1 к = { Г уи, I + +\.

В точках Р = ~~ и Р = СКОБ имеет второй порядок аппроксимации и обладает /2-устойчивостью. Но при попытке совместить эти качества схемы при переходе к действительным числам оказывается, что такой переход без потери какого-либо из свойств невозможен: /2-устойчивость не является достижимой в случае действительных коэффициентов, поэтому приходится делать выбор между вторым порядком аппроксимации и /1-устойчи-востью. Значения коэффициента р , при которых выполняется одно из описанных выше

свойств, равны соответственно Р = 1 (пересечение прямой Яе(Р) = -2 с осью действительных чисел) и Р = 1 (пересечение полуокружности

Р-2

1 1 = К-е(Р) > ^ с осью действительных чисел).

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

3. Четырехстадийный явный метод Рунге — Кутты. Этот метод обладает четвертым порядком точности. Переход к решению на новом слое происходит по следующим формулам [1]:

у и+1 = у п+тК;к = 6 (к1 + 2к 2+2к з+к 4),

где

к1 =Г (У п, гп);к 2 = г ( у п + \ к гп + £|;к з = г ( у п + 2 ^ гп + £|;

к4 = г(Уп + тk3, гп + т).

4. Метод конечных суперэлементов [11-14]. Пусть поставлена задача Коши (1) для линейной системы с постоянными коэффициентами Г (и) = Аи. Точное решение данной задачи имеет вид и(г) = и 0еАг,

0 < г < го [11].

Рассмотрим интервал Аг = , входящий во временной интервал всей задачи, 0 <Аг <г0. Пусть гк = к Аг. Тогда и(гк) = и0еАгк =

ААг ААг ААг ¡, \ ААг гр г

= и0 е е ...е 1 = и (гк-1)е . Таким образом, для нахождения знак

чения функции и(гк ) в конце данного интервала требуется знать значение функции и(гк-1) в конце предыдущего интервала и значение матричной экспоненты еААг. Для ее нахождения введем п (размерность исходной задачи) вспомогательных задач следующего вида:

ф; = АФт, гк-1 <г < гк;

Фт (0) = (0,0..,1,..,0)г,

т

где 1 < т < п. Эти вспомогательные задачи можно решать с малым шагом т « Аг стандартными методами, используемыми для решения систем ОДУ. Тогда решение поставленной задачи на слое гк можно

представить следующим образом: u(tk) = ^ u к-1, J&m, где uk-1, m — m-я

ш=1

компонента вектора Uk-1.

Данный метод называют методом конечных суперэлементов, так как в описанном виде он совпадает с МКСЭ Федоренко [11-14], называемым также методом матричной экспоненты [11].

Тестовые задачи. Рассмотрим примеры, использованные в данной работе для тестирования описанных методов [3]. Тест 1

Дана система

г

u1 =

u2 = (l0)u1 + (1 +V1 )-v1u3;

u3 = (o -V1) u1 + 2v1U2 + -V1)u3;

u4 = (( - I - V1)u1 + 2v1U2 + (| - V1 - I2)u3 + (I2 + V2 )u4 - V2u5; u5 = ( -^1 -V1 )u1 + 2V1u2 + ( -V1 2 -V 2 )u3 + 2v 2u4 + ( -V 2 )u5.

Пусть u2 (0) = u3 (0) и u4 (0) = u5 (0). Тогда точное решение данной системы линейных ОДУ (при 0 < t < 1) имеет следующий вид: u1(t) = u1(0)el0t;

u2 (t) = u1 (t) + (u2 (0) - u1 (0)) el1t cos V1t; < u3 (t) = u1 (t) + V2 (u2 (0) - u1 (0)) e11 sin v1 (t + n /4); u4 (t) = u3 (t) + (u4 (0) - u2 (0)) el2t cos V2t; u5 (t) = u3 (t) + V2 (u4 (0) - u2 (0)) el2t sin v2 (t + n /4). Собственные значения матрицы системы: X = 10; X23 =I1 ^z'v1;

X4,5 =12 ± iv2.

Рассмотрим пять случаев различных параметров исследуемой системы.

Случай 1: плохо обусловленная задача. Начальные условия: u1 (0) = 0,1; u2 (0) = u3 (0) = 1; u4 (0) = u5 (0) = 0,5, значения параметров:

|0 = 10; |1 = 4; v1 = 20n; |2 = 5; V2 = 100.

Определим число обусловленности M задачи с матрицей A сле-

дующим образом: M =

. -1

A [15]. Тогда число обусловленности

матрицы рассматриваемой системы имеет значение М = 241 (здесь и далее оценка проведена с использованием октаэдрической нормы [1]).

Собственные значения задачи: ^ = 10; Л23 = 4 ± 20л/; X 45 = 5 ± 100/.

Все их действительные части являются положительными, а значит, задача не является жесткой в смысле приведенного определения.

Случай 2: хорошо обусловленная задача. Начальные условия: и1(0) = 1; и2(0) = и3(0) = 1,5; и4(0) = и5(0) = 2,5. Значения параметров:

10 = -2; I =1; V =1; 12 = -1; ^ =10.

Число обусловленности матрицы М = 105. Собственные числа матрицы: Х1 =-2; Х23 = 1 ±/; X45 =-1 ± 10/. Не все их действительные части отрицательны, поэтому данная задача не является жесткой.

Случай 3: быстроосциллирующая задача. Начальные условия: и1(0) = 0,5; и2(0) = и3(0) = 0,8; и4(0) = и5(0) = 2. Значения параметров:

10 = -2; | = 1; V! = 1; ^ = -1; у2 = 1000.

Число обусловленности задачи М = 10499. Собственные числа: Х1 = -2; X 2 3 = 1 ± /; X 4 5 = -1 ± 1000/. Данный случай аналогичен предыдущим, задача не является жесткой.

Случай 4: жесткая задача. Начальные условия: и1(0) = 10;

и2 (0) = и3 (0) = 11; и4(0) = и5(0) = 111, значения параметров: |0 =-100; I = -1; V = 1; |2 = -10000; v2 = 10.

Число обусловленности матрицы М = 79956. Собственные числа: Х1 =-100; X23 =-1 ±/; X45 =-10000 ± 10/. Число жесткости задачи

£ = 104.

Случай 5: жесткоосциллирующая задача. Начальные условия: и1 (0) = 100; и2 (0) = и3(0) = 101; и4 (0) = и5 (0) = 201. Значения параметров: |0 = -10000; | = 1; v1 = 1; |2 = -100; v2 = 1000.

Число обусловленности задачи М = 175084. Собственные числа: ^ =-10000; X2,3 = 1 ±/; X4,5 =-100 ± 1000/. Тест 2

Рассмотрим следующую линейную систему ОДУ:

г

и - ^^ ;

и2 = л+№;

г

и3 --;

<

Г

и 4 — и3 + 12^; и5' = 2и4 + |2и5; и6' = 3и5 + |2и6.

Точное решение этой системы имеет вид (0 < 7 < 1): 'щЦ) = м1(0)е^1/;

и2(7) = ((0) + «1(0)7 );

«3( 7) = и3 (0)вЦ2/;

' «4(0 = («4(0) + и3(0)7

и5 (7) = («5 (0) + 2«4 (0) 7 + «3 (0) 72 ) в»2';

«6 (7) = («6 (0) + 3«5(0> + 3«4 (0) 72 + «3 (0) 73 .

При начальных условиях «1 (0) = «2 (0) = 1; «3 (0) = «4 (0) = = «5(0) = и6(0) = 1000 и значениях параметров =-1, »2 =-10 000

данная система линейных ОДУ является жесткой. Число обусловленности решаемой задачи М = 20 006, жесткость £ = 10 000. Тест 3

Дано: « = -а«.

Точное решение: и ( г) = и (0)в-а.

Начальное условие решаемой задачи (0 < I < 1): «(0) = 1, параметр а принимает четыре значения: а = 1; 10; 100; 1000 , от которых зависит сложность решаемой задачи. Тест 4

Рассмотрим систему

«1 = -а«1;

<

г

«2 — «2.

Точное решение задачи:

«1(7) = «1(0)в-а;

«2(0 = «2(0)в-/.

Начальные условия (0 < 7 < 1): «1(0) = «2(0) = 1.Параметр а, регулирующий жесткость задачи, принимает четыре значения: а = 1; 10; 100; 1000.

Зависимость чисел обусловленности и жесткости от параметра а выглядит следующим образом: М = а, £ = а. Тест 5

Рассмотрим систему

«1 = -а«2;

<

г

«^ — а.«1 «ъ.

Начальные условия (0 < г < 1): м1(0) = и2(0) = 1. Точное решение задачи:

и1 (г) = е 2

( Ьг Л

(1 - 2 а)81п- ы

-— + соб—

Ь 2

и2(Г) = е 2

У

^ • Ьг л

(2а-1)81пт ы

-— + соб—

Ь2

У

где Ь = л/4а2 -1.

Параметр задачи а принимает четыре значения: а = 1; 10; 100; 1000. Чем выше значения параметра а, тем сильнее осциллирует решение задачи. Случай а = 1000 соответствует жесткоосциллирующей задаче.

Результаты тестовых расчетов. Ниже для каждого из рассмотренных методов приведена абсолютная ошибка численного решения в зависимости от используемого при вычислениях временного шага. Абсолютная ошибка А вычислялась по следующей формуле: А = тах || и - уг ||м, где иг, уг — векторы точного и численного ре-

г

шений на г-м временном слое, ||. ||м — кубическая норма вектора [1].

1. (4,2)-метод. В табл. 1-5 приведена абсолютная погрешность А (4,2)-метода, полученная при использовании различных временных шагов т.

Таблица 1

г

г

Тест 1 Вариант 1 Вариант 2 Вариант 3 Вариант 4 Вариант 5

т = 1,00 • 10-5 2,78 • 10-10 9,29 •Ю-15 -7 1,71 •Ю ' 8,64 • 10-5 8,64 • 10-5

т = 4,00 • 10-5 6,94 •Ю-8 -14 1,22 • 10 14 4,35 •Ю- 5 2 1,48 •Ю 2 2 1,48 •Ю 2

т = 1,60 • 10-4 1,78 -10-5 9,42 •Ю-13 1,09 • 10-2 1,32 1,32

т = 6,40 • 10-4 4,54 • 10-3 2,39 •Ю-10 9,44 • 10- 1 9,84 1,01 • 101

т = 2,56 •Ю-3 1,11 6,09 •Ю-8 1,66 6,39 4,38 •Ю1

Из результатов для теста 1, представленных в табл. 1, видно, что (4,2)-метод на малых шагах (до т = 4,00 -10-5) адекватно работает для всех типов задач. При увеличении шага интегрирования в случаях жесткой и жесткоосциллирующей задач абсолютная ошибка становится велика, и численное решение начинает сильно отличаться от точного.

Таблица 2

Тест 2 т = 2,00 -10 5 т = 8,00 -10 5 _4 т = 3,20 -10 4 _3 т = 1,28 -10 3 3 т = 5,12 -10 3

А 1,20 -10 2 1,57 5,39 -101 9,39 -101 3,71 -101

Задача, поставленная в тесте 2, является жесткой. При ее решении с использованием малых шагов интегрирования численное решение не сильно отличается от точного (см. табл. 2), но при увеличении шага (от т = 8,00 -10" 5) наблюдается рост абсолютной ошибки до довольно большого значения ( ~ 10 ).

Таблица 3

Тест 3 а = 1 а = 10 а = 100 а = 1000

т = 1,00 -10_3 1,08 -10_ 14 9,87 -10_ 11 8,64 -10_4 3,34 -10_3

т = 1,00 -10"1 8,64 -10_7 3,34 -10_3 1,01 -10_ 1 2,05 -10_2

Таблица 4

Тест 4 а = 1 а = 10 а = 100 а = 1000

т = 1,00 -10_3 1,08 -10_ 14 9,87 -10_ 11 8,64 -10_4 3,34 -10_3

т = 1,00 -10"1 8,64 -10_7 3,34 -10_3 1,01 -10_ 1 2,05 -10_2

Таблица 5

Тест 5 а = 1 а = 10 а = 100 а = 1000

т = 1,00 -10_3 1,43 -10_ 14 2,28 -10_9 2,31 -10_4 1,24

т = 1,00 -10"1 1,48 -10_ 6 1,16 -101 1,15 1,28

При оценке результатов решения тестовых задач 3-5 (см. табл. 3-5) заметим, что (4,2)-метод успешно решает задачи со значениями параметра а = 1; 10; 100, при решении задач с параметром а = 1000 в тесте 5 наблюдается резкое увеличение численной ошибки метода.

2. СЯ08. В табл. 6-10 отображена зависимость абсолютной ошибки численного решения А метода СЯОБ от выбранного шага интегрирования т.

Таблица 6

Тест 1 Вариант 1 Вариант 2 Вариант 3 Вариант 4 Вариант 5

т = 1,00 -10_5 3 1,54 -10 3 _9 8,60-10 9 1,04 -10 2 5,69-10 2 5,69-10 2

т = 4,00 -10_5 2 2,47 -10 2 1,38 -10 ' 1,65 -101 7,26 -10_Х 7,29 -10_ 1

т = 1,60 -10 4 3,95 -10_Х 2,20 -10_6 1,57 5,58 5,60

т = 6,40 -10_4 6,32 3,51 -10_5 1,96 3,43 2,73 -101

т = 2,56 -10_3 8,73101 5,62 -10_4 1,67 3,0110_ 1 5,17 -101

По результатам, представленным в табл. 6, видно, что применение метода СЯОБ для решения плохо обусловленной задачи (вариант 1), жесткой (вариант 4) и жесткоосцилирующей (вариант 5) задач на шагах интегрирования от т = 1,60 -10—4 приводит к довольно большой абсолютной ошибке.

Таблица 7

Тест 2 т = 2,00 -10 5 т = 8,00 -10 5 -4 т = 3,20-10 4 -3 т = 1,28-10 3 3 т = 5,12 -10 3

А 2,12 2,24101 6,65 -101 1,04 -101 7,34 -10-1

При решении задачи с большой жесткостью (тест 2) даже на малых шагах интегрирования т применение СЯОБ приводит к большой ошибке численного решения (см. табл. 7).

Таблица 8

Тест 3 а = 1 а = 10 а = 100 а = 1000

т = 1,00 -10-3 6,13 -10-8 6,09 -10—6 5,69 -10—4 3,21 -10—2

т = 1,00 -10-1 5,66 -10—4 3,21 -10—2 1,63 -10—2 1,96 -10—4

Таблица 9

Тест 4 а = 1 а = 10 а = 100 а = 1000

т = 1,00 -10-3 6,13 -10-8 6,09 -10—6 5,69 -10—4 3,21 -10—2

т = 1,00 -10-1 5,69 -10 4 3,21 -10 2 2 1,63-10 2 4 5,69 -10 4

Таблица 10

Тест 5 а = 1 а = 10 а = 100 а = 1000

т = 1,00 -10-3 —7 1,10-10 7 4 1,39 -10 4 1,41 -10—1 1,46

т = 1,00 -10-1 —3 1,03-10 3 7,01 -10—1 1,30 1,30

При решении тестовых задач 3 и 4 даже на крупных шагах интегрирования т и при больших значениях параметра ошибка решения, полученного с помощью одностадийной комплексной схемы СЯОБ, не является значительной. Аналогичные результаты получаются и при решении тестовой задачи 5 со значениями параметра а = 1 и а = 10. В случаях а = 100 и а = 1000 возникает резкое увеличение абсолютной ошибки.

З.Четырехстадийный метод Рунге — Кутты. В табл. 11-15 приведены абсолютные ошибки А метода Рунге — Кутты в зависимости от используемого шага интегрирования т.

Таблица 11

Тест 1 Вариант 1 Вариант 2 Вариант 3 Вариант 4 Вариант 5

т = 1,00 -10—5 6,15 -10—16 1,52 -10—18 2,29 -10—10 4,70 -10—5 7,45 -10—5

т = 4,00 -10—5 —13 5,09 -10 13 8,80 -10—18 —7 2,23 -10 ' 2 1,52 -10 2 2 2,41 -10 2

Окончание табл. 11

Тест 1 Вариант 1 Вариант 2 Вариант 3 Вариант 4 Вариант 5

х = 1,60 • 10-4 6,73 -10- 10 5,98 -10- 15 1,55 -10-4 9,67 1,53 -101

-4 х = 6,40 -10 4 1,63 -10 7 - 12 1,50 -10 12 2 4,04 -10 2 44 1,37 -10 44 2,31 -10

х = 2,56 -10 3 3 4,65 • 10 3 9,27 • 10 9 1,92 1,9110567 4,73 -10567

Из результатов применения четырехстадийного метода Рунге — Кутты видно, что он является непригодным для решения быстроос-циллирующих задач (вариант 3), жестких задач (вариант 4) и жестко-осциллирующих задач (вариант 5).

Таблица 12

Тест 2 х = 2,00 -10 5 х = 8,00 -10 5 -4 х = 3,20 -10 4 х = 1,28 -10 3 3 х = 5,12 -10 3

А 2 1,16 -10 2 4,81 17 1,47 -10 2,03 -1041 7,16 -101060

Тест 2 является жесткой задачей, о чем свидетельствуют и результаты табл. 12. Можем наблюдать резкий рост абсолютной ошиб-

ки численного решения.

Таблица 13

Тест 3 а = 1 а = 10 а = 100 а = 1000

х = 1,00 -10-3 3,09 -10- 15 3,09 -10- 11 3,33 -10-7 7,12 -10-3

х = 1,00 -10-1 3,32 -10-7 7,12 -10-3 1,50 -1022 2,65 -1059

Таблица 14

Тест 4 а = 1 а = 10 а = 100 а = 1000

х = 1,00 -10-3 4,34 -10- 15 3,09 -10- 11 -7 3,33-10 7 3 7,12-10 3

х = 1,00 -10-1 4,69 -10 ' 3 7,12-10 3 22 1,50 -10 59 2,65 -10

Таблица 15

Тест 5 а = 1 а = 10 а = 100 а = 1000

х = 1,00 -10-3 4,56 -10- 15 6,98 -10- 10 7,13 -10-5 1,25

х = 1,00 -10-1 4,20-10 ' 2 6,77 -10 2 23 3,61 -10 59 5,33-10

Использование метода Рунге — Кутты четвертого порядка точности для решения тестовых задач 3-5 со значениями параметра а = 1 и а = 10 приводит к довольно малым ошибкам численного решения. При увеличении сложности решаемой задачи (а = 100, а = 1000) наблюдаем рост абсолютной ошибки численного решения на шагах, больших х = 1,00 -10-3.

4. МКСЭ. В табл. 16-20 приведены абсолютные ошибки МКСЭ в зависимости от выбранного шага интегрирования т, размер используемой ячейки АГ = 1-10 .

Решение вспомогательных задач получено с помощью метода Рунге — Кутты четвертого порядка точности и (4,2)-метода. Результаты, приведенные в таблицах первыми, получены с использованием (4,2)-метода, а результаты, приведенные вторыми, — с помощью метода Рунге — Кутты четвертого порядка точности (тоном выделены меньшие ошибки).

Таблица 16

Тест 1 АГ т Вариант 1 Вариант 2 Вариант 3 Вариант 4 Вариант 5

т = 1,00 -10_6 1,00 -104 2,74 -10 12 8,95 -10 15 1,76 -10_ 11 9 9,87-10 9 9 9,87-10 9

9,4510_12 5,4610_ 15 5,2110_12 9 3,0910 9 9 3,0910 9

т = 4,00 -10_ 6 2,50 -103 3,34 -10_П 9,40 -10_15 9 4,36 -10 9 2,41 -10 6 2,41 -10_6

1,0810_11 5,3210_15 9 1,3310 9 8,1110 7 8,1110_7

т = 8,00 -10_6 3 1,25 -103 1,25 -10_ 10 3,9010_11 9,33 -10_15 5,1010_15 6,97 -10_8 2,1310_8 3,64-10 5 1,3410_ 5 3,64 -10_5 1,3410_5

т = 2,00 -10_ 5 5,00 -102 4,33 -10 9 1,33 -10 9 9,41 -10_ 15 5,12 -10_15 2,72 -10_6 8,32 -10 ' 1,20 -10 3 4 5,80 -10 4 3 1,20 -10 3 4 5,80 -10 4

т = 8,00 -10_ 5 1,25 -102 1,11 -10_ 6 3,39 -10_7 6,66 -10_14 2,28 -10_ 14 6,94 -10_4 2,13 -10_4 1,57 -10^ 2,40 -10_1 1,57 -10^ 2,40 -10_1

Таблица 17

Тест 2 т = 2,00 -10 6 т = 5,00 -10 6 т = 1,00 -10 5 т = 4,00 -10 5 т = 8,00 -10 5

АГ т 5,00 -103 2,00 -103 3 1,00 -103 2,50 -102 2 1,25 -102

А 1,56 -10_6 4,9910_7 5,81 -10_5 2,0010_ 5 4 8,64-10 4 4 3,33-10 4 1,48 -10_ 1 1,08 -101 1,57 2,40

Таблица 18

Тест 3 АГ т а = 1 а = 10 а = 100 а = 1000

т = 1,00 -10_6 1,00 -104 8,65 -10_16 5,4710_ 16 1,18 -10_15 9,4210_ 16 1,18 -10_15 9,5010_16 1,08 -10_ 14 3,5610_15

т = 1,00 -10_5 3 1,00 -103 7,74 -10_16 1,0310_15 1,17 -10_ 15 9,5010_16 14 1,08 -10 14 3,5610_15 9,87 - 10_П 3,0910_ 11

Окончание табл. 18

Тест 3 At т а = 1 а = 10 а = 100 а = 1000

т = 1,00 • 10"4 2 1,00 • 102 8,22 -10" 16 5,12 •Ю" 16 " 14 1,08 •Ю 14 3,56 • 10"15 9,87 •Ю" 11 3,09 •Ю" 11 "7 8,64 • 10 7 3,32 • 10 '

т = 1,00 • 10"3 1,00 • 101 14 1,08 •Ю 14 3,56 • 10"15 9,87 •Ю" 11 3,09 •Ю" 11 "7 8,64 • 10 7 3,33 • 10 7 3 3,34 •Ю 3 3 7,12 •Ю 3

Таблица 19

Тест 4 At т а = 1 а = 10 а = 100 а = 1000

т = 1,00 •Ю-6 1,00 • 104 8,65 • 10"16 5,4710"16 1,17 •Ю"15 9,4210"16 1,17 •Ю"15 1,0710" 15 1,08 • 10"14 3,5610"15

т = 1,00 -10~5 1,00 • 103 7,74 • 10"16 1,0310" 15 1,17 •10"15 1,0310" 15 1,08 •Ю" 14 3,5610"15 9,87 •Ю" 11 3,09-Ю"11

т = 1,00 -10~4 2 1,00 • 102 8,22 •Ю" 16 5,01 • 10"16 1,08 •Ю" 14 3,56 • 10"15 9,87 •Ю" 11 3,09 •Ю" 11 8,64 • 10"7 7 3,33 • 10 7

т = 1,00 40"3 1,00 • 101 14 1,08 •Ю 14 3,56 • 10"15 9,87 •Ю" 11 3,09 •Ю" 11 7 8,64 • 10 7 3,32 • 10 ' 3 3,34 •Ю 3 3 7,12 •Ю 3

Таблица 20

Тест 5 At т а = 1 а = 10 а = 100 а = 1000

т = 1,00 • 10"6 1,00 • 104 1,44 •Ю"15 5,7010"16 1,19 • 10"14 5,0610"15 7,46 •Ю"14 5,07-10"15 2,41 •Ю" 11 7,1510"12

т = 1,00 • 10"5 3 1,00 • 103 1,00 •Ю" 15 1,2910" 15 6,93 •Ю"15 1,0510"14 12 2,31 •Ю 12 7,9810"13 7 2,34 • 10 7,1510"7

т = 1,00 • 10"4 2 1,00 • 102 1,08 • 10"15 5,09 •Ю" 16 13 2,36 •Ю 13 14 7,02 •Ю 14 2,33 •Ю"8 7,12 •Ю 9 3 2,32 • 10 3 4 7,15•10 4

т = 1,00 • 10"3 1,00 • 101 1,43 • 10"14 4,92 •Ю"15 2,28 •10"9 6,97 •Ю"10 2,31 •Ю"4 7,13 •Ю" 5 1.24 1.25

В табл. 21-24 показана зависимость абсолютной ошибки МКСЭ при фиксированном (малом) шаге интегрирования т от выбранного размера ячейки. Рассмотрены жесткие случаи (тест 2, а также тесты 3-5 при а = 1000). В таблицах через At обозначен размер ячейки, через т — шаг суммирования метода.

В приведенных результатах для всех тестов отсутствуют результаты для т = 5 -10_6 и АГ = 1 -10_6 из-за нецелесообразности использования МКСЭ при столь близких значениях размеров ячейки и шага суммирования.

Таблица 21

Тест 2 АГ = 1 -10 1 _2 АГ = 1 -10 _3 АГ = 1 -10 _4 АГ = 1 -10 АГ = 1 -10 5 АГ = 1 -10 6

_7 т = 5 -10 ' 1,92 -10 9 1,92 -10 9 1,92 -10 9 1,92 -10 9 1,92 -10 9 9 1,92-10 9

т = 1 -10_6 3,09 -10_8 3,09 -10_8 3,09 -10_ 8 3,09 -10_8 3,09 -10_ 8 3,09 -10_ 8

т = 5 -10_6 2,00 -10_ 5 2,00 -10_ 5 2,00 -10_5 2,00 -10_ 5 2,00 -10_5 -

Таблица 22

Тест 3 АГ = 1 - 10_ 2 АГ = 1 -10 _3 АГ = 1 -10 _4 АГ = 1 -10 АГ = 1 -10_5 АГ = 1 -10_6

т = 5 -10_7 8,23 -10_16 8,23 -10_16 8,23 -10_16 8,23 -10_16 8,23 -10_ 16 8,05 -10_ 16

т = 1 -10_6 3,56 -10_15 3,56 -10_15 3,56 -10_15 3,56 -10_15 3,56 -10_ 15 3,56 -10_ 15

т = 5 -10_6 1,92 -10_ 12 1,92 -10_12 1,92 -10_ 12 1,92 -10_ 12 1,92 -10_12 -

Таблица 23

Тест 4 АГ = 1 - 10_ 2 АГ = 1 -10 _3 АГ = 1 -10 _4 АГ = 1 -10 АГ = 1 -10_5 АГ = 1 -10 6

т = 5 -10 ' 7,17 -10_15 8,23 -10_16 8,23 -10_16 6,78 -10_15 5,72 -10_15 _14 1,89 -10 14

т = 1 -10_6 4,60 -10_15 3,56 -10_15 3,56 -10_15 3,56 -10_15 3,99 -10_ 15 8,03 -10_15

т = 5 -10_6 12 1,92 -10 12 12 1,92 -10 12 12 1,92 -10 12 12 1,92 -10 12 12 1,92 -10 12 -

Таблица 24

Тест 5 АГ = 1 -10_ _2 АГ = 1 -10 2 АГ = 1 -10 3 АГ = 1 -10_4 АГ = 1 -10_5 АГ = 1 -10_6

т = 5 -10_7 1,59 -10_ 11 4,47 -10_ 13 4,47 -10_13 1,59 -10_П 1,59 -10_ 11 1,59 -10_П

т = 1 -10_6 1,66 -10_ 11 7,15 -10_12 7,15 -10_ 12 1,66 -10_ 11 1,66 -10_11 1,66 -1011

т = 5 -10_6 9 4,47 -10 9 4,47 -10 9 4,47 -10 9 4,47 -10 9 4,47 -10 9 -

Табл. 21-24 демонстрируют очень слабую зависимость результатов расчетов от используемого размера ячейки АГ.

Анализ полученных результатов. Для сравнения результатов, полученных с помощью (4,2)-метода (см. табл. 1-5) и метода СЯОБ (см. табл. 6-10), рассмотрим два случая. Первый — решение задач с малым числом жесткости. В этом случае лучший результат показывает (4,2)-метод, поскольку обладает четвертым порядком точности, в то время как СЯОБ — только вторым.

Второй случай — решение жестких задач. Если допустимо использование малых шагов интегрирования, то предпочтительнее выбрать (4,2)-метод. При необходимости использования более крупных

шагов следует применять комплексную схему СЯОБ. Она является ¿2-устойчивой схемой, в то время как (4,2)-метод ¿-устойчив.

Для оценки возможности применения явного метода Рунге — Кутты четвертого порядка точности рассмотрим результаты, приведенные в табл. 11-15. При решении мягких задач (тест 1, вариант 2; тесты 3-5 при значениях параметра а = 1, а = 10) с использованием больших т (относительно рассматриваемого диапазона шагов интегрирования) ошибка численного решения довольно мала. Но при применении метода Рунге — Кутты к жестким или осциллирующим задачам (тест 1, варианты 1, 3, 4, 5; тест 2; тесты 3-5 при а = 100, а = 1000) наблюдается резкое возрастание ошибки численного решения даже при использовании средних шагов интегрирования т из рассматриваемого диапазона. Поэтому явный четырехстадийный метод Рунге — Кутты не является приемлемым для решения жестких задач.

Далее обсудим выбор метода, используемого в качестве вспомогательного для МКСЭ. В данной работе выбор совершался между явным методом Рунге — Кутты четвертого порядка точности и (4,2)-методом. Результаты расчетов представлены в табл. 16-20. Метод Рунге — Кут-ты на малых шагах применительно к задачам умеренной жесткости показывает лучшие результаты, чем (4,2)-метод, но его применение требует больше машинного времени. Это связано с тем, что при реализации метода Рунге — Кутты четвертого порядка для получения численного значения искомой функции на новом временном слое требуется четыре вычисления правой части решаемой задачи, а при реализации (4,2)-метода — всего два.

Однако при увеличении числа жесткости и шага интегрирования задачи результаты (4,2)-метода приближаются к результатам метода Рунге — Кутты, и разница между ними становится незначительной (тест 1, варианты 4, 5; тест 2; тесты 3-5 при а = 1000 ). На задачах с высоким уровнем жесткости метод Рунге — Кутты будет вести себя некорректно, в то время как (4,2)-метод благодаря ¿-устойчивости будет давать удовлетворительные результаты. Исходя из этих рассуждений, (4,2)-метод выбран в качестве вспомогательного метода для МКСЭ.

Проведение расчетов, результаты которых представлены в табл. 2124, предполагало определение оптимального соотношения между шагом интегрирования и размером ячейки МКСЭ. Однако при таких малых шагах точность вычислений напрямую зависит от размера шага интегрирования, поэтому зависимость ошибки метода от размера ячейки в исследуемом диапазоне не была установлена.

Применение МКСЭ для нелинейных задач. МКСЭ для линейных задач представляет собой довольно простой в реализации и лучший по получаемым результатам алгоритм. Попробуем обобщить МКСЭ на нелинейные задачи.

(2)

Для применения МКСЭ проведем линеаризацию в окрестности вектора и0:

где А — матрица Якоби поставленной задачи при и = и0, Ц =

Тогда вместо исходной задачи требуется решить вспомогательную:

ходим уже сам вектор решения задачи (3). Для ответа на вопрос, соответствует ли решение вспомогательной задачи (3) решению поставленной нелинейной задачи (2), построим алгоритм контроля точности МКСЭ для нелинейных задач.

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

где у — найденный вектор приближенного решения задачи (3) в текущий момент времени; в — некая наперед заданная точность вычислений; 80 — малая постоянная, 80 > 0. При выполнении условия (4) происходит переход на следующий временной шаг. Если же условие (4) не выполнено, то происходит очередная линеаризация задачи (2) и переход к новой временной ячейке At.

Тестирование нелинейного МКСЭ. Рассмотрим тестовый пример, иллюстрирующий работу МКСЭ применительно к нелинейным системам ОДУ.

Тест 6

Представлена система

Ц(и) = А(и " и0) + Ц(и0) = Аи + Ц

= Г(и0)" Аи0.

(4)

f 2 щ ^) = ащ ^)и2 ^);

f 2 < щ2 ^) = " )щ2 (t);

^(0) = щ2(0) = 1.

Точное решение данного тестового примера имеет следующий вид:

щЦ) = , ) = е-а.

Жесткость поставленной задачи зависит от выбора параметра а. Поэтому рассмотрим относительные ошибки численного решения, полученного МКСЭ для нелинейных задач, в зависимости от параметра а, шага метода т, а также выбранного размера ячейки At, приведенных в табл. 25-28).

В качестве вспомогательного метода для МКСЭ использован (4,2)-метод.

_3

Заданная точность линеаризации в = 10 , время исследования ^ = 1.

В каждой ячейке таблицы в первой строке указана относительная ошибка МКСЭ, во второй — средний размер реально используемой методом ячейки.

Таблица 25

а = 0,1 т = 110 1 _2 т = 110 2 _3 т = 110 3 _4 т = 110 4 т = 110 5

At = 110_ 1 3,32-Ю-6 1,00-10-1 3,65-Ю-6 1,00-10-1 3,33-Ю-6 1,00-10-1 3,66-Ю-6 1,00-10-1 3,33-10-6 1,00-10-1

At = 110 2 - 3,37-10-8 1,00-10-2 3,33-10-8 1,00-10-2 3,37-10-8 1,00-10-2 3,33-10-8 1,00-10-2

At = 110 3 - - 3,33-10-10 1,00-10-3 3,34-10-10 1,00-10-3 3,33-10-10 1,00-10-3

По полученным результатам, представленным в табл. 25, можно сделать вывод о том, что при решении мягких задач ошибка МКСЭ обусловлена выбранным размером ячейки; повторная линеаризация внутри самой ячейки не происходит.

Таблица 26

а = 1 т = 110 1 _2 т = 110 2 т = 110 3 _4 т = 110 4 т = 110 5

At = 110_1 2,92-10-3 1,00-10-1 _4 5,22-10 4,00-10-2 _4 3,38-10 3,25-10-2 -4 3,33-ю-4 3,20-10-2 -4 3,32-Ю-4 3,20-10-2

2 At = 110 2 - 3,32-10-5 1,00-10-2 3,28-10-5 1,00-10-2 3,32-10-5 1,00-10-2 3,28-10-5 1,00-10-2

3 At = 110 3 - - 3,33-10-7 1,00-10-3 3,33-10-7 1,00-10-3 3,33-10-7 1,00-10-3

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

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

Таблица 27

а = 10 т = 110 1 _2 т = 110 2 _3 т = 110 3 _4 т = 110 4 т = 110 5

Дt = 1-10"1 7 2,40-10' 1,00-10-1 1,01-10-2 1,00-10-2 3,58-10-3 4,00-10-3 2,67-10-3 3,30-10-3 2,55-10-3 3,22-10-3

Д = 110 2 - 1,01-10-2 1,00-10-2 3,59-10-3 4,00-10-3 2,65-10-3 3,30-10-3 2,55-10-3 3,22-10-3

_3 Д = 110 3 - - _4 3,06-10 1,00-10-3 3,06-10-3 1,00-10-3 _4 3,06-10 1,00-10-3

По данным табл. 27 можно сделать вывод о том, что при решении задач средней жесткости с помощью МКСЭ с использованием малого размера ячейки погрешность метода обусловлена размером ячейки, а выбранный шаг имеет лишь косвенное значение. И наоборот, при больших размерах ячейки погрешность метода напрямую зависит от используемого шага.

Таблица 28

а = 50 т = 110 1 _2 т = 110 2 _3 т = 110 3 _4 т = 110 4 т = 110 5

Д = 110_1 1,00 1,00-10-1 9,99-10-1 1,00-10-2 3,67-10-2 4,00-10-3 6,87-10-3 _4 7,00-10 4,24-10-3 _4 6,50-10

Д = 110_2 - 9,99-10-1 1,00-10-2 3,67-10-2 1,00-10-3 6,87-10-3 _4 7,00-10 4,24-10-3 _4 6,50-10

Д = 110_3 - - 3,67-10-2 1,00-10-3 6,87-10-3 _4 7,00-10 4,24-10-3 _4 6,50-10

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

Решение дифференциального уравнения Ван-дер-Поля. Одним из классических примеров нелинейных дифференциальных уравнений второго порядка, относящихся к жесткому типу при большом значении некоторого параметра д, является уравнение Ван-дер-Поля [16]. Это уравнение играет важную роль в прикладных задачах, так как к нему сводятся дифференциальные уравнения, описывающие динамику развития колебаний в различных колебательных системах, например автогенераторах на электронных лампах, биполярных и полевых транзисторах. Рассмотрим результаты, которые дает МКСЭ применительно к этой задаче.

Уравнение Ван-дер-Поля выглядит следующим образом:

и" + |(и2 _1)м' + ы = 0. Представим его в виде системы ОДУ:

ых = ы2;

ы^ = |(1 _ ы^) ы2 _ ых.

(5)

(6)

При исследовании численного решения задачи (6) выясняется, что период решения поставленной задачи возрастает с ростом д. Это усложняет численный анализ уравнения Ван-дер-Поля. Поэтому

масштабируем решение, используя замену 1 = —; zl(t') = ы^);

г2(/Г) = |ы2(1) [16]. В уравнениях для удобства заменим 2 на ы, 1 на 1 и получим конечный вид задачи Ван-дер-Поля:

ы{ = ы2;

ы2 =|((1 _ ы!2)ы2 _ ы1); ы1(0)=ыю; ы2(0) = ы20.

(7)

Время исследования работы системы (7) — = 20.

Жесткость задачи (7) увеличивается с ростом параметра Для ее численного анализа выберем специализированные алгоритмы: МКСЭ для нелинейных систем ОДУ и (4,2)-метод.

Рассмотрим результаты (представленные графически далее) работы алгоритмов при различных значениях |2, размерах ячейки (в случае использования МКСЭ) А— и шага т методов, а также оценим реальный размер ячейки 5—, используемый программой, реализующей МКСЭ.

Для сравнения получаемых двумя методами численных результатов введем ошибку 5 МКСЭ относительно (4,2)-метода следующим

образом: 5 =

тах

з

тах N

где л

вектор численного решения,

полученного (4,2)-методом на г-м шаге алгоритма; у^ — вектор численного решения, полученный на г-м шаге с помощью МКСЭ.

Из физического смысла уравнения Ван-дер-Поля очевидно, что его решение имеет периодический характер. Поэтому введем еще

тах

з

один тип исследуемой ошибки: 5 =-й—г—, где у г — вектор

тах у г

г

численного решения на первом периоде колебаний (выбран из-за наименьшего уровня накопления численной ошибки); у ^ — вектор

численного решения на всех последующих периодах.

Ниже приведена серия задач с различными значениями ц2, описаны входные параметры и приведены графические представления решения задачи (7) (фазовые траектории численного решения). Константа контроля точности линеаризации 8 = 10_3 .

Результаты решения ряда задач приведены в табл. 29.

Таблица 29

Задача Ц2 Аг т Реально используемый МКСЭ размер ячейки 5г Ошибка МКСЭ относительно 4,2-метода 5 Ошибка периодичности 5 для МКСЭ Ошибка периодичности 5 для 4,2-метода

1 100 10 2 10 3 2,4-10-3 1,02-10-3 2,67-10-1 1,65-10-1

2 1000 10 3 10-4 _4 6,87-10 6,71-10-1 4,66-10-3 3,36-10-1

3 10000 10 3 5-10-5 _4 2,39-10 1,12 1,02 1,04

4 15000 10 3 5-10-5 _4 1,95-10 1,24 1,03 1,00

5 17500 10 3 5-10-5 _4 1,75-10 1,08 1,03 1,10

6 20000 10 3 5-10-5 _4 6,24-10 1,52-Ю1 2,67-Ю1 1,00

Решения задачи 1 приведены на рис. 1, задачи 6 — на рис. 2.

Рис. 1. Решение задачи 1 с помощью МКСЭ (а) и (4,2)-метода (б)

Рис. 2. Решение задачи 6 с помощью МКСЭ (а) и (4,2)-метода (б)

На основании полученных при численном решении результатов можно сделать следующий вывод: МКСЭ (в его использованном варианте) для нелинейных задач работает корректно лишь при решении задач умеренной жесткости. Это объясняется тем, что при каждой очередной итерации в МКСЭ при выходе ошибки предыдущей линеаризации за установленный предел происходит очередная линеаризация поставленной нелинейной задачи. Но при большом числе жесткости системы ОДУ его решение имеет резко изменяющийся характер. Поэтому даже при выполнении линеаризации на каждой итерации МКСЭ ее погрешность превышает допустимую, что влечет за собой и ошибку численного решения. Следовательно, при решении нелинейных задач с большой жесткостью желательно применять (4,2)-метод.

Заключение. При решении систем ОДУ с большим числом жесткости лучшим выбором является СЯОБ.

МКСЭ является точным методом решения линейных систем ОДУ в том смысле, что если решение вспомогательных задач метода является точным, то и итоговое решение задачи точно. При нахождении соответствующих решений численно в качестве вспомогательного метода следует применять (4,2)-метод из-за его высокой точности.

Предложенная реализация МКСЭ для нелинейных задач корректно работает на задачах умеренной жесткости, для решения задач высокой жесткости необходимо применять альтернативные методы.

Ряд дополнительных подробностей о выполненных тестовых расчетах представлен в [17].

ЛИТЕРАТУРА

[1] Галанин М.П., Савенков Е.Б. Методы численного анализа математических моделей. Москва, Изд-во МГТУ им. Н.Э. Баумана, 2010, 591 с.

[2] Калиткин Н.Н. Численные методы. Москва, Наука, 1978, 512 с.

[3] Арушанян О.Б., Залеткин С.Ф., Калиткин Н.Н. Тесты для вычислительного практикума по обыкновенным дифференциальным уравнениям. Вычислительные методы и программирование, 2002, т. 3, с. 11-19.

[4] Новиков Е.А. Исследование (т,2)-методов решения жестких систем. Вычислительные технологии, 2007, т. 12, № 5, с. 103-115.

[5] Калиткин Н.Н. Полуявные схемы для задач большой жесткости. ЭНТП, сер. Б, т. VII-1, ч. 1. Москва, Янус-К, 2008, с. 153-171.

[6] Альшин А.Б, Альшина Е.А, Калиткин Н.Н, Корягина А.Б. Схемы Розен-брока с комплексными коэффициентами для жестких и дифференциально-алгебраических систем. Журнал вычислительной математики и математической физики, 2006, т. 46, № 8, с. 1392-1414.

[7] Альшин А.Б., Альшина Е.А., Лимонов А.Г. Двустадийные комплексные схемы Розенброка для жестких систем. Журнал вычислительной математики и математической физики, 2009, т. 49, № 2, с. 270-287.

[8] Ширков П. Д. Оптимально затухающие схемы с комплексными коэффициентами для жестких систем ОДУ. Математическое моделирование, 1992, т. 4, № 8, с. 47-57.

[9] Калиткин Н.Н., Панченко С.Л. Оптимальные схемы для жестких неавтономных систем. Математическое моделирование, 1999, т. 11, № 6, с. 52-81.

[10] Rosenbrock H. H. Some general implicit processes for the numerical solution of differential equations. Comput J., 1963, vol. 5, № 4, рр. 329-330.

[11] Федоренко Р.П. Введение в вычислительную физику. Москва, Изд-во МФТИ, 1994, 528 с.

[12] Галанин М.П., Милютин Д.С., Савенков Е.Б. Разработка, исследование и применение метода конечных суперэлементов для решения бигармони-ческого уравнения. Препринт ИПМ им. М. В. Келдыша РАН, 2005, № 59, 26 с.

[13] Галанин М.П., Савенков Е.Б. Метод конечных суперэлементов для задачи о скоростном скин-слое. Препринт ИПМ им. М. В. Келдыша РАН, 2004, № 3, 32 с.

[14] Галанин М.П., Савенков Е.Б. К обоснованию метода конечных суперэлементов. Журнал вычислительной математики и математической физики, 2003, т. 43, № 5, с. 713-729.

[15] Самарский А.А., Гулин А.В. Численные методы. Москва, Наука, 1989, 432 с.

[16] Hairer E., Norsett S., Wanner J. Solving ordinary differential equations. Part 2. Stiff and Differential-Algebraic Problems. New York etc., Springer, 1991.

[17] Галанин М.П., Ходжаева С.Р. Методы решения жестких обыкновенных дифференциальных уравнений. Результаты тестовых расчетов.

Препринт ИПМ им. М.В. Келдыша, 2013, № 98, 29 с.

Статья поступила в редакцию 24.12.2014

Ссылку на эту статью просим оформлять следующим образом: Галанин М.П., Ходжаева С.Р. Разработка и тестирование методов решения жестких обыкновенных дифференциальных уравнений. Математическое моделирование и численные методы, 2014, № 4, с. 95-119.

Галанин Михаил Павлович родился в 1956 г., окончил МГУ им. М.В. Ломоносова в 1979 г. Д-р физ.-мат. наук, зав. отделом ИПМ им. М.В. Келдыша РАН, профессор кафедры «Прикладная математика» МГТУ им. Н.Э. Баумана. Автор двух монографий и более 200 статей в научных рецензируемых журналах. e-mail: galan@keldysh.ru

Ходжаева Софья Рубеновна родилась в 1993 г., студентка МГТУ им. Н.Э. Баумана. Автор научной публикации. e-mail: dosd@yandex.ru

Development and testing for methods of solving stiff ordinary differential equations

© M.P. Galanin1, S.R. Khodzhaeva2

1 Keldysh Institute of Applied Mathematics of the Russian Academy of Sciences,

Moscow, 125047, Russia 2 Bauman Moscow State Technical University, Moscow, 105005, Russia

The paper is aimed at research of the (m,k)-method, CROS, finite superelement method and 4-stage explicit Runge-Kutta method for solving stiff systems of ordinary differential equations. Analysis of tests results showed that the best choice for systems with high stiffness is CROS. The finite superelement method is the «precise» method for solving linear systems of ordinary differential equations, the best supporting method for its implementation is (4,2)-method. The variation of the finite superelement method has been built and tested for solving nonlinear problems, this method proved to be unsuitable for problems with high stiffness.

Keywords: ordinary differential equations, numerical solution. stiff systems, finite superelement method, (4,2)-method, CROS.

REFERENCES

[1] Galanin M.P., Savenkov E.B. Metody chislennogo analiza matematicheskikh modeley [Methods of numerical analysis of mathematical models]. Moscow, BMSTU Publ., 2010, 591 p.

[2] Kalitkin N.N. Chislennye metody [Numerical methods]. Moscow, Nauka Publ., 1978, 512 p.

[3] Arushanyan O.B., Zaletkin S.F., Kalitkin N.N. Vychislitel'nye Metody i Pro-grammirovanie — Numerical Methods and Programming, 2002, vol. 3, pp. 11-19.

[4] Novikov E.A. Vychislitel'nye tekhnologii — Computational Technologies, 2007, vol. 12, no. 5, pp. 103-115.

[5] Kalitkin N.N. Poluyavnye skhemy dlya zadach bolshoy zhestkosti [Semi-explicit schemes for the problems of high rigidity]. Entsiklopediya nizkotemperaturnoy plazmy, seriya B [Encyclopedia of low-temperature plasma, series B], vol. VII-1, part. 1. Moscow, Yanus-К Publ., 2008, pp. 153-171.

[6] Al'shin A.B., Al'shina E.A., Kalitkin N.N., Koryagina A.B. Zhurnal Vychislit-el'noi Matematiki i Matematicheskoi Fiziki — Computational Mathematics and Mathematical Physics, 2006, vol. 46, no. 8, pp. 1392-1414.

[7] Al'shin A.B., Al'shina E.A., Limonov A.G. Zhurnal Vychislitel'noi Matematiki i Matematicheskoi Fiziki — Computational Mathematics and Mathematical Physics, 2009, vol. 49, no. 2, pp. 270-287.

[8] Shirkov P.D. Matematicheskoe Modelirovanie — Mathematical Models and Computer Simulations, 1992, vol. 4, no. 8, pp. 47-57.

[9] Kalitkin N.N., Panchenko S.L. Matematicheskoe Modelirovanie — Mathematical Models and Computer Simulations, 1999, vol. 11, no. 6, pp. 52-81.

[10] Rosenbrock H. H. Some general implicit processes for the numerical solution of differential equations. Comput J., 1963, vol. 5, no. 4, pp. 329-330.

[11] Fedorenko R.P. Vvedenie v vychislitelnuyu fiziku [An Introduction to Computational Physics]. Moscow, MIPT Publ., 1994, 528 p.

[12] Galanin M.P., Milyutin D.C., Savenkov E.B. Razrabotka, issledovanie I primenenie metoda konechnykh superelementov dlya resheniya bigarmonich-eskogo uravneniya [The development, research and application of the finite super element method for solving biharmonic equation]. Keldysh Institute of Applied Mathematics RAS, Preprint, 2005, no. 59, 26 p.

[13] Galanin M.P., Savenkov E.B. Metody konechnykh superelementov dlya zadachi o skorostnom skin-sloe [The method of finite superelement for the problem of high-speed skin layer]. Keldysh Institute of Applied Mathematics RAS, Preprint, 2004, no. 3, 32 p.

[14] Galanin M.P., Savenkov E.B. Zhurnal Vychislitel'noi Matematiki i Ma-tematicheskoi Fiziki — Computational Mathematics and Mathematical Physics, 2003, vol. 43, no. 5, pp. 713-729.

[15] Samarsky A.A., Gulin A.V. Chislennye metody [Numerical methods]. Moscow, Nauka Publ., 1989, 432 p.

[16] Hairer E., Norsett S., Wanner J. Solving ordinary differential equations. Part 2. Stiff and Differential-Algebraic Problems. New York, Springer, 1991.

[17] Galanin M.P., Khodzhaeva S.R. Metody resheniya zhestkikh obyknovennykh differentsialnykh uravneniy. Rezultaty testovykh raschetov [The methods of solving stiff ordinary differential equations. The results of test calculations]. Keldysh Institute of Applied Mathematics RAS, Preprint, 2013, no. 98, 29 p.

Galanin M.P. (b. 1956) graduated from Lomonosov Moscow State University in 1979. Dr. Sci. (Phys.&Math.), head of the Department at Keldysh Institute of Applied Mathematics of the Russian Academy of Sciences, professor of the Applied Mathematics Department at Bauman Moscow State Technical University. Author of two monographes and of more than 200 papers published in scientific peer-reviewed journals. e-mail: galan@keldysh.ru

Khodzhaeva S.R. (b. 1993), a student of the Bauman Moscow State Technical University. Author of one paper. e-mail: dosd@yandex.ru