Разработка (4,2)-метода третьего порядка точности для решения жестких задач

При моделировании динамических процессов в электрических сетях, радиоэлектронике, химической кинетике и других важных приложениях возникает необходимость численного решения задачи Коши для жестких систем обыкновенных дифференциальных уравнений [1; 2]. Стремление к все более точному описанию физических процессов приводит к росту размерности и жесткости соответствующей системы уравнений. В этом случае эффективность алгоритмов интегрирования существенно зависит от свойств устойчивости числен- ной формулы [3; 4]. Здесь получены коэффициенты L-устойчивого (4,2)-метода третьего порядка точности для численного решения жестких задач.

Класс (m,k)-методов. Рассмотрим задачу Коши для системы обыкновенных дифференциальных уравнений y′= f(y), y(t0)= y0, t0 ≤t≤tk, (1) где y и f - вещественные N-мерные вектор-функции; t - независимая переменная. Ниже будем предполагать, что задача (1) жесткая.

Зададимся значениями целых чисел m и k , k <  m и введем в рассмотрение следующие множества:

M m = {1 , 2 ,..., m },

M k = { m_t e M m | 1 = m 1 <  m ₂ <  • • <  m_k <  m }, M m — k = M m \ M k ,

Ji = { m - i ^{€ M}m I j >  1 , m j e M k , m <  i }, 2 <  i <  m .

Для решения задачи (1) будем применять ( m , k )-схемы вида [5]

m

• У п + 1 = У п ⁺ Е Pk ’ D n = E ^- ahf n ’ i = 1 i - 1

d a = hf ⁽ У п ⁺ Е р j k j ) ⁺ Е ^а y^kj , ^{i e} M k ,   ⁽²⁾

, ■= i                j ^e Ji

Dk, = k i + У а,, i e M„ t. n i       i — 1              ij j "            m—k j e Ji

Здесь E - единичная матрица; f_n' = dfy_n )/d y - матрица Якоби системы (1); h - шаг интегрирования; a , p_i , a j и ву - вещественные константы, определяющие свойства точности и устойчивости методов (2). В отличие о других известных одношаговых численных схем, в ( m , k )-методах правая часть задачи (1) вычисляется не для всех стадий. За счет этого значительно упрощается проблема замораживания матрицы Якоби, т. е. проблема использования одной матрицы D_n на нескольких шагах интегрирования. Для описания традиционных одношаговых методов достаточно одной константы m (число стадий). В данных схемах для описания затрат на шаг необходимо введение двух постоянных m и k . Вычислительные затраты на шаг интегрирования в методах (2) следующие: один раз вычисляется матрица Якоби и осуществляется декомпозиция матрицы D n , k раз вычисляется функция f и m раз осуществляется обратный ход метода Гаусса.

Устойчивость численных методов обычно исследуют на скалярной тестовой задаче

• У ' = ^У У , У ⁽ t о ) = У 0 , t ^ 0,

где X - произвольное комплексное число, Re(X) < 0. Применяя численную схему (2) для решения данной задачи, имеем у_{п +} 1 = Q ( x ) у_п , x = h X. Схема называется устойчивой при x = h X, если | Q ( x )| < 1. Область R комплексной плоскости называется областью устойчивости схемы, если она устойчива в каждой точке области. Схема называется A -устойчивой, если ее область устойчивости включает всю полуплоскость Re( h X) < 0. Схема называется L -устойчивой, если она A -устойчива и | Q ( h X)|^0 при Re( h X) ^ -да.

Максимальный порядок ( m ,2)-методов. Рассмотрим ( m ,2)-схемы следующего вида:

m

• У п + 1 = У п ⁺ Е P i k i , D n = E ^— a ^hf

i = 1

D n k 1 = hf ( У п ) , D n k l = k — 1 , 2 <  i <  5 1 — 1 ,       (3)

$ 1 ^{— 1}

D ^k = ^{hf (} У п ⁺ Е P $ 1 , j k j ) ^{+ a} $ 1 , $ 1 — 1 4 — 1 , j = 1

D„k, = kj i + a, „ ik, i, $i +1 < i < m, п i         i—1          i, $1 —1 $1 — 1’ 1

где m и $ i , $ i <  m - произвольные целые постоянные. Нетрудно видеть, что схемы (3) описывают всевозможные варианты ( m ,2)-методов.

Теорема. При любом выборе множеств M_k и J_i и при любом m нельзя построить ( m ,2)-метод выше четвертого порядка.

Без потери общности для простоты доказательство проведем для скалярной задачи (1), точное решение у ( t_{n +1} ) которой можно записать в виде

