Параметризация задачи Коши для нелинейных дифференциальных уравнений с контрастными структурами

 ных задач Коши, но не могут сравниться по быстродействию с явными методами, т. к. на каждом шаге неявного метода приходится решать систему нелинейных уравнений. Этот недостаток удается устранить при помощи метода продолжения решения33, который заключается в замене исходного аргумента задачи Коши на новый. Можно выбирать различные аргументы продолжения решения, но наиболее часто используют наилучший аргумент. Данный аргумент отсчитывается по касательной к интегральной кривой рассматриваемой задачи и обладает рядом исключительных свойств. Для задачи (3)-(4) наилучший аргумент X записывается в скалярной форме:

dX¹ = dy 1 +-----+ dy 2 + dt ².      (9)

Будем полагать, что все переменные у1,у2,...,yn и аргумент t зависят от X. Дополняя систему (3) соотношением (9) и разрешая полученную систему относительно производных по λ, получим преобразованную систему dyi _ fi (t,У1,^,Уп)

d^   Q(t,Ух,-,Уп) ’ dt _        1

d^   Q (t, У1,^, Уп) ’ i _ 1,..., n,

где

_____________ Q ( t , У 1 , ^ , y „ )= ________________ = 7^{1 +} ft ( t , yv-, У п ) + - ^{+ f} 2 ( t , yv-, У п ) .

Начальные условия (4) для системы (10) перепишутся в виде

Задача (10)-(11) обладает рядом важных для численного решения свойств34:

• 1.    Квадратичная норма правой части системы (10) равна единице, т. е. устраняются все вычислительные трудности, связанные с неограниченным возрастанием правых частей системы (3).

• 2.    Обусловленность системы (10) является наилучшей.

• 3.    Показатель жесткости преобразованной системы (10) меньше, чем у исходной.

Все отмеченные свойства дают возможность решать преобразованную задачу (10)-(11) любыми численными методами, в том числе и явными.

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

Результаты исследования

В статье А. А. Белова и Н. Н. Калит-кина [3] для проверки эффективности численных методов при решении задач с контрастными структурами предложен ряд тестовых задач. Проанализируем решение одной из этих задач, называемой степенным тестом.

Постановка задачи

Рассмотрим начальную задачу для уравнения следующего вида [3]:

du __ £(t)(u2 -a2)2

dt        ( u ² + a ² )

с начальным условием u (0) = 0.

При ^ ( t ) = ^ ₀ • cos t решение задачи (12)-(13) имеет вид:

u

1 + ^/1 + 4 a ^E²( t)

