Numerical methods for super stiff ODE

УДК 519.622.2

ЧИСЛЕННЫЕ МЕТОДЫ ИНТЕГРИРОВАНИЯ СВЕРХЖЁСТКИХ ОДУ

Власов Сергей Александрович

Магистр техники и технологии по направлению «Системный анализ и управление», фрилансер;

E-mail: vlasovs2@yandex.ru

Работа посвящена исследованию методов интегрирования задачи Коши для сверхжёстких систем

обыкновенных дифференциальных уравнений. Исследовались следующие методы: полностью неявные и

диагонально неявные методы Рунге-Кутты, методы Розенброка, коллокационные методы. Поиск решения

производится на основе вычислений повышенной точности.

Было выполнено экспериментальное сравнение десятка различных методов, и полученные результаты

представлены в виде графиков зависимостей "настраиваемая точность — точность полученного решения" и

"точность полученного решения — трудоёмкость".

В работу входит описание методов поиска коэффициентов непрерывных методов Розенброка до

точности 30 и более значащих цифр.

Ключевые слова: сверхжёсткие обыкновенные дифференциальные уравнения, непрерывные методы

Розенброка, вычисления повышенной точности.

NUMERICAL METHODS FOR SUPER STIFF ODE

Vlasov Sergey Aleksandrovich

Master of technology and technologies, freelancer;

E-mail: vlasovs2@yandex.ru

This work is devoted to the study of integration methods for the Cauchy problem for super stiff systems of

ordinary differential equations. The following methods were investigated: FIRK and DIRK methods, Rosenbrock

methods, and collocation methods. The search for a solution is based on calculations of high precision

computations.

An experimental comparison of a dozen different methods was carried out, and the results obtained are

presented in the form of graphs of dependencies "adjustable accuracy — accuracy of the obtained solution" and

"accuracy of the obtained solution — computations time".

The work includes a description of methods for finding the coefficients of multirate Rosenbrock methods up to

an accuracy of 30 or more significant digits.

Keywords: super stiff ordinary differential equations, multirate Rosenbrock methods, high precision

computations.

Введение

Численное интегрирование систем обыкновенных дифференциальных уравнений (ОДУ) — прикладная

задача, относящаяся к традиционному математическому моделированию.

Жёсткой системой называется такая система ОДУ, численное решение которой явными методами

(например, методами Рунге – Кутты или Адамса) является неудовлетворительным из-за резкого увеличения

числа вычислений. Такие системы часто решают неявными методами или явно-неявными схемами типа

Розенброка, которые дают обычно несравненно лучший результат, чем явные методы. Сверхжёсткими

называют системы ОДУ со степенью жёсткости более 10

[1].

Стандартная точность чисел с плавающей точкой современных компьютеров: от 15 до 20 цифр в

мантиссе. Этого бывает недостаточно, особенно в плохо обусловленных задачах и в задачах с внутренними

пограничными слоями (контрастными структурами) [2]. В них ошибки начальных данных, а также ошибки

округления, неизбежно возникающие в ходе расчёта, приводят к огромной потере точности. Для них нередко

приходится использовать расширенную точность (30 и более цифр), чтобы единичные округления были как

можно меньше.

При использовании методов Розенброка можно столкнуться с проблемой, когда коэффициенты

приведены с малым числом значащих цифр, рассчитанных на вычисления со стандартной точностью. При

этом приходится обращаться к условиям порядка и устойчивости, чтобы уточнить коэффициенты методов.

Иногда требуется, чтобы решение было непрерывным. В отличие от разностного решения непрерывное

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

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

непрерывный метод можно воспользоваться интерполяцией узлов, или построить решение на основе стадий

численного метода, выбранного для решения системы.

В данной работе описан метод уточнения коэффициентов некоторых методов Розенброка, а так же даны

сами уточнённые коэффициенты, в том числе и для непрерывных методов на основе стадий.

Так же в работе проводится экспериментальное тестирование и сравнение десятка различных методов, и

полученные результаты представлены в виде графиков зависимостей "настраиваемая точность — точность

полученного решения" и "точность полученного решения — трудоёмкость". Тестирование было проведено на

классических сверхжёстких задачах с использованием расширенной точности с получением непрерывного

решения.

1. Современные вычисления

Многие вычисления удобно делать в какой-нибудь системе компьютерной математики, такой как Maple.

Система Maple предназначена для символьных вычислений, обладает развитыми графическими средствами и

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

произвольной точности с различным количеством значащих цифр.

Ещё активно использовался язык программирования C/C++ и С++ Builder 6 как программная среда

разработки.

1.1. Длинная арифметика

Основная масса библиотек длинной арифметики написана под Unix. К ним относится и библиотеки gmp

(GNU Multi-Precision Library) и mpfr (multiple-precision floating-point computations with correct rounding),

которые позволяют делать вычисления над числами произвольной точности. Библиотека gmp считается самой

быстрой библиотекой чисел произвольной точности.

Под Windows не существует библиотек подобных gmp и mpfr. Но можно найти описание методов,

позволяющих считать числа точности double double [3](как минимум 31 значащая цифра) и quad double [4]

(как минимум 62 значащие цифры). Числа double double работают быстрее, чем числа mpfr той же точности. В

основе обычно используется 64-х битный тип double, но его можно заменить 80-ю битным типом long double

(extended) в 32-х битном приложении, что даст больше значащих цифр. При этом нужно изменить константу