У ( t n + 1 ) = У ( t n ) + hf + ¹ h ² ff + ¹ h ³[ f f + ff ²] + 2       6

+ 24 h ⁴[ f ³ f + 4 ff'f ² + f ' f ³ ] + 1^0 h ⁵[ f 4 f + 4 ff f ³ +

+5 f,2f‘f2 + f "2f3 + fIVf4] + O (h6), где элементарные дифференциалы вычислены на точном решении y(tn). Учитывая, что ряд Тейлора для Dn-1 имеет вид

D;1 = E + ahf' + a2 h2 f'2 + a3 h3 f‘ +..., получим, что второе вычисление функции f(yп,c) будет осуществляться в точке

$ ₁ — 1

Уп,с = Уп +ЕР $1, jkj = j=1

= У п + Е C i h'f : ( ^{i — 1} ) f„ + O ( h ⁵) , i = 1

где числа су, 1 < i < 4 определяются через коэффициенты схемы (3), а элементарные дифференциалы вычислены на приближенном решении Уп. Учитывая ряд для точного решения у(tn + 1), для доказательства теоремы достаточно показать, что в представлении функции f(yп, c) в виде ряда Тейлора по степеням h не содержится слагаемое h5f„"2{п3. Разлагая f(yп,c) в ряд Тейлора в окрестности точки Уп до членов с h5 включительно, имеем f (Уп,с ) = hfn + С1 h 2 f; fn + h 3[ c 2 f‘ fn +

+ 2 Л"Л ²] + h ⁴[ c 3 f ;³ f + С 1 c 2 f ‘ ff ² +

+ ¹ c 3 ff ³] + h ⁵[ c 4 Г f + С 1 С 3 f² ff ² +

+2 c2 c 2 fXf3 + 24 c fIVf4]+O (h6), что завершает доказательство.

Исследование (4,2)-метода. Для решения задачи (1) рассмотрим численную формулу

Уп+1 = У п + Е Piki, D = E — ahf, i =1

D n k 1 = hf ( У п ), D n k 2 = k 1 ,                  (4)

Dk3 = hf (У п+1) + a32 k 2 , Dk4 = k3 + a42 k 2 , где схему

У п + 1 = У п +Р 31 k 1 +Р 32 k 2              (5)

называют внутренней схемой метода (4).

Для построения численной формулы третьего порядка точности разложим стадии k i , 1 <  i < 4 в ряды Тейлора в окрестности точки y n , т. е.

k = hf„ + ah2 f‘fn + a2 h3 f,'2 fn + a3 h4 f‘ f„ + O (h5), k2 = hfn + 2 ah2 f‘ fn + 3 a2 h3 fn2 fn + 4 a3 h4 f3 fn + O (h5), k3 = (1 + a32) hfn + (a + 3 a a32 +Рз1 +Рз2)h 2 ffn +

+ ( a ² + 6 a 2 a 32 + 2 a P 31 + 3 a P 32 ) h ³ f n ² f n +

+ 2( р 31 +P 32 )² h³ ff + ( a ³ + 10 a 0 + 3 a ^ +

+ 6 a ²P 32 ) h⁴ f n'³ f n + 2 a ( P 31 + P 32 )² h⁴ f n' f n f n ² +

+ a ( P 31 +P 32 )( P 31 + 2 P 32 ) h ⁴ f n ff n ² +

+ 2( р 31 +P 32 )³ h⁴ ff ³ + O ( h ⁵),          (6)

k 4 = (1 + a 32 + a 42 ) hf n + (2 a + 4 a 0 32 + 3 a 0 42 +

+ В,, + В, _э) h ² f’ f + (3 a ² + 10 a ²a_3? + 6 a ²a_4? +

¹ 3 1     ¹ 3 2 ^z n fl n fi                                32                42

+ 3 a P 31 + 4 a P 32 ) h³ f ‘ ² f n + 2( в 31 + P 32 )² h³ f n" f n ² +

+ (4 a + 20 a 0 32 + 10 a 0 42 + 6 a P 31 +

+ 10 a X) h ⁴ f n ³ f n + a ( P 31 +P 32 )² h ⁴ Л'ЛУп² +

+ a ( P 31 +P 32 )( P 31 + 2 P 32 ) h ⁴ Л"/,'/ ,² +

+ 2( p 31 +P 32 )³ h ⁴ f n" f n ³ + O ( h ⁵).

Подставим разложения стадий (6) в первую формулу (4). В результате запишем