Z-(⁰⁾= Ую, t⁽⁰⁾= 0, i = 1, —,ⁿ. (11) где H(t) = ^_osint.

В дальнейших расчетах a = п. Жесткость задачи характеризуется множителем ^₀ при периодической функции. Условно по значению параметра £₀ можно разбить задачи вида (12)-(13) на следующие классы [5]: при малых значениях ^₀ (порядка нескольких единиц) задача является нежесткой, при ξ₀ ≥ 10 – жесткой и при ξ₀ ≥ 1000 – сверхжесткой. Это определяется видом правой части: при приближении значения косинуса к нулю скорость изменения решения, даже при большом параметре £₀, стремится к нулю, а при стремлении значения косинуса к единице значение правой части становится наибольшим. Таким образом, происходит смена плавно меняющихся компонент решения с быстро меняющимися переходными зонами, которые тем короче, чем больше значение ^₀. Данные особенности видны на рис. 1.

Явные методы решения

Результаты решения задачи явным методом Эйлера и явным методом Рунге-Кутты четвертого порядка точности с постоянным шагом, схемы которых

MORDOVIA UNIVERSITY BULLETIN даны выше, представлены на рис. 1 и в табл. 1.

На рис. 1 изображены кривые, соответствующие решениям задачи (12)-(13) для значений ^₀ = 1; 10; 100, полученные методом Рунге-Кутты четвертого порядка точности с постоянным шагом h = 0,001. Полученные численные решения визуально совпадают с аналитическими. Для явного метода Эйлера кривые аналогичны изображенным на рис. 1.

В табл. 1 приведены данные о средней погрешности е, вычисленной с использованием аналитического решения (14), и времени счета t_c решений задачи (12)-(13), полученных явными методами с постоянным шагом интегрирования. Из нее видно, что метод Рунге-Кутты четвертого порядка выигрывает у явного метода Эйлера в точности, хотя и требует большего времени счета. Стоит также отметить, что при повышении значения параметра $₀ невозможно получить решение при малых шагах интегрирования. Так, для $₀ = 10 невозможно получить решение при

Р и с. 1. Численное решение задачи (12)-(13) при $₀ = 1; 10; 100, метод Рунге-Кутты четвертого порядка с постоянным шагом h = 10^-3

F i g. 1. Numeric solution of problem (12)-(13) for $₀ = 1; 10; 100, fourth order Runge-Kutta method with constant step size h = 10^-3

Т а б л и ц а 1

T a b l e 1

Погрешность решения и время счета для задачи (12)-(13), явные методы с постоянным шагом

Computational error and computational time for problem (12)-(13), explicit methods with constant step size

Параметры / Parameters Явный метод Эйлера / Explicit Euler method Метод Рунге-Кутты четвертого порядка / Fourth order Runge-Kutta method $0 h ɛ tc, c ɛ tc, c 1 0,1 0,4548 0,009 2,1332∙10-4 0,0081 0,01 0,0361 0,0095 2,2776∙10-8 0,0656 0,001 0,0035 0,0255 2,1461∙10-12 0,1732 10 0,1 – – – – 0,01 0,0826 0,0109 2,3893∙10-5 0,0168 0,001 0,008 0,0301 2,5136∙10-9 0,0759 100 0,1 – – – – 0,01 – – – – 0,001 0,0127 0,0381 2,4843∙10-6 0,1011 h = 0,1, а при $0 = 1 000 нельзя получить решение при h ≥ 10-3.

Здесь и далее расчет проведен на персональном компьютере с процессором DualCore Intel Core i5-660 3477 MHz, оперативной памятью 2,00 ГБ DDR3 SDRAM Kingston «ValueRAM» KVR1333D3N9/2G (2 шт.), видеокартой ATI Radeon HD 5770 (Juniper) 1,00 ГБ и 64-разрядной операционной системой Windows 10. Используемые методы реализованы в вычислительной среде Matlab R2012b.

Замечание 3

В рассматриваемом диапазоне шагов интегрирования не удалось получить решение задачи (12)-(13) для $₀ > 1 000.

Анализируя полученные результаты, отметим, что для решения задачи с фиксированным $₀ шаг интегрирования h не может превышать величину 1/$₀. Это говорит о малоэффективности применения явных методов для решения сверхжестких задач при больших значениях $₀.

Неявный метод Эйлера

Явные методы с постоянным шагом интегрирования позволяют получить ре-498

шение задачи (12)-(13) только при умеренных значениях $₀. Уже при $₀ = 10³ в рассматриваемом диапазоне шагов интегрирования не удается найти решение. Можно ожидать, что неявные методы за счет расширения области устойчивости позволят получить решение при больших шагах интегрирования.

Используем неявный метод Эйлера. На k-ом шаге решение задачи (12)-(13) неявным методом Эйлера находится как решение нелинейного уравнения вида:

uk+1 = uk^{- h} •

£oco^s(tk+1) ((uk+1 )^{2 -}a²)

(uk+1 )²⁺a²

. (15)

Уравнение (15) решается методом простых итераций и методом Ньютона. Итерационный процесс метода простых итераций задается соотношением:

uk +11) = uk — ^h •

(j)\²       2 2

^0cos(tk+1) ((uj ) — a )

⁽j ) uk +1 )

,