разбиения 134217729 на 8589934593, а так же заменить некоторые формулы (для экспоненты, синуса и гамма-

функции) на более точные. Однако и приложение при этом будет 32-х битное, так как в 64-х битном

приложении нет 80-ти битного типа.

1.2. Получение условий порядка

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

подобные слагаемые и группировать слагаемые относительно матричных операторов. Риск ошибиться при

этом значительно уменьшается. Так получены условия порядка для всех методов Розенброка, которые

включены в данную работу.

Рисунок 1. Программа получения условий порядка.

1.3. Метод Ньютона для решения некорректных задач

Имеются в виду задачи нахождения решения системы нелинейных уравнений, в которых количество

переменных превышает количество уравнений.

Пусть дана система уравнений:

xxxF















0),...,,(

...

0),...,,(

212

211

(1)

Эта система дополняется недостающими уравнениями, например, вида 0=0. Далее в методе Ньютона

выбирается начальное приближение и вычисляется матрица Якоби. Матрица Якоби вырождена и имеет

определитель равный нулю, поэтому нельзя перейти к следующей точке по стандартной схеме:

iii

xxx

xFxx







1

),()(

(2)

где )(



— матрица Якоби.

Воспользуемся регуляризацией системы (1) по Тихонову.

Пусть матрица Якоби равна A, тогда приближённое решение (2) можно найти из следующей системы

уравнений:

*)( vfAvEAA





(3)

где



— параметр регуляризации,

— единичная матрица,

xv  , )(

xFf  , *v — желаемое

x .

Примем 0*



v . Параметр регуляризации



влияет на скорость сходимости получившегося метода

Ньютона, чем он больше, тем сходимость медленнее и может вообще остановиться. С другой стороны, если

взять очень малый параметр, то система (3) может дать очень крупные значения и метод Ньютона разойдётся.

Возможное значение параметра







1.4. Непрерывный метод Розенброка на основе стадий

Чтобы метод был непрерывным, добавляют дополнительные условия порядка относительно параметра θ,

зависящего от времени [5]. Решение ищется в виде:

 

 

 



jii

pky

1 1



(4)

где

k — стадии метода Розенброка,

— параметры задаваемые методом,









1..0



 nnn

tttt



n — количество стадий метода, m — максимальная степень θ.

Максимальная степень m задаётся вручную. Условия порядка записываются относительно параметра θ.

Например, должны выполняться четыре условия, что даст 3-й порядок для автономных систем. Тогда

необходимо выполнения m



4 условий. Если параметров в (4) больше, то ещё

условий должно

выполняться, для того, чтобы (4) совпала с основной формулой при 1





. Если параметров больше, чем

)4( nm





, то можно воспользоваться методом Ньютона с регуляризацией СЛАУ.

Для методов с комплексными коэффициентами решение ищется в виде:

 

 

 















jii

pky

1 1



(5)

где

k — комплексные стадии,

— комплексные параметры, 1..0





2. Примеры уточнённых методов Розенброка

Для взятых из литературы методов было проведено уточнение коэффициентов до точности 30 и более

значащих цифр. Так же приводятся коэффициенты непрерывных методов основанных на стадиях данных

методов.

2.1. (4,2)-метод

Данный метод взят из [6]. Это L-устойчивый метод четвертого порядка. Метод задается следующими

формулами:

iinn

fhaEDkpyy











)(

1 nn

yfhkD  ,

kkD

 (6)

2322321313

)( kkkyfhkD





24234

kkkD





где k

— стадии метода, E — единичная матрица, f(y) — правая часть системы, h — шаг интегрирования.

Уточнённые параметры даны в таблице 1. В ней так же приводится непрерывный метод, основанный на

формуле (4).

Таблица 1.

Параметры (4,2)-метода

a=0.5728160624821348554080013849767683409315

p1=1.278369390124472506000782151942886989045

p2=-1.007386809804384747838093468437643128000

p3=0.9265539109395042110093604870343782187148

p4=-0.3339613183469116184167678944417856261224

β31=1.009004690299215025588082799604327399592

β32=-0.2590046902992150255880827996043273995917

α32=-0.4955220641657818341715530456206093602956

α42=-1.287776482339217217685184223893417502206

Коэффициенты для непрерывного метода 3-го порядка

p11=2.807940918382931280268331126232341095630

p12=-2.283352377926245217075188710786266412373

p13=0.7537808496677864428076397364968123057875

p21=-3.136230519398356699219158918564434846756

p22=4.602228694452115124287372736689335424252

p23=-2.473384984858143172906307286562543705497

p31=1.031459744165087582629789670058674398639

p32=-1.800682298774149742460551733666582983217

p33=1.695776465548566370840122550642286803293

p41=-1.031459744165087582629789670058674398639

p42=1.800682298774149742460551733666582983217

p43=-1.103183872955973778247529958049694210700

2.2. CROS2_1

Данный метод взят из [7]. Это двухстадийный L1-устойчивый метод с комплексными коэффициентами

4-го порядка. Метод задается следующими формулами:





22111

Re kbkbyy









)()(

11 nny

yfkyfE 



(7)





))Re(())Re((

1212

kcyfkkayfE

nny





где k

— комплексные стадии метода, E — единичная матрица, f(y) — правая часть системы, τ — шаг

интегрирования.

Уточнённые параметры даны в таблице 2. В ней так же приводится непрерывный метод, основанный на

формуле (5).

Таблица 2.

Параметры CROS2_1

α1= 0.09705048233513176225872820658157096372762+0.1441824711215367963888908188125269341956i

α2= 0.1886638033791539520269860791327147505581+0.06177441689689081715947141717106675420030i

b1= 0.04833419895509438083025798424226646350104-0.3205959705202491275579925812180370610420i

b2= 0.9516658010449056191697420157577335364990-1.696774337833590152863879033397701636829i

c= 0.1730887968652108268236509028841256629989-0.1694095699539015207491928002966606686529i

a= 0.5359744564304914728722121662633122018534-0.9665922748484189311614005696233999104339i

Коэффициенты для непрерывного метода 3-го порядка

p11=3.619032671469712553690516593029702603963-0.2004225284051934625395768416824377885643i

p12=-8.414446230029770923182842133180291069678-0.6613406243058764343155995103891857483354i

p13=4.843747757515152750322583524392854929215+0.5411671821908207692971837708535864758580i

p21=-2.619032671469712553690516593029702603963-9.183665759393742328007731132919548994127i

p22=8.414446230029770923182842133180291069678+29.50534469034316892770669484128541508572i

p23=-4.843747757515152750322583524392854929215-22.01845326878301675256284274176356772842i

2.3. CROS2_Optimal

Данный метод взят из [8]. Это двухстадийный метод с комплексными коэффициентами 4-го порядка. Он

имеет оптимальное L – затухание как для функции устойчивости, так и для функции внутренней

устойчивости. Метод задается следующими формулами:





22111

Re kkyy











)()(

11 nny

yfkyfE 



(8)





12112112122

))Re(())Re(()( kkyfkyfkyfE

nynny





где k

— комплексные стадии метода, E — единичная матрица, f(y) — правая часть системы, τ — шаг

интегрирования.

Уточнённые параметры даны в таблице 3. В ней так же приводится непрерывный метод, основанный на

формуле (5).

Таблица 3.

Параметры CROS2_Optimal

α1=0.4573733434972975729886712378597810560027-0.2351004879985426732069962332538443381790i

α2=0.04262665650270242701132876214021894399728-0.3946329531721133798109289632446957436482i

γ21=0.9147466869945951459773424757195621120054+0.6546908396281089262421198339615985775839i

δ21=0.7451903446015528669012176569684448709433+3.042955553310975129001697189423158388626i

π21=-0.2317220370463119616442110173716549151421+0.07081125027258135212790263920701540770149i

ω1=0.2784065608064581775941260630978839609622-0.9196225438536243933605777235086716993250i

ω2=0.7215934391935418224058739369021160390378+0.2014773772275646399542954916451697899548i

Коэффициенты для непрерывного метода 3-го порядка

p11=1.662352335632077153904266284657282585433+3.058378731739278701194655437881970784617i

p12=-1.145501425793116439658894990826412127076-6.527071235918501156496073309360094210372i

p13=-0.2384443490325025366512452307329864973956+2.549069960325598061940840147969451726429i

p21=-0.6623523356320771539042662846572825854333-3.084120463545428936697275176270423611394i

p22=1.145501425793116439658894990826412127076+5.333814343596499665083881944741250473414i

p23=0.2384443490325025366512452307329864973956-2.048216502823506088432311276825657072065i

2.4. RODAS

Это метод Розенброка из [9]. Он 6-ти стадийный, вложенный, жёстко точный.



















jjnn

jij

jijni

kbyy

sikJhkyfhk

..1,)(



, (9)

где J — матрица Якоби, h — шаг интегрирования, k

— стадии метода.

Данный метод можно распространить и на неавтономные задачи с теми же коэффициентами

(Подробности см. в [9]).

Уточнённые параметры даны в таблице 4. В ней так же приводится непрерывный метод, основанный на

формуле (4). При этом:





ijijij



 ,





— вложенный метод 3-го порядка,





,bb

— основной метод 4-го порядка.

Построение непрерывного метода рассматривается в [5].

Таблица 4.

Параметры метода RODAS

γ=0.25

α21=0.75

α31=0.08612040081415230810139697878760241798691

α32=0.1238795991858476818986030212123975820131

α41=0.7749345355073239835959759197344993729071

α42=0.1492651549508681389613302285878078259563

α43=-0.2941996904581921225573061483223071988633

α51=5.308746682646141097685244508085578815512

α52=1.330892140037268841059079254823776347387

α53=-5.374137811655560533686603914331557047750

α54=-0.2655010110278494050577198485777981151496

α61=-1.764437648774481929301295441577120895355

α62=-0.4747565572063022542415003747471682133343

α63=2.369691846915801042868069955427380559080

α64=0.6195023590649831287129733353848142171276

α65=0.2500000000000000119617525255120943324826

β21=0.0

β31=-0.04939199999999976379720604242479557432309

β32=-0.0141120000000002362027939575752044256769

β41=-0.4820494693877553295338356343641135140289

β42=-0.1008795555555553188031270166574940582696

β43=0.9267290249433106483369626510216075722987

β51=-1.764437648774484950984492409473177979061

β52=-0.474756557206303164236822915996783772739

β53=2.369691846915804789575079569305906277052

β54=0.6195023590649833256462357561640554747468

β61=-0.080368370789111680412472246275789336523

β62=-0.0564906135924470382928821126158819718130

β63=0.488285630042796820218631292635299492216

β64=0.5057162114816190413438659233992289589771

β65=-0.1071428571428571428571428571428571428569

Коэффициенты для непрерывного метода 3-го порядка

p11=0.4195946799873075119058720335388621401299

p12=-0.7335246808281686702714286385831914932339

p13=0.2677412801216726048848534454489047125024

p14=-0.03417965006992312693176908668036469592144

p21=0.2547354250795645221587858725988736232911

p22=-0.5928167839085814444337143623649121710159

p23=0.2900798161885619002781701285533800486149

p24=-0.008489070951992016296123751403223472703150

p31=1.051349424207045854163844526228495743638

p32=-1.025245760555795036484891571562738176784

p33=0.4239971915887983504197795780790293924265

p34=0.03818477480274765211989875989051253293510

p41=-0.8887851816599769967853215338953422550436

p42=3.271289083573963740495968694592938222440

p43=-1.879097079219423171177154225473452234601

p44=0.002309388787055468810372988175085226182488

p51=-0.0007361278248275885787333063783016651992557

p52=-0.5586980766585664375101099181839054056919

p53=0.2721840821423380149400326795532119142342

p54=0.1801072651981988682916676878661380138000

p61=0.1638417802108866971355524079074124131835

p62=-0.3610037816228521517958242038981909757140

p63=0.6250947091780523006543183938389261668234

p64=-0.1779327077660868459940465978481476042930

2.5. RODAS5

Это метод Розенброка из [7], основанный на тех же формулах (9), что и предыдущий. Он 8-ми стадийный,

вложенный, жёстко точный. Построение данного метода можно найти в [10]. При этом:





ijijij



 ,





— вложенный метод 4-го порядка,





,bb

— основной метод 5-го порядка.

Уточнённые параметры даны в таблице 5. Рабочий непрерывный метод построить не получилось.

Таблица 5.

Параметры метода RODAS5

γ=0.19000000000000000000000000000000

α21=0.38000000000000003358592000427743

α31=0.18998096323048065165930209720647

α32=0.18621536159480624834069790279353

α41=0.083489136955726187158204690432470

α42=0.69567870317296510697170912192462

α43=-0.27150977576330429412991381235709

α51=0.019749954685320738490462476406881

α52=1.1527680537682619629021698267018

α53=-0.58172139560641273372994044079831

α54=-0.11569326366863669068176690949660

α61=0.83441493842466974700123546967678

α62=1.4341150127667629897852468982053

α63=-2.7463052995987880016901240308337

α64=4.6363567009884130000000000000000

α65=-3.1585813525810578672248231125632

α71=3.4243826682568668406794981276729

α72=-9.0682947038850302879842814456849

α73=3.6910912705994596427861075132406

α74=5.9760755264645538521206406306885

α75=-3.2132547614358500195546952967799

α76=0.18999999999999997195273047086282

α81=0.37542948780457440998356143171437

α82=-0.46222173916448809079146613083026E-16

α83=-0.19854963795383124505629117057112

α84=1.3845456784046262637970502755724

α85=-0.69859676359770365574433060778458

α86=-0.052828764657665387985155187551153

α87=0.18999999999999966122733917506885

β21=-0.00384213367066316653924755891158

β31=-0.06587508126936369036787904015459

β32=-0.06331442359835370063212095984541

β41=0.53487908306342626417680465966545

β42=-0.47157078286588050578677711751138

β43=0.44331311780172184160997245784594

β51=0.73026015628341180828986744008974

β52=-0.2532267377201464585976956105773

β53=0.37740086754994634597527595020071

β54=-0.098347581382224195667447779713122

β61=3.4243826682568670302497950472480

β62=-9.0682947038850299925361899111267

β63=3.6910912705994601310908139407904

β64=5.9760755264645538704774440136053

β65=-3.2132547614358500001592354994989

β71=0.3754294878045746740337485322662

β72=-0.781123376542690E-16

β73=-0.1985496379538313447723579795563

β74=1.3845456784046259182088969764613

β75=-0.6985967635977040646521761136264

β76=-0.05282876465766510470577376127586

β81=0.37542948780457467127309153966693

β82=-0.71947341385872637662088537334586E-16

β83=-0.19854963795383134239547860424177

β84=1.3845456784046259068336921183261

β85=-0.69859676359770406041453472873944

β86=-0.051105182319437447937208832186395

β87=-0.00172358233822765541222010695281

2.6. MyRosenbrock

Это авторский метода, основанный на тех же формулах (9), что и предыдущий.

Метод 4-х стадийный, имеет 4-й порядок и он жёстко точный. При этом:





44332211

Уточнённые параметры даны в таблице 6. В ней так же приводится непрерывный метод, основанный на

формуле (4).

Таблица 6.

Параметры метода MyRosenbrock

γ=0.2204284102592123180415434126292381629422

α21=-0.5205432483291658413100814077709160752376

α31=0.6075130779487283175846362337845170732959

α32=0.1009846139765253545887270473713338894896

α41=0.1641108001779867014074129836895915640195

α42=0.1293432375977119743143904143048741327283

α43=0.7065459622243013242781966020055343032522

γ21=1.041916673708566646332986903249560906166

γ31=-0.01906306867082466862784767220621635695740

γ32=-0.4883505673094258925645673717663805516965

γ41=0.3208618397035608796737063652926844741367

γ42=0.02205957636947700893965888505141667291981

γ43=-0.5633498263322502066549086629733393099986

b1=0.4849726398815475810811193489822760381562

b2=0.1514028139671889832540492993562908056481

b3=0.1431961358920511176232879390321949932536

b4=0.2204284102592123180415434126292381629422

Коэффициенты для непрерывного метода 3-го порядка

p11=1.614220168082677875270156338410046476717

p12=-1.773522416520713007296954629873264838965

p13=0.6442748883195827131079176404454944004042

p21=-0.7730850179266677645012163439058016182438

p22=2.000378477754902478764580585880475653431

p23=-1.075890645861045731009314942618383229539

p31=-0.1016288971935796221654533100265442181350

p32=0.6328462020633125972007704371496734160304

p33=-0.3880211689776818574120291880909342046418

p41=0.2604937470375695113965133155222993596616

p42=-0.8597022632975020686683963931568842304964

p43=0.8196369265191448753134264902638230337770

3. Тестирование

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

экстраполяционное решение на различных сетках, решение на основе контроля погрешности, сравнивать

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

Рисунок 2. Программа для тестирования численных методов решения ОДУ

В ней реализовано около 130 численных методов.

Рисунок 3. UML диаграмма классов численных методов решения ОДУ

Все тесты проводились при повышенной вычислительной точности (38-39 цифр).

3.1. Выбранные методы

Коэффициенты всех методов были уточнены до 30 и более значащих десятичных цифр, и все

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

выбраны следующие методы:

 Multirate (4,2) — непрерывный на основе 4-х стадийного с 2 вычислениями правой части метода

Розенброка, порядок 4 (Новиков Е.А.[6]);

 MRCROS2_1 — непрерывный на основе двухстадийного L1-устойчивого метода Лимонова А. Г.

[7] с комплексными коэффициентами 4-го порядка;

 MRCROS2_Optimal — непрерывный на основе двухстадийного оптимально затухающего метода

Ширкова П.Д. [8] из статьи 1992 года с комплексными коэффициентами, порядок 4(3);

 MRRos — авторский непрерывный жёстко точный 4-х стадийный A-устойчивый метод

(Власов С.А.), порядок 4(3);

 MRRODAS — непрерывный жёстко точный 6-ти стадийный вложенный A-устойчивый метод

Розенброка, порядок 4(3) (V.Savcenco [5]);

 RODAS5 — жёстко точный 8-ми стадийный вложенный A-устойчивый метод Розенброка, порядок

5(4) (Giovanna A. Di Marzo [10]), с интерполяцией Ньютона 4-й степени по пяти ближайшим

точкам 5-го порядка;

 GearAdapt — многошаговый метод Гира [11] на основе формул дифференцирования назад с

разгоном от 1-го до 4-го порядка, с интерполяцией Эрмита 3-й степени по двум ближайшим

точкам 4-го порядка;

 SkvorcovAdapt — явный многошаговый адаптивный метод с разгоном от 1-го до 5-го порядка

(Скворцов Л.М. [12]), с интерполяцией Эрмита 5-й степени по трем ближайшим точкам 6-го

порядка;

 ESDIRK86 — классический метод, основанный на 9-ти стадийном однократно диагонально-

неявном методе Рунге-Кутты из статьи Скворцова Л.М. [13] 6-го порядка и интерполяции

Ньютона 5-й степени по шести ближайшим точкам 6-го порядка или интерполяции Эрмита 5-й

степени по трем ближайшим точкам 6-го порядка;

 Lobatto IIIC8 — классический метод, основанный на полностью неявном 5-ти стадийном методе

Рунге-Кутты Lobatto IIIC8 [9] и интерполяции Ньютона 7-й степени по восьми ближайшим точкам

8-го порядка или интерполяции Эрмита 7-й степени по четырем ближайшим точкам 8-го порядка;

 RADAU IIA15 — классический метод, основанный на полностью неявном 8-ми стадийном методе

Рунге-Кутты RADAU IIA15 [9] и интерполяции Эрмита 15-й степени по восьми ближайшим

точкам 16-го порядка;

 Spline7 — авторский непрерывный коллокационный A-устойчивый метод (Власов С.А. [14]),

сплайн 7-й степени, который имеет 8-й порядок;

 RIIA_B2 — линейно неявная схема LN-эквивалентная неявному методу RadauIIA (Семейство Б)

3-го порядка (Зубанов А.М. [15]), с интерполяцией Ньютона 2-й степени по трем ближайшим

точкам 3-го порядка;

 BORK4 — L4-устойчивый обратный оптимальный неявный метод Рунге-Кутты из диссертации

Пошивайло И.П. [16], порядок 4, с интерполяцией Эрмита 3-й степени по двум ближайшим

точкам 4-го порядка;

3.2. Стратегия выбора шага

Размер шага выбирается следующим образом:

kof

err

Tol





10,

, (10)

где

1i

h — предыдущий шаг,

1i

err — абсолютная погрешность на предыдущем шаге, Tol — задаваемая

погрешность, p — порядок метода, 9,0



kof — коэффициент, который в случае интерполяции решения

позволяет дополнительно дробить шаг для методов высокого порядка.

Далее оценивается текущая погрешность либо по правилу Рунге, либо способом, предусмотренным

решаемым методом, и проверяется условие:

Tolerr



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

шага до тех пор, пока условие не выполнится, в соответствии с формулой:

kof

err

Tol

jiji





1,,

, (11)

где j — номер попытки.

В случае, когда решение интерполируется и используется метод порядка больше седьмого, имеет смысл

дополнительно дробить шаг на число равное количеству интерполируемых точек, с помощью коэффициента

kof в формулах выше.

3.3. Задача Капса

Это автономная сверхжёсткая (из-за параметра) тестовая задача [17], содержащая две компоненты

решения и не имеющая запаздываний:

















],0(,

2212

211

xxxx

xxx





(12)

где

10





Начальная задача Коши:









1)0(

. (13)

Задача Капса имеет плавное решение x

(t) = exp(-2t), x

(t) = exp(-t), не зависящее от параметра жёсткости

μ (собственные значения якобиана при больших μ примерно равны – μ, –1).

В [17] отмечается, что данная задача — это пример сингулярно возмущенной задачи Коши.

Несмотря на то, что задача Коши не имеет запаздываний, было принято решение интерполировать

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

Следует так же заметить, что матрица Якоби для методов Розенброка считается аналитически, то есть точно.

В качестве оценки погрешности методов была выбрана норма C: для разностей приближённого решения и

точного решения (для двух компонент) в точках на всем отрезке интегрирования (от 0 до 1) с постоянным

интервалом (0,01). Количество шагов у всех методов ограниченно одним миллионом. Результаты

тестирования показаны на рисунках 4, 5, 6, 7.

Рисунок 4. Достигнутая точность — трудоёмкость.

Точность на рисунках отображает количество значащих цифр, а трудоемкость считается как





time

log ,

где time — время выполнения программы в миллисекундах.

BORK4 не показан на рисунках, так как он не справился с задачей, не смотря на то, что используется

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

нелинейных функций системы, которую приходится решать на итерациях Ньютона. Данный метод работает

на менее жёсткой задаче с параметром

10



. При текущей же жёсткости работает BORK3 и BORK2 [16],

но они имеют 3-й и 2-й порядок соответственно.

Неожиданный результат получился у RODAS5: до точности в 15 значащих цифр он имеет 5-й порядок, а

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

Свою работоспособность показал авторский непрерывный жёстко точный 4-х стадийный A-устойчивый

метод Розенброка 4-го порядка для разностного решения и 3-го порядка для интерполяционного.

По результатам тестирования все методы можно разбить на 3 класса точности: те, с помощью которых

можно получить решение с инженерной точностью (5-6 значащих цифр), те, которые рассчитаны на

получение решения высокой точности (10-12 значащих цифр), и те, которые могут считать со сверх высокой

точностью (20-30 значащих цифр). Разбиение выглядит следующим образом:

1. Multirate (4,2), MRCROS2_1, MRCROS2_Optimal, MRRos, RIIA_B2;

2. MRRODAS, GearAdapt, RODAS5, ESDIRK86;

3. SkvorcovAdapt, Lobatto IIIC8, RADAU IIA15, Spline7.

В классе инженерной точности наилучшее соотношение ( lg(||Tol||) - lg(||error||)) трудно выявить.

В классе высокой точности наилучшее соотношение ( lg(||Tol||) - lg(||error||)) у RODAS5.

В классе сверх высокой точности наилучшее соотношение ( lg(||Tol||) - lg(||error||)) у SkvorcovAdapt, он же

и самый точный на данной задаче.

Рисунок 5. Задаваемая точность - достигаемая точность для класса инженерной точности.

Рисунок 6. Задаваемая точность - достигаемая точность для класса высокой точности.

Рисунок 7. Задаваемая точность - достигаемая точность для класса сверх высокой точности.

3.4. Пример Крайсса

Это неавтономная сверхжёсткая (из-за параметра) тестовая задача [17], содержащая две компоненты

решения и не имеющая запаздываний:

       

,30,10,





















ttytEtEty





(14)

где

-12

10



 









   

















cossin

sincos

Начальная задача Коши:









3)0(

1)0(

. (15)

Имеет плавное решение, которое почти не зависит от параметра, отвечающего за жёсткость.

Задача диссипативна в евклидовой норме, так как матрица E(t) является ортогональной при всех t. Точнее,

задача принадлежит и классу F

v=-1

, где v=-1 — постоянная в одностороннем условии Липшица.

Ван Вельдхёйзен ввел класс задач P, к которому можно отнести и данную задачу [17]. Эти задачи

допускают одновременное вхождение в решение медленных и быстрых компонент, содержат малый параметр,

позволяющий перейти к произвольно высокой степени жёсткости, и имеют зависящие от времени

собственные векторы. Вводится данный класс вместе с определением о D-устойчивости.

Эталонное решение было построено явным многошаговым адаптивным методом Скворцова с контролем

погрешности



eps . Внешний вид решения показан на рисунке 8.

Рисунок 8. Эталонное решение примера Крайсса.

Функции компонент решения, получаемые в разностном виде, были интерполированы, чтобы посмотреть

точность интерполяции. Матрица Якоби для методов Розенброка считается аналитически, то есть точно. Так

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

автономности. В качестве оценки погрешности методов была выбрана норма C: для разностей приближённого

решения и полученного эталонного решения (для двух компонент) в точках на всем отрезке интегрирования

(от 0 до 3) с постоянным интервалом (0,01). Количество шагов у всех методов ограниченно одним миллионом.

Результаты тестирования показаны на рисунках 9, 10, 11, 12.

Рисунок 9. Достигнутая точность — трудоёмкость.

BORK4, не показан на рисунках, так как он не справился с задачей. Но он работает на менее жёсткой

задаче с параметром

-6

10



RIIA_B2, RODAS5 работают, но слишком медленно, поэтому они тоже не показаны на рисунке.

По результатам тестирования все методы можно разбить на 3 класса точности: те, с помощью которых

можно получить решение с низкой точностью (2-3 значащие цифры), те, которые рассчитаны на получение

решения высокой точности (10-12 значащих цифр), и те, которые могут считать со сверх высокой точностью

(20-30 значащих цифр). Разбиение выглядит следующим образом:

1. Multirate (4,2), MRCROS2_1, MRCROS2_Optimal, MRRos, MRRODAS, Spline7.

2. GearAdapt и ESDIRK86.

3. SkvorcovAdapt, Lobatto IIIC8, RADAU IIA15.

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

параметром

-6

10



В классе инженерной точности наилучшее соотношение ( lg(||Tol||) - lg(||error||)) у MRRODAS.

В классе высокой точности наилучшее соотношение ( lg(||Tol||) - lg(||error||)) у ESDIRK86.

В классе сверх высокой точности наилучшее соотношение ( lg(||Tol||) - lg(||error||)) у SkvorcovAdapt, он же

и самый точный на данной задаче.

Рисунок 10. Задаваемая точность - достигаемая точность для класса инженерной точности.

Рисунок 11. Задаваемая точность - достигаемая точность для класса высокой точности.

Рисунок 12. Задаваемая точность - достигаемая точность для класса сверх высокой точности.

3.5. Осциллятор ван дер Поля

Это автономная сверхжёсткая (из-за параметра) тестовая задача, содержащая две компоненты решения и

не имеющая запаздываний. Осциллятор с нелинейным затуханием подчиняется следующей системе

уравнений:

  













],0(,

xxyay



(16)

где

10a ,



Начальная задача Коши:









0)0(

2)0(

. (17)

Решение было построено сплайн методом [14] 7-й степени, порядок 8, с контролем погрешности



eps (Основная сетка: 23 485 узлов).

Ещё два решения были получены на двух сетках, которые сгустили в 2 и в 4 раза, для оценки локальной

погрешности.

Правило Рунге-Ричардсона позволяет численно оценить локальную погрешность и порядок метода:

log



 ,

error



 , (18)

erroruy 

где p— порядок, error — погрешность, u

, u

— решения на сгущающихся сетках, y — эталонное решение.

Результаты экстраполяции показаны на рисунках 13, 14, 15, 16.

Рисунок 13. Эталонное решение (1-я компонента).

Рисунок 14. Эталонное решение (2-я компонента).

Рисунок 15. Эталонное решение (фазовый портрет).

Рисунок 16. Локальная точность (– lg(||error||)).

Функции компонент решения, получаемые в разностном виде, были интерполированы, чтобы посмотреть

точность интерполяции. Матрица Якоби для методов Розенброка считается аналитически, то есть точно. В

качестве оценки погрешности методов была выбрана норма C: для разностей приближённого решения и

полученного эталонного решения (для двух компонент) в точках на всем отрезке интегрирования (от 0 до 2) с

постоянным интервалом (0,01). Количество шагов у всех методов ограничено одним миллионом. Результаты

тестирования показаны на рисунках 17, 18, 19.

Рисунок 17. Достигнутая точность — трудоёмкость.

BORK4 и RIIA_B2 не показаны на рисунках, так как они не справились с задачей. Но они работают на

менее жёсткой задаче ван дер Поля с параметром

10a .

Так как решение содержит быстро изменяющиеся компоненты, GearAdapt работает плохо, поэтому он не

показан на рисунках.

По результатам тестирования все методы можно разбить на 2 класса точности: те, с помощью которых

можно получить решение с инженерной точностью (5-6 значащих цифр), и те, которые рассчитаны на

получение решения высокой точности (10-12 значащих цифр). Разбиение имеет следующий вид:

1. Multirate (4,2), MRCROS2_1, MRCROS2_Optimal, MRRos, MRRODAS, SkvorcovAdapt, ESDIRK86;

2. RODAS5, Lobatto IIIC8, RADAU IIA15, Spline7.

В первом классе наилучшее соотношение ( lg(||Tol||) - lg(||error||)) у MRCROS2_Optimal.

Во втором классе наилучшее соотношение ( lg(||Tol||) - lg(||error||)) у RODAS5.

Уравнение ван дер Поля применяется и в физике, и в биологии. Так, например, в биологии создана

модель ФитцХью-Нагумо. Данное уравнение также было использовано в сейсмологии для моделирования

геологических разломов.

Уравнение ван дер Поля используется в радиотехнике для описания не только триодного генератора, но и

генератора на туннельном диоде или транзисторного генератора [18].

Один из ярких примеров: ван дер Поль (в соавторстве с ван дер Марком) обсуждает эффективность

применения этого уравнения для описания кардиоритмов.

Рисунок 18. Задаваемая точность - достигаемая точность для класса инженерной точности.

Рисунок 19. Задаваемая точность - достигаемая точность для класса высокой точности.

Список литературы

1. Андронов А. В. Математическое моделирование систем, процессов и явлений во временной

области. Машиностроение и смежные отрасли. CAD/CAM/CAE Observer #7 (59) / 2010

//Интернет ресурс http://backend.pa10.rk6.bmstu.ru/pa10/doc/math-mod-pa10.pdf , 09.02.2021.

2. Жолковский Е.К., Белов А.А., Калиткин Н.Н. Решение жестких задач Коши явными схемами с

геометрически-адаптивным выбором шага // Препринты ИПМ им. М.В.Келдыша. 2018. № 227.

20 с. doi:10.20948/prepr-2018-227 //Интернет ресурс

http://library.keldysh.ru/preprint.asp?id=2018-227 , 09.02.2021.

3. The doubledouble homepage//Интернет ресурс

https://boutell.com/fracster-src/doubledouble/doubledouble.html , 23.01.2018.

4. David H. Bailey, Xiaoye S. L., Yozo Hida. Quad-Double Arithmetic: Algorithms, Implementation, and

Application: October 30, 2000// Интернет ресурс

www.davidhbailey.com/dhbpapers/quad-double.pdf , 23.01.2018.

5. Savcenco V. Construction of a multirate RODAS method for stiff ODEs: Journal of Computational

and Applied Mathematics Volume 225, Issue 2, 15 March 2009, Pages 323-337.//Интернет ресурс

http://www.sciencedirect.com/science/article/pii/S0377042708003701, 29.12.2017.

6. Новиков E.A. Исследование (m,2)-методов решения жёстких систем — Институт

вычислительного моделирования СО РАН, Красноярск, Россия: Ж. Вычислительные

технологии, том 12, № 5, 2007.//Интернет ресурс

http://cyberleninka.ru/article/n/issledovanie-m-2-metodov-resheniya-zhestkih-sistem.pdf , 29.12.2017.

7. Лимонов А.Г. Разработка двухстадийных схем Розенброка с комплексными коэффициентами и

их применение в задачах моделирования образования периодических наноструктур —

Автореферат диссертации на соискание ученой степени кандидата физико-математических

наук: Уральский государственный университет им. А.М. Горького — Екатеринбург,

2010.//Интернет ресурс

http://elar.urfu.ru/bitstream/10995/3112/2/urgu0811s.pdf ,29.12.2017.

8. Ширков П.Д. Оптимально затухающие схемы с комплексными коэффициентами для жёстких

систем ОДУ — Институт математического моделирования РАН — Москва: Матем.

моделирование, 1992, том 4, номер 8, 47–57.//Интернет ресурс http://mi.mathnet.ru/mm2101 ,

29.12.2017.

9. Хайрер Э., Ваннер Г. Решение обыкновенных дифференциальных уравнений. Жёсткие и

дифференциально-алгебраические задачи. — М.: Мир, 1999.

10. Giovanna A. Di Marzo RODAS5(4) Méthodes de Rosenbrock d’ordre 5(4) adaptées aux problémes

différentiels-algébriques — Mémoire de diplôme en Mathématiques: Université de Genéve Faculté

des Sciences Section de Mathématiques — Genéve, Mars 1993.//Интернет ресурс

http://cui.unige.ch/~dimarzo/papers/DIPL93.pdf 29.12.2017.

11. Хайрер Э., Нёрсетт С., Ваннер Г. Решение обыкновенных дифференциальных уравнений.

Нежёсткие задачи. — М.: Мир, 1990.

12. Скворцов Л. М. Диагонально-неявные методы Рунге–Кутты для жёстких задач: Ж. вычисл.

матем. и матем. физ., 2006, том 46, номер 12, С. 2209–2222.//Интернет ресурс

http://mi.mathnet.ru/zvmmf367 , 29.12.2017.

13. Скворцов Л. М., Явный многошаговый метод численного решения жёстких дифференциальных

уравнений. — вычисл. матем. и матем. физ., 2007, том 47, номер 6, 959–967//Интернет ресурс

http://mi.mathnet.ru/zvmmf4593 , 22.02.2018.

14. Власов С.А., Ширков П.Д. Сплайн-интегрирование дифференциальных уравнений с

запаздыванием: Системный анализ в науке и образовании. Университет "Дубна", eISSN: 2071-

9612. Дубна 2010.

15. Зубанов А. М., Кутрухин Н. Н., Ширков П. Д. О построении линейно неявных схем, LN-

эквивалентных неявным методам Рунге–Кутты: Компьютерные исследования и моделирование

2012 Т. 4 № 3 С. 483-496.//Интернет ресурс

http://crm.ics.org.ru/uploads/crmissues/crm_2012_3/483-496.pdf , 29.12.2017.

16. Пошивайло И.П. Жёсткие и плохо обусловленные нелинейные модели и методы их расчета —

Диссертация на соискание ученой степени кандидата физико-математических наук: Институт

прикладной математики им. М.В. Келдыша — Mосква, 2014.//Интернет ресурс

http://keldysh.ru/council/3/D00202403/poshivaylo_diss.pdf , 29.12.2017.

17. Деккер К., Вервер Я. Устойчивость методов Рунге-Кутты для жёстких нелинейных

дифференциальных уравнений. — М.: Мир, 1988.

18. Кузнецов А. П., Селиверстова Е. С., Трубецков Д. И., Тюрюкина Л. В. Феномен уравнения ван

дер Поля. — Обзор актуальных проблем нелинейной динамики: Изв. вузов «ПНД», т. 22, № 4,

2014. УДК 517.91, 517.938, 51.73 // Интернет ресурс http://sgtnd.narod.ru/papers/2014PND3.pdf ,

26.03.2018.