Уп + 1 = Уп ⁺ [ Р 1 ⁺ Р 2 ^{+ (1 + a} 32 ) Р 3 ^{+ (1 + a} 32 ⁺

+ a ₄₂) p ₄] hf n + [ ap 1 + 2 ap ₂ + ( a + 3 a a ₃₂ +P ₃₁ +

+ P ₃₂) P ₃ + (2 a + 4 a a ₃₂ + 3 a a ₄₂ + P ₃₁ +

+ P 32 ) Р 4 ] h ² f n f n + [ a ² Р 1 + 3 a ² p 2 + ( a ² + 6 a ^ +

+ 2 a P ₃₁ + 3 a P ₃₂) p ₃ + (3 a ² + 10 a ² a ₃₂ + 6 a ² a ₄₂ +

+ 3 a p 31 + 4 a P 32 ) p ₄] h³ fn² f + ^( p 3 + p 4 )( p 31 + (7)

+ P 32 )² h ³ f n f ² + [ a ³ p + 4 a ³ p 2 + ( a ³ + 10 a ^ +

+ 3 a P 31 + 6 a P 32 ) p 3 + (4 a + 20 a 0 32 +

+ 10 a ³a₄ 2 + 6 a К + 10 a Ю p 4 ] h ⁴ fn³ f +

+ 2 a ( P 31 +P 32 )²( p 3 + 2 p 4 ) h⁴ fff ² +

+ a ( p 3 + p 4 )( p 31 +P 32 )( P 31 + 2 P 32 ) h⁴ f n f'f n ² +

+ 2( P 31 +P 32 )³( p 3 + p 4 ) h⁴ff ³ + O ( h ⁵) . 6

Условия порядка. Сравним полученное соотношение (7) с разложением точного решения у(tn + 1) в ряд Тейлора в окрестности точки tn до членов с h3 включительно. В результате получим условия третьего порядка точности схемы (4), которые имеют вид p1 + p 2 + (1 + 032) p 3 + (1 + 032 +042) p 4 = 1, ap1 + 2 ap 2 + (a + P31 + P32 + 3 a a32) p 3 +

+ (2 a + P ₃₁ +P ₃₂ + 4 a a ₃₂ + 3 a a ₄₂) p ₄ = 2 ,      (8)

a ² p 1 + 3 a ² p ₂ + ( a ² + 2 a P ₃₁ + 3 a P ₃₂ + 6 a ² a ₃₂) p ₃ +

+ (3 a ² + 3 a P₃₁ + 4 a P₃₂ + 10 a ²a₃₂ + 6 a ²a₄₂) p ₄ =2,

31         32            32          42   4     6

^{( в} 31 ⁺ P 32 ) ^{( p} 4 ^{+ p} 3 ) = 3 .

Неравенство для контроля точности вычислений построим по аналогии [3]. Для этого потребуем, чтобы локальная ошибка 5 _n метода (4) представлялась в виде

§ n =^ h ⁴ f n³ f n + O ( h ⁵) ,               (9)

где ^ - некоторая постоянная. С применением (6) и (7) нетрудно видеть, что требование (9) приводит к следующим дополнительным соотношениям на коэффициенты метода:

Ф31 + Р32) (p4 + p3) = 4, a (Р31 +р32)(р31 + 2р32)( p4 + p 3) = 1,

a ^{( в} 31 ^+в 32 ) ^{(2 p} 4 ^{+ p} 3 ) = 12 , при этом постоянная ^ вычисляется по формуле

^ = ^^ - a³p 1 - 4 a ³ p ₂ - ( a ³ + 3 a ²P₃₁ +

+ 6 a 2 P 32 + 10 a ³a 32 ) p 3 — (4 a ³ + 6 a ² ^ 31 +

+ 10 a 2 P ₃₂ + 20 a ³a ₃₂ + 10 a ³ a ₄₂) p ₄ .

Устойчивость. Исследуем устойчивость схемы (4) на линейном скалярном уравнении y' = Z у . Применяя (4) для решения у' = Z у , получим y_{n +} 1 = Q ₄ ( x ) y_n , x = h Z, где Q ₄ ( x ) есть дробно-рациональная функция. Мы не приводим Q 4 ( x ) в силу ее громоздкости. Из вида Q 4 ( x ) следует, что для L -устойчивости (4) необходимо выполнение соотношения

a ( a - p 1 ) + ( P 31 - a ) p 3 = 0 .            (12)

Теперь исследуем устойчивость промежуточной численной формулы (5). Применяя (5) для решения уравнения у' = Z у , получим у, + 1 = Q 2 ( x ) у_п , где Q 2 ( x ) имеет вид

Q 2 ( x ) =