а для метода Ньютона ‒

u

(j+1) k+1

_ u (j) - ^{F (}tk+1, uk+1)

^k+1F4t    u⁽ j))

Fu⁽tk+1, uk+1)

.

В формуле (16)

^₀cost (u²- a²)

F (t, u) = u - uk + h------—— u2 + a2

FU( t, u) = 1- h •

2^₀cost (u²- a²)u (u²+ 3a²) (u² + a² )2

dF (t, u) где F,'(t, u) = —^.

du

Кривые, полученные при решении задачи (12)-(13) неявным методом Эйлера, имеют вид, аналогичный представленному на рис. 1. В табл. 2 даны значения средней погрешности е и времени счета tc решений задачи (12)-(13), полученные неявным методом Эйлера с постоянным шагом интегрирования при использовании метода простых итераций и Ньютона.

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

Т а б л и ц а 2

T a b l e 2

Погрешность решения и время счета для задачи (12)-(13), неявный метод Эйлера с постоянным шагом

Computational error and computational time for problem (12)-(13), implicit Euler method with constant step size

Параметры / Parameters Неявный метод Эйлера / Implicit Euler method Метод простой итерации / Fixed-point iteration method Метод Ньютона / The Newton’s method 5о h ɛ tc, c ɛ tc, c 1 0,1 0,3137 0,0662 0,3138 0,0255 0,01 0,034 0,1774 0,0345 0,1457 0,001 8,6629∙10-4 1,188 0,0035 0,6134 10 0,1 – – – – 0,01 0,074 0,1864 0,0787 0,0878 0,001 0,0029 1,1653 0,008 0,7302 100 0,1 – – – – 0,01 – – – – 0,001 0,0078 0,9977 0,0126 1,5954 при этом большее время счета. Для повышения точности можно, как и в статье А. А. Белова и Н. Н. Калиткина «Численные методы решения задач Коши с контрастными структурами» [5], рекомендовать использование более сложных методов, например, неявных обратных схем Рунге-Кутты [37] или явно-неявных комплексных схем Розенброка [38].

Наилучшая параметризация

Из приведенных выше расчетных данных видно, что ни явные методы, ни неявный метод Эйлера не позволяют получить решение задачи (12)-(13) без значительных вычислительных затруднений. Применим к этой задаче метод продолжения решения по наилучшему аргументу35, описанный выше.

Наилучший аргумент X для рассматриваемой задачи запишется в скалярном виде:

dX¹ = du²+ dt 2.        (17)

Преобразуя к аргументу (17) систему (12), получим:

du dn dt d^

-§(t) (u²- a²)²

( u² + a²) + £²( t)( u² - a²)

( u²+ a²)

( u²+ a²) + ^²( t)( u²- a²)

Начальные условия (13) для системы (18) примут вид:

u(0) = 0, t(0) = 0.      (19)

Как отмечалось выше, преобразованная задача (18)-(19) обладает рядом преимуществ при численном решении. В табл. 3 приведены результаты расчета преобразованной задачи (18)-(19) методом Рунге-Кутты четвертого по- рядка точности с постоянным шагом (качественно для метода Эйлера результаты аналогичны), где l – шаг по аргументу λ. Кривые, полученные при численном решении, аналогичны приведенным на рис. 1. Для сравнения в табл. 3 приведены результаты решения исходной задачи (12)-(13) с шагом, подбираемым из условия равенства порядков средней ошибки решения исходной и преобразованной задач.

Видно, что для непреобразованной задачи (12)-(13) для достижения того же порядка ошибки приходится брать шаг до нескольких порядков меньше по сравнению с преобразованной. При этом, несмотря на то, что интервал изменения значений аргумента λ больше, чем у аргумента t, для задачи (18)-(19) удается сократить время счета от 0,4 (при ^0 = 1) до 10 раз (при $₀ = 1 000).

