CC BY
44
УДК 519.622
('А.Е. Новиков, В.В. Шайдуров
АЛГОРИТМ ИНТЕГРИРОВАНИЯ ПЕРЕМЕННОЙ КОНФИГУРАЦИИ НА ОСНОВЕ ЯВНО-НЕЯВНЫХ СХЕМ ЧЕТВЕРТОГО ПОРЯДКА1
Построены /.-устойчивый (4,2)-метод и явная пятистадийная схема типа Рунге-Кутты четвертого порядка точности. Создан алгоритм интегрирования переменного шага, в котором выбор эффективной численной схемы осуществляется на каждом шаге с применением неравенства для контроля устойчивости. Приведены результаты расчетов, подтверждающие эффективность построенного алгоритма.
Ключевые слова: жесткие задачи, явный и неявный методы, контроль точности и устойчивости.
А.Е. Novikov, V V. Shaidurov
AN INTEGRATION ALGORITHM OF VARIABLE CONFIGURATION BASED ON EXPLICIT-IMPLICIT SCHEMES OF 4™ ORDER OF ACCURACY
An /.-stable (4,2)-method of 4th order of accuracy and an explicit Runge-Kutta scheme of 4th order of accuracy are constructed. An integration algorithm of variable step size is formulated. The most effective numerical scheme is chosen for each step by means of stability control. The numerical results which confirm the effectiveness of the algorithm are given.
Key words: stiff problems, explicit and implicit methods, stability and accuracy control
Введение
Во многих важных приложениях возникает проблема численного решения задачи Коши для обыкновенных дифференциальных уравнений. В современных методах решения жестких задач при вычислении стадий применяется ¿^-разложение матрицы Якоби. В случае достаточно большой размерности исходной системы быстродействие алгоритма интегрирования фактически полностью определяется временем декомпозиции этой матрицы. Для повышения эффективности расчетов в ряде алгоритмов используется замораживание матрицы Якоби, то есть применение одной матрицы на нескольких шагах интегрирования [1]. Наиболее успешно этот подход применяется в многошаговых методах [2]. Не вызывает эта проблема особых трудностей и при построении алгоритмов интегрирования на основе других численных схем, если в них стадии вычисляются с участием матрицы Якоби в некотором итерационном процессе. В алгоритмах интегрирования на основе известных безытерационных численных схем, к которым относятся методы типа Розенброка [3] и их различные модификации [1], проблема замораживания более трудная. Следует отметить, что с точки зрения реализации безытерационные методы существенно проще алгоритмов на основе численных формул, в которых стадии вычисляются с применением итераций. Некоторым аналогом замораживания матрицы Якоби является применение в расчетах алгоритмов интегрирования на основе явных и /.-устойчивых методов с автоматическим выбором численной схемы. В этом случае эффективность алгоритма может быть повышена за счет расчета переходных участков явным методом [4]. В качестве критерия выбора эффективной численной формулы естественно применять неравенство для контроля устойчивости [5].
Здесь на основе /.-устойчивого (4,2)-метода и схемы Мерсона четвертого порядка точности построен алгоритм переменной структуры. Приведены результаты расчетов задачи проникновения помеченных радиоактивной меткой антител в пораженную опухолью ткань живого организма.
1. Исследование (4,2)-метода
Рассмотрим задачу Коши вида
у' = /(у), y(fo) = yo, *Q<t<tk, (!)
где у и/ — вещественные TV-мерные вектор-функции, t — независимая переменная. Для решения задачи (1) рассмотрим численную формулу [6]
Уп+1=Уп+Р\к1+'" + Р\кА, Dn=E~ahfn,
Dnk\ = hf(yn ),D„k2=ki, (2)
'Работа выполнена при финансовой поддержке РФФИ (проекты 11-01-00106 и 11-01-00224)
Опк3 = к/ (уп + [Ь\к\ +/¿¡2^2) + аЪ2к2 > 1)пк4 = к3 + Щ2к2 > где к - шаг интегрирования, а. р,. Д; и «(; — числовые коэффициенты, к,. 1</<4, — стадии метода, Е — единичная матрица, /п=д/!ду - матрица Якоби задачи (1). Подставим разложения к/ в виде рядов Тейлора в первую формулу (2). Полагая уп= у(!п) и сравнивая полученное представление с рядом Тейлора для точного решения, запишем условия четвертого порядка точности
1) Р1 + Р2 + 0 + аЪ2)РЪ + 0 + а32 + «42)^4 =1,
2) а^+2а^2+(а + /%1+/%2+3а«32)Л +
+(2а + у6?5^ + /?32 +4ассз2 +За«42)/>4 =0.5 ,
2 2 2 2 3) а /?1+За Р2+(а + 2а/% ^ + За/%2 + 6а 032)^3 +
+(3а +3а/?^ +4а/?з2 +Юа «32+ 6а а42)РА =1/6 ,
4) с?р\ + 4а3/?2 + (а3 + За2/% ^ + 6а2рт)2 + (3)
+10а3«з2)/?з + (4а3 + 6а2/%1 +10а2/%2 +
20а3о:з2 + 10а3о:42) = 1 / 24,
5) (/%1 + Рзг)2(Рз + ра) = 1/3,
6) а(/%1 +А2ХА1 +2Рз2)(рз +ра)=\1%,
7) а(^1+/?з2)2(0.5Л+^4) = 1/24,
8) (/%1+А2)3(/?3 + ) = 1 / 4 .
Применяя метод (2) для решения тестовой задачи
У = Лу, у(0) = У0^>0,
имеем условие ¿-устойчивости численной формулы (2) вида
а(а- Р]) + {Р31 ~а)ръ = 0.
Исследуя совместность этого соотношения и (3), запишем
76а2 - 29а + 3 146а2 + 89а -12 32а-4
Р\ о ’ ^*2 2 ’ РЗ ^-7 ’
27а2 27а2 27а
4-16а Л 48а-9 Л 9-24а
/>4= —, Рз\= ——, РЪ2=——> (4)
27а 32 а 32а
-54а2 + 57а -12 -864а3 + 828а2 - 288а + 36
«32 =-----------ô---> «42 =-------------------л-------,
8а-32а а(4-16а)
где а определяется из необходимого условия /.-устойчивости
24а4 - 96а3 + 72а2 - 16а +1 = 0 .
Данное уравнение имеет четыре вещественных корня
а\ = 0.10643879214266, а2 =0.22042841025921, а3 =0.57281606248213, а4 =3.10031673511599.
Для расчетов рекомендуется а=0.57281606248213, потому что в этом случае схема (2) дополнительно является ^-устойчивой.
Для контроля точности вычислений метода (2) четвертого порядка будем применять метод третьего порядка вида
Уп+1Д =Уп +hh +b2k2 +ЬЗкЗ + Ь4к5 , где DJc5=k4, а к,. 1</<3, определены в (2). Нетрудно видеть, что требования третьего порядка имеют вид
h +1>2 + (1 + «32)^3 + (1 + а32 + «42)^4 = 1 ; аЪу + 2аЬ2 + (а + +>632 +3a«32)¿3 +
(3а + ¡3^1 + /%2 +5ааз2 + Ааа42)Ъ4 =1/2 ,
111 1 a Ьу+За b2+(a + 2a(h,i+3a/332+6a «32)^3 +
+(6а2 + 4 a foi + 5а/332 + 15а^а32 + 10а^а42)Ь4 =1/6,
(/%1+/%2)2(*3 +Ь^) = ИЗ, где коэффициенты а, /?зь /А2. а32 и а42 заданы в (4). Данная система линейна относительно параметров Ьь 1</<4. При применении коэффициентов (4), имеем
by =1.203100567018353, Ъ2 =-6.552116304144386-10“1,
Ьъ = 7.115271884598151-10_1, Ь4 = -1.189345958672225-10“1.
Теперь оценку ошибки £„о можно вычислить по формуле г,„ „=|[у„ \- уп 1.11|. а при выборе шага проверять неравенство еп,о<£, где ||-|| - некоторая норма в IV. е - точность расчетов.
Оценку максимального собственного числа v„ 0=/2|A„ max| матрицы Якоби системы (1), необходимую для перехода на явную формулу, оценим через ее норму по формуле
vn,0 =11 3//ду ||= тах^хдд £^=1| dfi(уп)/dyj\.
2. Контроль точности и устойчивости метода Мерсона
Одним из самых эффективных явных методов типа Рунге-Кутта четвертого порядка точности является метод Мерсона [7]
27
Уп+1 - Уп +7^1 +Т^4 +7^5 >
6 3 6
h = hf(yn), k2 = hf(yn + \kx), k3 = hf(yn + + -U2) » (5)
3 6 6
13 13
k4 - ¥(Уп +~k\ +7^3 ). *5 = hf(yn +~kl - 7^3 + 2k4) ■
00 2 2
Пятое вычисление функции f не дает увеличение порядка до пятого, но позволяет расширить интервал УСТОЙЧИВОСТИ ДО 3.5 И оценить величину локальной ошибки ¿„4 с помощью kh то есть
5п^4 = (2k\ - 9к3 + 8Æ4 - к$ ) / 30 .
Для контроля точности будем применять неравенство ||^,4||<5е5/4, которое получено с учетом накопления глобальной ошибки из локальных ошибок [5]. Несмотря на то, что обоснование неравенства проведено на линейном уравнении, оно с достаточно высокой надежностью использовалось для решения нелинейных задач.
Теперь построим неравенство для контроля устойчивости. Применяя к разности k-—ki формулу Тейлора с остаточным членом в форме Лагранжа первого порядка, имеем
к3~к2 =h[df(jun)/dy](k2-kl)/6, где вектор ип вычислен в некоторой окрестности решения >’(/„). Учитывая, что
^-к^^/^/З + Оф3),
для контроля устойчивости (5) можно использовать неравенство
v„ 4 = 6 max | (к3 - к2\ / (к2 - k\)t |< 3.5 ,
\<i<N
где числу 3.5 равна длина интервала устойчивости. Отметим, что по мнимой оси область устойчивости также ограничена числом 3.5 (рис. 1). Введем обозначения r,„4 ¿„.VS. Тогда для контроля точности схемы (5) можно применять неравенство е„4<е5 4, а для контроля устойчивости следующее - v„4<3.5.
Рис. 1. Область устойчивости метода (5)
Так как оценка максимального собственного числа г^„=4=/г|Я„=тах| является грубой, то контроль устойчивости используется как ограничитель на размер шага интегрирования. В результате прогнозируемый шаг /?и+1 вычисляется следующим образом. Новый шаг И" по точности определим по формуле hac=q\hn, где кп есть последний успешный шаг, а q\, учитывая соотношение е„4=0{к5„), задается уравнением д5\епЛ<е. Шаг И1 по устойчивости зададим формулой hst=q2hn, где д2, учитывая соотношение Уп4=0(Ип), определяется из равенства <72^4=3.5. Тогда прогнозируемый шаг /2„+1 вычисляется по формуле
Ип+1 =тах[/2„,тт(/2а'с,/2^)].
Данная формула стабилизирует шаг на участке установления решения, где определяющую роль играет устойчивость.
3. Алгоритм интегрирования
На основе построенных методов легко сформулировать алгоритм переменного шага, в котором на каждом шаге выбирается наиболее эффективная численная схема. Расчеты всегда начинаются явным методом потому, что в нем не используется матрица Якоби. Нарушение неравенства у„4<3.5 вызывает переход на ¿-устойчивую схему. Передача управления явным методам происходит в случае выполнения неравенства у„0<3.5, где оценка V,,,, вычисляется через норму матрицы Якоби. Норма ||^|| в неравенствах для контроля точности вычисляется по формуле
|| £||= шах {]<§ \/(\угп |+г)}, (6)
1</<Аг
где 7 - номер компоненты, г - положительный параметр. Если по 7-й компоненте решения выполняется неравенство |[уп ||<г, то контролируется абсолютная ошибка ге, в противном случае - относительная ошибка £. Ниже построенный алгоритм переменного порядка и шага, а также с автоматическим выбором явной или ¿-устойчивой численной схемы будем называть МСМК4.
4. Результаты расчетов
Модель описывается системой двух дифференциальных уравнений в частных производных с начальными и граничными условиями. Это задача проникновения помеченных радиоактивной меткой антител в пораженную опухолью ткань живого организма [8]. Рассматривается система одномерных уравнений реакции-диффузии
ди д^и , ду
-кт— = -киу, (7)
а <Эх &
которые возникают из химической реакции А+В^С с константой скорости реакции к, где А - антитело с радиоактивной меткой, реагирующее с субстратом В - тканью, пораженной опухолью. Концентрации А и В обозначены через и и V соответственно. При выводе уравнений (7) предполагалось, что кинетика реакции описывается законом действующих масс, причем реагент А подвижен, тогда как реагент В неподвижен. Изучается полубесконечная пластина, внутри которой равномерно распределен субстрат В. Реагент А, попадая на поверхность пластины, начинает в нее проникать. Для моделирования проникновения уравнения (7) рассматриваются в полосе ЛУ {(х,/): 0<х<оо, 0 / 7} с начальными и(х, 0)=0, у(х. 0) у„. х>0 и граничными и(0, /) Ф(1). 0 / Т условиями, где V,, константа. Для численного решения переменная х преобразуется так, чтобы полубесконечная пластина преобразовалась в конечную. Такое преобразование обеспечивает специальное
семейство преобразований Мебиуса £=х!(х+с), с>0. Это замена переменных преобразует 8Т в прямоугольник {(С 1)'. 0<С< 1. О / 7 ). Через переменную С задача (7) переписывается в виде
ди (С-1)4 д2и 2(£-\)3 ди 8у ,
— = ,------+ ,---------киу,— = -кш (8)
Ы с д£2 с дt
с начальными и(£0)=0, у(С0) у„. С 0 и граничными и( 0, /) Ф(1). ди(\.1)/дС 0. 0 / /' условиями. Последнее граничное условие получено из соотношения би(/.1)1дх 0. Дискретизация производных по пространственным переменным с использованием метода прямых приводит к задаче Коши для системы обыкновенных дифференциальных уравнений. Для дискретизации применяется равномерная сетка {С]}, С=]- \С \С= 1 /Л'. 1</<Ж Через и у,- обозначены аппроксимации пЦ'.!) и уЦ'./). соответственно. Очевидно, что и} и у,- являются функциями от I. Дискретизации производных первого и второго порядков по пространственной переменной соответственно имеют вид
ди] и]+1~и]-1 д\ _и]_х-2и]+и]+х
“ 2лС • д(2 - (лС)2 ’ ■
Значения и0 и и у , получены из граничных условий, они имеют вид и.,, Ф(1). иы+1=иы. Полагая
у=(иь V], и2, у2, ..., м№ у^)т и Т=20, эта задача имеет вид
У=Ж,у), У(0) = 8, уеяш, 0<?<20, (9)
где N — параметр, а функция/определяется формулами
У2]+1~У2]-3 , 0 У2]-3~2У2]-1+У2]+1
/2 у-1 =а]----^---------+ Р]--------------о----------
2д^Г (д^)
-ку2¡-\У2у > /гу = -^2У-1,
где £=(0,Уо,0,Уо,...,0,Уо)1’, а/=2(/-АС-1)3с2, Д=(/-АС-1)4с2, Л<Г=Ш, 1</<ЛГ, У-\(0=^(0, Уш+\= Уш-\■ Функция Ф(0=2 при 0 /<5 и Ф(0=0 при 5 /<20. то есть Ф(() имеет разрыв первого рода в точке t=5. Согласно [8] подходящими значениями для параметров к, у0и с являются к 100, у0=1 и с=4. Расчеты проводились при N 400. то есть система (9) состоит из 800 уравнений. Задача о нахождении разрыва функции Ф(0 при / 5 возлагалась на алгоритм управления шагом. Решение задачи (9) приведено на рис. 2.
Рис. 2. Решение задачи (7)
Расчеты проводились с двойной точностью. Параметр г в формуле (6) выбирался таким обозом, чтобы фактическая точность вычисления решения совпадала с задаваемой точностью. Сравнение проводилось с приближенным решением, которое получено при задаваемой точности вычислений £=10 9 различными методами. Матрица Якоби вычислялась численно. В табл. 1 приведены результаты расчетов при различных значениях задаваемой точности вычислений е. Вычислительные затраты приведены в форме //(//). где через // и у обозначены соответственно число вычислений правой части и декомпозиций матрицы Якоби на интервале интегрирования. В табл. 1 используются обозначения: МС4 - метод Мерсона, МС4з1 - метод Мерсона с контролем устойчивости, МК4 - Ь-устойчивый метод, ЯКМК4 - алгоритм с автоматическим выбором численной схемы.
Таблица 1
Результаты расчетов
Метод / s RK4 RK4st MK4 RKMK4
10" 1 017 921 883 561 39 560(49) 35 286(34)
10~3 1 018 024 886 914 62 918(78) 44 923(46)
1£+ 1 018 177 889 604 76 717(95) 66 845(74)
10~5 1 018 701 882 836 92 165(114) 83 984(86)
10~6 1 019 268 916 280 103 687(128) 96 205(99)
Из сравнения результатов расчетов RK4 и RK4st следует, что на данной задаче контроль устойчивости приводит к сокращению вычислительных затрат более чем на 10%. Это является следствием устранений некоторых повторных вычислений решения, возникающих за счет неустойчивости численной схемы. На других задачах преимущество RK4st перед RK4 может достигать 50% [5]. Из сравнения результатов расчетов МК4 и RKMK4 следует сокращение числа декомпозиций матрицы Якоби, что является следствием расчета некоторых переходных участков по явной численной формуле.
Заключение
В RKMK4 с помощью признака можно задавать различные режимы расчета: 1) явным методом с контролем или без контроля устойчивости; 2) ¿-устойчивым методом с аналитической или численной матрицей Якоби; 3) с автоматическим выбором численной схемы. Все это позволяет применять данный алгоритм для решения как жестких, так и нежестких задач. При расчетах с автоматическим выбором численной схемы вопрос о том, является ли задача жесткой или нет, перекладывается на алгоритм интегрирования.
Литература
1. Хайрер Э., Ваннер Г. Решение обыкновенных дифференциальных уравнений. Жесткие и дифференциально-алгебраические задачи. - М.: Мир, 1999. - 685 с.
2. Byrne G.D., Hindmarsh А.С. ODE solvers: a review of current and coining attractions // J. of Comput. Physics. 1987. №70. P. 1-62.
3. Rosenbrock H.H. Some general implicit processes for the numerical solution of differential equations // Computer. 1963. №5. P. 329-330.
4. Новиков A.E., Новиков E.A. Численное решение жестких задач с небольшой точностью // Математическое моделирование. 2010. Т.22, №1. С. 46- 56.
5. Новиков Е.А. Явные методы для жестких систем. - Новосибирск: Наука, 1997. - 197 с.
6. Новиков Е.А., Шитов Ю.А., Шокин Ю.И. Одношаговые безытерационные методы решения жестких систем//ДАН СССР. 1988. Т. 301, №6. С. 1310-1314.
7. Merson R.H. An operational methods for integration processes // Proc. of Symp. on Data Processing. Weapons Research Establishment, Salisbury, Australia. 1957. P. 331.
8. Mazzia Iavemaro F. Test Set for Initial Value Problem Solvers // Department of Mathematics. University of Bari, August 2003. URL: http://www.dm.uniba.it/~testset.
Новиков Антон Евгеньевич, ассистент кафедры Сибирского федерального университета, тел. (391)2495979, e-mail: [email protected]
Шайдуров Владимир Викторович, член-корреспондент РАН, директор Института вычислительного моделирования СО РАН, тел. (391)290-74-76, e-mail: [email protected]
Novikov Anton Evgen ’evich, teaching assistant of the Siberian Federal University, tel. +7 (391)2495979, e-mail: aeno vikov@bk. ru
Shaidurov Vladimir Viktorovich, corresponding member of RAS, director of Institute of Calculation Mathematics of SB RAS, tel. (391)290-74-76^ e-mail: [email protected]