1 + 0,25(3 - 8 a ) x + a ( a - P ₃₁) xx (1 - ax )²

Отсюда следует, что промежуточная схема (5) будет L-устойчивая, если в31 = a. Учитывая (12), имеем, что основная схема (4) и промежуточная (5) будут L-устойчивые, если p1 =Р31 = a.                     (13)

Коэффициенты метода. Исследуем совместность нелинейной системы алгебраических уравнений (8)

и (10) при условии (13). Из четвертого соотношения (8) и первого уравнения (10) следует, что β 31 + β 32 = 3/4. Учитывая (13), имеем β 32 = 3/4 -a . Из четвертого равенства (8) следует p 3 + p 4 = 16/27. Тогда из третьего уравнения (10) имеем p 4 = a^{- 1} (4 - 16 a )/27 и p 3 = a^{- 1} (32 a- 4)/27. Подставляя полученные значения во второе равенство (10), получим квадратное уравнение относительно коэффициента a , т. е.

^a 42

32 a — 48 a + 9 — 0 .

Данное уравнение имеет два вещественных корня: a ₁ = ¾ + 3√2/8 и a ₂ = 3/4 - 3√2/8. Отметим, что для известных четырехстадийных численных формул уравнение относительно a имеет четвертую степень. Требование L -устойчивости промежуточной схемы (5) позволило понизить порядок до второго.

Оставшиеся коэффициенты p 2 , α 32 и α 42 определим из трех первых уравнений (8). В результате получим следующие коэффициенты метода (4):

P1 — P31 — a, P32 — 4 — a,

3 — 57 a + 312 a ² — 528 a³ 54 a ²(4 a — 1)²

a (4 a — 1)²( — 162 a ² + 196 a — 11) 54 a ²(4 a — 1) 2

32 a — 4      4 —16 a p 3 —       , p 4 — ~

27 a         27 a

^a 32

432 a ³ — 1040 a ² + 392 a — 24

64 a (4 a — 1)

a 42

— 216 a ⁴ + 284 a ³ — 154 a ² + 44 a — 3

8 a (4 a — 1)²

Подставляя сюда корни уравнения (14), получим два набора коэффициентов вида

33 a — —+ —

65 281

p₂ —----

2    324 648

64 16

81 81

16 16        33

p4 —----V2. P31 ——i— V2,

4    81 81     31 48

85 395

—

32       8

—

49 —---1--

32     16 128

2121 3095

a 42 — ""64"—

2461 145 2

--1--,

3072   256

P 32   p

^a 32

85 395 ---n 2.

16 128

2121 3095

--1--

64    128

72, 2

2461 145 2

3 072    256 ’

при которых схема (4) имеет третий порядок точности и является L -устойчивой вместе с промежуточной численной формулой (5).

Контроль точности вычислений. Согласно [6], для контроля точности вычислений метода (4) с локальной ошибкой (9) можно использовать ошибку ε n вида

£ n - h ³ f;² f n + о ( h ⁴) ,

где постоянная ξ определена формулой (11). Величину ε _n оценим следующим способом. Определим коэффициенты b_i , 1 ≤ i ≤ 4, из условия выполнения соотношения

S bk — h ³ f' f n + O ( h ⁴). i — 1

Используя разложения (6) стадий ki, 1 ≤ i ≤ 4, в ряды Тейлора в окрестности точки yn до членов с h3 включительно и приравнивая соответствующие соотношения на коэффициенты нулю, получим относительно bi, 1 ≤ i ≤ 4, линейную систему алгебраических уравнений b1 + b 2 + (1 + a32) b3 + (1 + a32 + a 42) b4 — 0, ab1 + 2 ab 2 + (a + P31 +P32 + 3 a a32) b3 +

+ (2 a + P ₃₁ +P ₃₂ + 4 a a ₃₂ + 3 a a ₄₂) b ₄ — 0 ,      (17)

a b + 3 a b 2 + ( a + 2 a P 31 + 3 a P 32 + 6 a a 32 ) b 3 +

+ (3 a ² + 3 a P ₃₁ + 4 a P ₃₂ + 10 a ² a ₃₂ + 6 a ² a ₄₂) b ₄ — 1 ,

( P 31 +P 32 ) 2 ( b 4 + b 3 ) — 0 .

Исследуя совместность (17), получим

b 4 — —2-------2-------

8 a a ₃₂ + 4 a a ₄₂ + 3 a

b3 — — b4,

b ₂ — — (1 + a ₃₂ + 2 a ₄₂) b ₄ , b 1 — (1 + a ₃₂ + a ₄₂) b ₄.