В табл. 4 приведены результаты расчета задач (12)-(13) и (18)-(19) методом Рунге-Кутты четвертого порядка с переменным шагом. Процедура изменения шага описана выше. Для преобразованной задачи (18)-(19) параметр точности для правила Рунге 9 = 10^-12. Для исходной задачи (12)-(13) параметр 6 выбирался из условия равенства порядков средней погрешности для обеих задач.

Результаты показывают, что и при переменном шаге интегрирования использование наилучшего аргумента X позволяет получить вычислительные преимущества: удается сократить время счета минимум в два раза. Для сверхжестких задач уменьшение времени счета может достигать нескольких порядков. Отметим также, что, в отличие от преобразованной задачи, решение исходной задачи зависит от начального шага; в особенности это характерно для больших значений $₀.

Однако стоит указать и на ряд недостатков наилучшей параметризации. Преобразованная задача (18)-(19)

34 Там же.

Т а б л и ц а 3

T a b l e 3

Средняя ошибка и время счета задач (12)-(13) и (18)-(19), метод Рунге-Кутты четвертого порядка с постоянным шагом

Computational error and computational time for problems (12)-(13) and (18)-(19), fourth order Runge-Kutta method with constant step size

ξ0	Параметризованная задача / Parameterized problem			Непараметризованная задача / Nonparametrized problem
	l	ɛ	tc, c	h	ɛ	tc, c
1	0,1	5,4369∙10-7	0,0213	0,02	3,6307∙10-7	0,0418
	0,01	2,9799∙10-11	0,183	0,002	3,6427∙10-11	0,2405
	0,001	5,0522∙10-14	0,5505	0,0002	6,4592∙10-13	2,0987
10	0,1	2,4647∙10-4	0,0236	0,02	3,0375∙10-4	0,0422
	0,01	3,2723∙10-9	0,0969	0,001	2,5136∙10-9	0,4414
	0,001	3,7533∙10-12	0,5886	2∙10-4	3,662∙10-12	2,2979
100	0,1	0,3525	0,0289	0,004	0,0012	0,1303
	0,01	4,834∙10-6	0,1178	0,001	2,4843∙10-6	0,4484
	0,001	5,0063∙10-11	0,5289	5∙10-5	1,7076∙10-11	8,2501
1000	0,1	–	–	–	–	–
	0,01	0,0047	0,2086	4,7∙10-4	0,0012	0,9134
	0,001	4,0885∙10-8	0,5288	7∙10-5	5,9638∙10-8	5,9466

Т а б л и ц а 4

T a b l e 4

Средняя ошибка и время счета задач (12)-(13) и (18)-(19), метод Рунге-Кутты четвертого порядка с переменным шагом

Computational error and computational time for problems (12)-(13) and (18)-(19), fourth order Runge-Kutta method with variable step size

Параметры / Parameters Параметризованная задача / Parameterized problem Непараметризованная задача / Nonparametrized problem ξ0 h0 / l0 θ ɛ tc, c θ ɛ tc, c 1 0,1 10-12 1,1561∙10-10 0,2388 10-10 2,1774∙10-10 0,4892 0,01 10-12 1,6079∙10-10 0,2592 3∙10-10 1,1338∙10-10 0,571 0,001 10-12 1,7008∙10-10 0,1889 10-10 2,1458∙10-12 1,3038 10 0,1 10-12 2,365∙10-9 0,2065 10-7 1,5164∙10-9 1,6533 0,01 10-12 2,2218∙10-9 0,2812 4∙10-8 1,1642∙10-9 2,025 0,001 10-12 1,6625∙10-9 0,2691 10-7 2,5134∙10-9 1,3089 100 0,1 10-12 3,497∙10-8 0,2982 10-4 5,8346∙10-8 3,2128 0,01 10-12 3,1933∙10-8 0,3346 10-5 2,3953∙10-8 4,0536 0,001 10-12 3,0792∙10-8 0,3015 3∙10-5 2,3106∙10-8 4,7691 1000 0,1 10-12 4,2953∙10-7 0,2729 – – – 0,01 10-12 3,2001∙10-7 0,2318 10-3 9,1802∙10-8 15,794 0,001 10-12 8,7845∙10-7 0,1772 10-2 6,4254∙10-7 9,8967 имеет более сложный вид и большую размерность по сравнению с исходной (12)-(13). Для рассматриваемой задачи это несущественно, но при решении задач большой размерности преимущества, даваемые наилучшей параметризацией, могут перекрываться сложностью преобразованной задачи.

Еще одним недостатком является заниженная оценка погрешности по правилу Рунге для преобразованной задачи. Уже при $0 = 1 оценка погрешности на два порядка ниже, чем реальная средняя погрешность. Для сверхжестких задач разница возрастает до пяти порядков. Диаметрально противоположная ситуация характерна для исходной задачи. Здесь правило Рунге дает завышенную оценку погрешности начиная с $0 = 10. Если посмотреть на характер изменения шага интегрирования, изображенный на рис. 2 в полулогарифмических координатах, то видно, что уже при значении $0 = 30 значение шага интегрирования не поднимается выше 10-5, снижаясь до 10-12 в окрестности пограничных и внутреннего слоев. При этом значение шага изме- няется плавно, т. е. после прохождения внутреннего слоя шаг начинает плавно возрастать до максимального значения, а при подходе к внутреннему слою – плавно убывать до минимума (при этом ширина переходного слоя уменьшается при увеличении ξ0 ). Данный механизм изменения шага интегрирования и приводит к уменьшению погрешности решения на горизонтальных участках до значений меньших, чем задаваемая точность θ.

Для преобразованной задачи кривая решения в зависимости от аргумента λ изменяется более плавно, как это отмечено в статьях А. А. Белова и Н. Н. Ка-литкина [3; 5]. По этой причине такого же падения погрешности не происходит; наоборот, в узловых точках, соответствующих внутренним слоям, концентрируется максимальная погрешность. В заключение отметим, что, исходя из характера изменения шага интегрирования, можно объяснить и увеличение времени счета для исходной задачи. Помимо усложнения процесса решения, связанного со сменой шага, уже при относительно малом значении $₀ = 30 шаг

Р и с. 2. Шаг интегрирования в полулогарифмических координатах для задачи (12)-(13) при $₀ = 30, метод Рунге-Кутты четвертого порядка, 9 = 10-4

F i g. 2. Step size in semilogarithmic coordinates for problem (12)-(13),

$₀ = 30, fourth order Runge-Kutta method, 9 = 10-4

интегрирования не поднимается выше 10-5, при повышении же значения $₀ шаг интегрирования будет уменьшаться еще сильнее, а вместе с этим будет возрастать время счета.

Замечание 4

В отличие от непреобразованной задачи (12)-(13), переход к наилучшему аргументу позволяет получить хорошие результаты при параметре жесткости $₀ большем, чем рассмотренные в табл. 3 и табл. 4. При задаваемой точности θ = 10^-12 и начальном шаге h = 0,001 удается построить решение вплоть до $₀ = 106.

Замечание 5

Использование неявного метода Эйлера с переменным шагом не позволило получить результаты значительно лучшие, чем приведенные в табл. 2. Получить решение со средней погрешностью е порядка 10^-2 удается, но время счета по сравнению с результатами табл. 2 возрастает до 10 раз при $₀< 100. При большей точности (ɛ ≤ 10^-3) или большей жесткости ($₀ > 1000) результаты либо не удавалось получить, либо время счета составляло до 1 000 с.

Обсуждение и заключение

В статье рассмотрено применение традиционных явных и неявных методов решения задачи Коши с контрастными структурами. Получены численные решения задачи (12)-(13) с постоянным и переменным шагом, изменяемым по правилу Рунге. Анализируя результаты, приведенные на рис. 1 и в табл. 1‒2, отметим следующее.