Подставляя сюда коэффициенты (15) и (16), получим два набора значений b i , 1 ≤ i ≤ 4 вида

3 968 960 — 2 745 760^

53 105 312V2 — 76 687 712

и следующего

33 a —---V2

65 281

64 16

,   698 36872 — 217 088

b —--------------------,

3       2987901

324 648

16 16

81 81

3 81 81

—---V2

31 48

,   217 088 — 698 36872

b —-------------------,

4       2987901

и следующего

,   3 968 960 + 2 745 76072

b —--------------------------,

1          331989

,    76 687 712 + 53105 31272

b =--,

2          2987901

217 088 + 69836872

2 987 901

217 088 + 698 36872

2987901     .

Теперь, согласно [6], в неравенстве для контроля точности вычислений можно применять оценку ошибки вида в n =7Е bk. (18) i=1

Отметим особенность построенной оценки. Для метода (4) в силу L-устойчивости выполняется соотношение Q4(x) → 0 при x → –∞. Так как для точного решения y(tn + 1) = exp(x)y(tn) задачи y′ = λy, y(t0) = y0 выполняется аналогичное свойство, то естественным будет требование стремления к нулю оценки (18) при x → –∞. В силу того что оценивается главный член ошибки, т. е. первый член при разложении ошибки в ряд Тейлора, для (18) это требование не выполняется. Поэтому вместо (18) будем контролировать оценку вида в n (j„)=d в n, 1 < к < 2, (19)

при этом имеет место ε n (2) → 0 при x → –∞. В результате для контроля точности вычислений и при выборе шага можно применять неравенство

II в n ( j n ) ||<в, 1 <  j n <  2 , (20)

где ε – требуемая точность расчетов. Нетрудно видеть, что в смысле главного члена оценки ε n и ε n ( j n ) совпадают при любом значении j n . Отметим, что применение (19) вместо (18) не приводит к значительному увеличению вычислительных затрат. При x → 0 оценка ε n (1) = ε n правильно отражает поведение ошибки и нет смысла проверять (20) при других значениях j n . При резком увеличении шага поведение ошибки ε _n может оказаться неудовлетворительным, что проявляется в повторных вычислениях решения. Поэтому если требуемая точность не выполняется, имеет смысл проверить (20) при j n = 2.

Заметим также, что применение (20) для контроля точности вычислений не приводит к дополнительным вычислениям правой части и матрицы Якоби, потому что в данном неравенстве используются ранее вычисленные стадии k i , 1 ≤ i ≤ 4.

Результаты расчетов. Численный эксперимент проводился для десяти тестовых примеров [7]. В качестве критерия эффективности алгоритмов выбрано if – число вычислений правой части и ij – число обращений матрицы D_n на интервале интегрирования. Время вычисления не приводится, потому что для некоторых оно настолько мало, что не может служить объективной характеристикой эффективности метода.

Целью данного численного эксперимента является сравнение эффективности методов вида (4) с допол- нительным требованием L-устойчивости промежуточной численной формулы (5) и без данного требования [8]. В конкретных расчетах левая часть ||εn(jn)|| неравенства (20) вычислялась по формуле

|| в n (jn )ii= max rj, 1< i < N | yn i + r где N - размерность задачи (1); r - некоторая положительная постоянная. Если по i-й компоненте решения выполняется неравенство |yni| < r, то контролируется абсолютная ошибка εr, в противном случае – относительная ошибка ε.

Из анализа результатов расчетов следует, что алгоритм интегрирования с коэффициентами (15) эффективнее алгоритма с коэффициентами (16). Поэтому ниже будем рассматривать только метод (4), (15), который назовем ODE_A. Алгоритм интегрирования на основе схемы (4) без требования L -устойчивости (5) назовем ODE_B.

Для решения десяти примеров [7] алгоритму ODE_A потребовалось 748 вычислений правой части и 345 декомпозиций матрицы Якоби при точности ε = 10 ^{- 2} и if = 4 007 и ij = 2 231 при точности ε = 10 ^{- 4} . Для алгоритма ODE_B вычислительные затраты следующие: if = 916, ij = 471 при ε = 10 ^{- 2} и if = 4743, ij = 2534 при ε = 10 ^{- 4} . В конце интервала интегрирования фактическая точность вычислений не хуже задаваемой точности. По суммарным затратам эффективность обоих алгоритмов отличается не сильно. Однако на некоторых задачах, на которые приходится небольшой процент вычислительных затрат и жесткость которых достаточно велика, алгоритм ODE_A эффективнее ODE_B примерно в 1,3 раза.