• 1.    Все полученные расчетные данные хорошо согласуются как с аналитическим решением, так и с результатами, полученными А. А. Беловым и Н. Н. Калиткиным [3; 5].

• 2.    При малых значениях параметра жесткости $₀ явные методы с постоянным шагом позволяют получить решение с приемлемой погрешностью, но при больших $₀ явные методы становятся малоэффективными, т. к. при фиксированном $₀ шаг интегрирования

• 3.    Используемый в статье неявный метод Эйлера дает возможность уменьшить погрешность полученного решения для малых $₀, но не позволяет построить решение при $₀ > 1 000. При этом, в силу особенностей реализации, время счета для неявного метода Эйлера значительно превышает аналогичные значения для явных методов.

• 4.    Метод Ньютона, используемый при решении нелинейных уравнений в неявном методе Эйлера, при малых $₀ позволяет сократить время счета по сравнению с методом простых итераций, но уступает ему в точности.

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

h не может превосходить значения 1/$₀. Уже при $₀ > 1 000 решение в рассматриваемом диапазоне шагов интегрирования построить не удалось.

Для устранения отмеченных недостатков исходная задача (12)-(13) преобразуется к наилучшему аргументу λ. Преобразованная задача (18)-(19) обладает рядом вычислительных преимуществ (см. табл. 3‒4).

• 1.    При решении задачи (18)-(19) с постоянным шагом удается получить большую точность. При этом для исходной задачи, в зависимости от значения $₀, аналогичные результаты можно получить с меньшим (от нескольких раз до нескольких порядков) шагом. Время счета для преобразованной задачи меньше от 0,4 (при $₀ = 1) до 10 раз (при $₀ = 1 000).

• 2.    При использовании переменного шага интегрирования задача (18)-(19) также выигрывает в быстродействии. Переход к аргументу λ позволяет при одинаковых значениях средней погрешности сократить время счета от двух раз для малых значений $₀ до нескольких порядков при $₀ = 1 000. Более того, в отличие от исходной, преобразованную задачу можно решать вплоть до значений $₀ = 10⁶.

• 3.    Для преобразованной задачи (18)-(19) процедура смены шага слабо зависит от начального шага 1₀, в то время как для исходной задачи решение существенно зависит от выбора начального шага h₀.

Помимо значительных преимуществ задача (18)-(19) имеет и ряд недостатков.

• 1.    Правило Рунге для преобразованной задачи дает заниженную оценку локальной погрешности. Уже при $₀ = 1 оценка погрешности на два порядка ниже, чем реальная средняя погрешность. Для сверхжестких задач разница возрастает до пяти порядков. Диаметрально противоположная ситуация отмечается с исходной задачей. Для нее правило Рунге дает завышенную оценку погрешности начиная с ξ = 10. Это необходимо учитывать при построении решения.

• 2.    В отличие от исходного уравнения (12), преобразованная система (18) имеет большую размерность и усложненный вид. В ряде случаев это может

  Том 28, № 4. 2018 перекрывать преимущества наилучшей параметризации.

Несмотря на отмеченные недостатки, наилучшая параметризация является одним из наиболее эффективных методов решения плохо обусловленных задач, к которым можно отнести и задачи с контрастными структурами. Показано, что традиционные явные и неявные методы имеют существенные ограничения в использовании при больших $₀ и малоэффективны. В то же время переход к наилучшему аргументу позволяет получить решение рассматриваемой задачи даже при больших $₀ быстро и с приемлемой точностью. В дальнейшем для подтверждения полученных результатов необходимо рассмотреть решение еще ряда тестовых и прикладных задач данного класса. Однако, опираясь на уже полученные результаты решения различных задач^36-38 [33], можно предположить, что метод продолжения решения по наилучшему аргументу будет также эффективен и для других задач данного класса.

Поступила 05.07.2018; принята к публикации 03.09.2018; опубликована онлайн 28.12.2018

Об авторах:

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