ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ, 2021, том 57, № 7, с.976-987

ЧИСЛЕННЫЕ МЕТОДЫ

УДК 519.622.2

АЛГОРИТМ АДАПТИВНОЙ ИНТЕРПОЛЯЦИИ

НА РАЗРЕЖЕННЫХ СЕТКАХ ДЛЯ ЧИСЛЕННОГО

ИНТЕГРИРОВАНИЯ СИСТЕМ ОБЫКНОВЕННЫХ

ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ

С ИНТЕРВАЛЬНЫМИ НЕОПРЕДЕЛЁННОСТЯМИ

Рассматриваются теоретические аспекты обобщения алгоритма адаптивной интерполяции

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

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

строении адаптивного иерархического разбиения области неопределённости на подобласти,

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

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

с заданной точностью. Основным недостатком алгоритма является его экспоненциальная

зависимость от количества интервальных параметров. Уже при наличии пяти-шести ин-

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

этим предлагается вместо интерполяции на регулярной сетке использовать интерполяцию

на адаптивных разреженных сетках, что позволяет в ряде случаев существенно расширить

область применения алгоритма. Получена оценка глобальной погрешности алгоритма на

разреженных сетках. Показана линейная зависимость глобальной погрешности алгорит-

ма от уровня сетки. На представительном ряде задач продемонстрирована эффективность

предложенного подхода.

DOI: 10.31857/S0374064121070104

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

точной информации о параметрах изучаемого явления или процесса [1]. По этой причине

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

уравнений (ОДУ) с интервальными параметрами и интервальными начальными условиями.

Решением таких систем является интервальная оценка значений фазовых переменных в зави-

симости от времени. Существует несколько подходов к решению данного класса задач: одни

основаны на использовании классической интервальной арифметики [2-4], другие - на ап-

проксимации множества решений с помощью геометрических примитивов [5-7], а в основе

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

параметров [8-11]. В общем случае наблюдается некоторый баланс между вычислительной

сложностью методов и качеством получаемых интервальных оценок. Методы, основанные на

классической интервальной арифметике, зачастую могут давать чрезмерно завышенные ин-

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

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

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

интервальные оценки. С практической точки зрения интерес представляют методы, которые

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

сложности соответствующих методов.

Ранее авторами был разработан и апробирован алгоритм адаптивной интерполяции для

моделирования динамических систем с интервальными параметрами [11, с. 16-19]. Его суть

заключается в адаптивном иерархическом разбиении области неопределённости на подобласти,

для каждой из которых выполняется построение полиномиальной функции [12], интерполи-

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

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

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

976

АЛГОРИТМ АДАПТИВНОЙ ИНТЕРПОЛЯЦИИ НА РАЗРЕЖЕННЫХ СЕТКАХ

977

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

ров, определяемых положением узла в пространстве. Количество узлов в одной сетке оцени-

вается как (p + 1)^m, где p - степень интерполяционного полинома по каждому измерению,

а m - количество интервальных параметров. При увеличении m время работы алгоритма

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

применения.

В работах [13, 14] дано теоретическое обоснование алгоритма и выполнено его тестирова-

ние как на модельных задачах, так и на прикладных задачах химической кинетики. В [15]

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

ру графических процессоров с использованием технологии CUDA [16]. В [17] продемонстри-

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

анализа получающегося адаптивного разбиения области неопределённости в процессе работы

алгоритма. В работе [18] выполнен анализ вычислительных затрат алгоритма и получены ре-

комендации по выбору оптимальной степени интерполяционного полинома p. В [19] показана

возможность использования разреженных сеток в алгоритме адаптивной интерполяции для

моделирования динамических систем, содержащих большое количество интервальных пара-

метров, и продемонстрирована эффективность данного подхода. Настоящая работа посвящена

теоретическим аспектам применения разреженных сеток для численного интегрирования сис-

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

параметрами.

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

представления, интегрирования и интерполяции многомерных функций [20]. За счёт незначи-

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

зить “проклятие размерности” (экспоненциальный рост числа узлов при увеличении размерно-

сти). Классические разреженные сетки являются результатом оптимизации вычислительных

затрат для аппроксимации функций с ограниченными смешанными производными [21, с. 3]. В

[22] даётся введение в разреженные сетки и технику работы с ними, а также приводится про-

граммный код на языке программирования Python с их реализацией. Работа [23] посвящена

алгоритмам распараллеливания, а работа [24] - различным приложениям. В [25] рассматри-

ваются вопросы повышения эффективности разреженных сеток.

Существует адаптивный вариант разреженных сеток, который позволяет динамически учи-

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

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

переходами - её разрежение. Кроме того, сетки учитывают разное влияние подмножеств пе-

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

комбинацией функций от определённых подгрупп переменных, то построение сетки будет про-

исходить только на подмножествах, соответствующих этим подгруппам.

Решение системы ОДУ с интервальными начальными условиями и параметрами в каждый

момент времени можно рассматривать как некоторую многомерную вектор-функцию, которая

в общем случае является “чёрным ящиком”. В начальный момент времени эта функция три-

виальная, её компоненты являются либо константами, либо линейными функциями. Основная

цель - построить для каждого момента времени соответствующий интерполянт. От шага к

шагу интерполяционная сетка может изменяться, так как происходит её адаптация. Для то-

го чтобы вычислить значение исходной вектор-функции в произвольной точке на некотором

шаге, нет необходимости интегрировать систему ОДУ из начального состояния; вместо этого

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

сится на текущий шаг. При этом важно понимать как в результате поведёт себя глобальная

погрешность. Этот вопрос рассматривается в данной работе.

1. Постановка задачи. В общем виде рассматривается система ОДУ, состоящая из n_e

уравнений с n_p параметрами:

du_i(t)

= g_i(u₁(t),u₂(t),... ,une(t),v₁,v₂,... ,vnp,t), i = 1,n_e, t ∈ [t₀,t_N],

u_i(t₀) ∈ [u0i,u0i], i = 1,m_e,

ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 57

№7

2021

978

МОРОЗОВ, РЕВИЗНИКОВ

u_i(t₀) = u0i, i = m_e + 1,n_e,

v_i ∈ [v_i,v_i], i = 1,m_p,

(1)

где m_e - количество интервальных начальных условий, m_p - количество интервальных пара-

метров. Вектор-функция g = (g₁, g₂, . . . , gne )^т удовлетворяет всем условиям, обеспечивающим

существование и единственность решения при всех u_i(t₀) ∈ [u0i, u0i], i = 1, m_e, и v_i ∈ [v_i, v_i],

i = 1,m_p.

Сгруппируем вместе интервальные неопределённости и выполним преобразование систе-

мы (1) к автономной системе ОДУ без параметров путём добавления в систему фиктивных

уравнений:

dy_i(t)

= f_i(y₁(t),y₂(t),... ,y_n(t)), i = 1,n,

y_i(t₀) ∈ [y0i,y0i], i = 1,m,

y_i(t₀) = y0i, i = m + 1,n, t ∈ [t₀,t_N ],

(2)

где n = n_e + n_p + 1 - количество уравнений, m = m_e + m_p - количество интервальных

начальных условий,

y_i(t) = u_i(t), y0i = u0i, y0i = u0i; f_i = g_i, i = 1,m_e,

yme+i(t) = vi, yme+i = vⁱ, yme+i=vi;^fme+i≡0,i=1,mp^,

ym+i(t) = ume+i(t), ym+i = ume+i; fm+i = gme+i, i = 1,ne - me,

yne+mp+i(t) = vmp+i, yⁿe+mp+i=vmp+i;fne+mp+i≡0,i=1,np^-mp^,

y_n(t) = t, y0n = t₀; f_n ≡ 1.

Предположим, что решение системы (2) существует при всех указанных начальных дан-

ных на всём рассматриваемом отрезке и является функцией от независимой переменной t и

точечных значений интервальных начальных условий:

y(y⁰¹, y⁰², . . . , y0m, t), y0i ∈ [y0i, y0i], i = 1, m.

Рассматривается N + 1 точка из интервала интегрирования [t₀, t_N ]:

t₀ < t₁ < t₂ < ... < t_k < ... < t_N .

Требуется для каждого значения t_k, k = 0, N , построить вектор-функцию P^k(y⁰¹, y⁰², . . . , y0m)=

= (Pk1 (y⁰¹, . . . , y0m), . . . P^kn(y⁰¹, . . . , y0m)), где y0i ∈ [y0i , y0i ], i = 1, m, интерполирующую значение

y^k(y⁰¹,y⁰²,... ,y0m) = y(y⁰¹,y⁰²,... ,y0m,t_k).

2. Алгоритм адаптивной интерполяции на разреженных сетках. Выполним после-

довательное построение решения системы (2). В начальный момент времени t₀ компоненты

P0j(y⁰¹,... ,y0m), j = 1,n, вектор-функции P⁰(y⁰¹,y⁰²,... ,y0m) будут являться линейными функ-

циями

P0i(y⁰¹,y⁰²,... ,y0m) = y0i, i = 1,m,

и константами

P0i(y⁰¹,y⁰²,... ,y0m) = y0i, i = m + 1,n.

Предположим, что на k-м шаге известен интерполяционный полином P^k(y⁰¹, y⁰², . . . , y0m).

В процессе построения Pk+1 необходимо многократно вычислять исходную вектор-функцию

ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 57

№7

2021

АЛГОРИТМ АДАПТИВНОЙ ИНТЕРПОЛЯЦИИ НА РАЗРЕЖЕННЫХ СЕТКАХ

979

yk+1 в определённых точках (ŷ01, ŷ02, . . . , ŷ0m). Для этого с помощью какого-либо численного

метода интегрирования решается следующая система ОДУ:

dŷ_i(t)

= f_i(ŷ₁(t), ŷ₂(t),... , ŷ_n(t)),

ŷi(tk) =

(ŷ⁰¹, ŷ⁰²,... , ŷ0m), i = 1,n, t ∈ [t_k,tk+1],

(3)

и в качестве yk+1 берётся ŷk+1.

На каждом шаге возникают две локальные погрешности: первая - от численного метода

интегрирования системы ОДУ, а вторая - от интерполяции при вычислении P^k(ŷ⁰¹, ŷ⁰², . . . , ŷ0m).

Будем предполагать, что первая несущественно мала по сравнению со второй.

Согласно постановке задачи интерполяции [12] интерполяционный полином обязательно

должен проходить через все узловые точки. Следовательно, погрешность интерполяции в точ-

ках, соответствующих узлам, тождественно равна нулю. Поэтому при построении интерпо-

ляционной сетки на (k + 1)-м шаге погрешность от интерполяции будет возникать только в

узлах, которых не было в сетке на k-м шаге.

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

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

интерполяции.

Далее приводится описание классических разреженных сеток согласно работам [21, 22].

В основе базиса лежит функция шляпы

{

1 - |x|, если x ∈ [-1,1],

ϕ(x) =

(4)

в противном случае.

Сначала рассматривается построение интерполяционного полинома для гладкой функции

одной переменной f(x), x ∈ [0, 1], f(0) = f(1) = 0.

Сетка уровня l определяется следующим образом:

Gl = {xl,i : xl,i = i · hl,

1 ≤ i ≤ 2^l - 1, h_l = 2^-l},

и её i-му узлу ставится в соответствие базисная функция ϕ_l,i, определяемая через функцию

шляпы (4) равенством

⁽x-i·h_l)

ϕl,i(x) = ϕ

= ϕ(2^l · x - i).

h_l

Интерполяционный полином записывается в виде

∑

P (x) =

a_l,iϕ_l,i(x),

(5)

i=1

где a_l,i = f(x_l,i). Здесь мы имеем классическую кусочно-линейную интерполяцию.

Выполняется переход к иерархическому базису. Для этого базисные функции c чётными

индексами уровня k выражаются через базисные функции уровня k - 1:

ϕk,2i(x) = ϕk-1,i(x) - (ϕk,2i-1(x) + ϕk,2i+1(x))/2,

1 ≤ i ≤ 2k-1 - 1.

(6)

Полученное выражение (6) подставляется в функцию (5):

∑

P (x) =

a_k,iϕ_k,i(x),

(7)

k=1

i=1

- нечётное

где a_k,i = f(x_k,i) - (f(xk,i-1) + f(xk,i+1))/2.

ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 57

№7

2021

8^∗

980

МОРОЗОВ, РЕВИЗНИКОВ

Коэффициенты a_k,i характеризуют вторую производную функции f(x) и при возрастании

k к бесконечности стремятся к нулю как O(2-2k). В отличие от (5) слагаемые, входящие в (7),

не равнозначны: чем больше k , тем меньше вклад соответствующих слагаемых в результат.

D-мерный базис конструируется при помощи произведения одномерных базисов:

∏

ϕl,i(x) =

ϕl_j ,i_j (xj ),

1 ≤ i_j ≤ 2lj - 1, i_j - нечётное,

j=1

где l = (l₁, l₂, . . . , l_d), i = (i₁, i₂, . . . , i_d) - мультиуровень и мультииндекс базисной функции,

x = (x₁,x₂,...,x_d) - вектор переменных.

В зависимости от ограничений на мультиуровень l базисных функций, использующихся

в интерполяционном полиноме, получается или полная сетка, или разреженная сетка. Для

полной сетки уровня n_g имеем ограничение ∥l∥_∞ ≤ n_g, где ∥l∥_∞ = max l_j, а для разреженной

j=1,d

сетки уровня n_g - ограничение

∥l∥₁ ≤ n_g + d - 1,

∑_d

где ∥l∥₁ =

l_j. Сравним асимптотические характеристики сеток. Обозначим p = 2ng - 1.

j=1

Для полной сетки количество узлов оценивается как O(p^d), а погрешность - O(h2ng ); для

разреженной сетки количество узлов - O(p(log₂ p)d-1), погрешность - O(h2ng (log₂ p)d-1) [21].

Интерполяционный полином:

∑

P (x) = al,iϕl,i(x),

∥l∥₁ ≤ n_g + d - 1,

l,i

1 ≤ i_j ≤ 2lj -1, i_j - нечётное, j = 1,d, (8)

∑

∑_d

al,i =

(-1/2)

j=1

|δ_j| ×

δ₁,...,δ_d

× f(xl1,i₁+δ₁, xl₂,i₂+δ₂, . . . , xl_d,i_d+δ_d),

-1 ≤ δ_j ≤ 1, j = 1,d,

(9)

где f(x₁, x₂, . . . , x_d) - интерполируемая глад-

кая функция, определённая на единичном d-

мерном кубе и равная нулю на его границе.

Для учёта граничных значений в базис до-

бавляются базисные функции нулевого уров-

ня. Если l_j = 0, то в (8) i_j = 0, 1, а в (9) δ_j =

= 0. Учёт граничных значений эквивалентен

построению разреженных сеток на всех гранях

области определения функции меньшей раз-

мерности. В частности, для трёхмерного куба

имеем 6 граней размерности 2, 12 граней (рё-

бер) размерности 1 и 8 граней (вершин) раз-

мерности 0. В общем случае для d-мерного

куба количество граней меньшей размерности

равняется 3^d - 1.

На рис. 1 снизу продемонстрирован одно-

мерный иерархический базис, а сверху показан

пример интерполяции функции

f (x) = 0.5 + 5x - 4x²

Рис. 1. Иерархический базис.

на интервале [0, 1].

ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 57

№7

2021

АЛГОРИТМ АДАПТИВНОЙ ИНТЕРПОЛЯЦИИ НА РАЗРЕЖЕННЫХ СЕТКАХ

981

На рис. 2, а показаны узлы, соответствующие двумерным базисным функциям разного

уровня, а на рис. 2, б - все узлы разреженной сетки одновременно.

Рис. 2. Разреженная сетка третьего уровня: (а) узлы, соответствующие базисным функциям одного уровня

по отдельности; (б) все узлы сетки одновременно.

Классические разреженные сетки являются результатом оптимизации вычислительных за-

трат для аппроксимации функций с ограниченными смешанными производными. Определим

соответствующее пространство Соболева [21, 22]:

{

}

∂∥α∥1 f

H_mix(Ω) = f : Ω → R |

∥α∥_∞ ≤ 2, f|∂Ω = 0

...∂xαd ∈L2(Ω),

∂xα11 ∂x22

Согласно работе [21, с. 3], если f ∈ H_mix(Ω), то

∥f(x) - P (x)∥₂ = O(h2n

ngd-1).

Для построения решения исходной задачи область неопределённости с помощью линей-

ных преобразований трансформируется в единичный m-мерный куб. Соотношения (8) и (9)

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

∑

Pk(y⁰¹, y⁰², . . . , y0m) =

akl,iϕl,i(y⁰¹,y⁰²,... ,y0m),

l,i

|l|₁ ≤ n_g + m - 1,

1 ≤ i_j ≤ 2lj - 1, i_j - нечётное, j = 1,m,

(

)∑m

∑

j=1

|δ_j|

akl,i =

ŷ^k(y01,l

,y02,l

,...,y⁰

1,i1+δ1

2,i2+δ2

m,lm,im+δm

δ₁,...,δm

-1 ≤ δ_j ≤ 1, j = 1,m.

(10)

Предположим, что на k-м шаге решения системы (2) фиксируется локальная погреш-

ность eps:

max

|P^kj(ŷ⁰¹, ŷ⁰², . . . , ŷ0m) - ŷ^kj(ŷ⁰¹, ŷ⁰², . . . , ŷ0m)| ≤ eps,

(11)

ŷ0i∈[y0i,y0i]

i=1,m

j=1,n

ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 57

№7

2021

982

МОРОЗОВ, РЕВИЗНИКОВ

где

ŷ^k = (ŷk1,... , ŷ^kn)^т - решение соответствующей системы (3) c интервалом интегрирова-

ния t ∈ [tk-1, t_k]. Другими словами, уровень разреженной сетки n^kg выбирается минимально

возможным, так, чтобы выполнялось условие (11).

Утверждение. Пусть компоненты y^kj(y⁰¹, y⁰², . . . , y0m), j = 1, n, k = 0, N , решения

y^k(y⁰¹,... ,y0m) системы (2) принадлежат пространству H_mix([y⁰¹,y⁰¹]×[y⁰²,y⁰²]×...×[y0m,y0m]).

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

max

|P^Nj (y⁰¹, y⁰², . . . , y0m) - y^Nj (y⁰¹, y⁰², . . . , y0m)| = O(n^Ng eps),

(12)

y0i∈[y0i,y0i]

i=1,m

j=1,n

где n^Ng - уровень разреженной сетки в конечный момент времени.

Доказательство. Локальная погрешность от интерполяции возникает только на тех ша-

гах, на которых происходило уплотнение сетки. При этом, если уровень сетки на некотором

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

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

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

ню сетки в конечный момент времени n^Ng .

Одним из способов оценки глобальной погрешности является суммирование всех оценок

для локальных погрешностей, перенесённых на конец интервала интегрирования [26, с. 23].

Так как решение системы (2) удовлетворяет условию Липшица, то локальные погрешности

переносятся в конец интервала интегрирования с константой Липшица [14]:

∥y(y^k + δ, t_N - t_k) - y(y^k , t_N - t_k)∥ ≤ L∥δ∥ ≤ L eps,

где L - константа Липшица, δ - локальная погрешность.

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

сверху как n^Ng L eps, и в результате получается оценка глобальной погрешности (12). Утвер-

ждение доказано.

Можно обобщить утверждение на случай ненулевых граничных значений, для этого необ-

ходимо рассмотреть отдельно границы области неопределённости.

На практике используется адаптивный вариант разреженных сеток. Если |al,i| > ϵ, то в

базис добавляются базисные функции на один уровень выше:

ϕ(l₁+1,l₂,...,l_d),(2i₁-1,i₂,...,i_d),

ϕ(l₁+1,l₂,...,l_d),(2i₁+1,i₂,...,i_d),

ϕ(l₁,l₂+1,...,l_d),(i₁,2i₂-1,...,i_d),

ϕ(l₁,l₂+1,...,l_d),(i₁,2i₂+1,...,i_d),

ϕ(l₁,l₂,...,l_d+1),(i₁,i₂,...,2i_d-1),

ϕ(l₁,l₂,...,l_d+1),(i₁,i₂,...,2i_d+1).

Коэффициенты al,i характеризуют смешанные производные и с увеличением уровня убывают,

поэтому процесс адаптации сходится.

3. Результаты. Выполняется численное интегрирование нескольких автономных систем

ОДУ с интервальными начальными условиями и параметрами.

Вначале рассматривается система ОДУ с двумя интервальными начальными условиями,

которая соответствует модели Лотки-Вольтерры [27, с. 30]:

x^′ = -0.9x + 0.5xy, y^′ = 0.7y - 0.8xy, t ∈ [0,50],

x(0) ∈ [0.5, 0.7], y(0) ∈ [0.3, 5.0].

(13)

Параметр ϵ выбирается равным 10-3.

ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 57

№7

2021

АЛГОРИТМ АДАПТИВНОЙ ИНТЕРПОЛЯЦИИ НА РАЗРЕЖЕННЫХ СЕТКАХ

983

На рис. 3, а показано множество решений системы (13) в различные моменты времени, а

на рис. 3, б - узлы сетки. Точки на двух рисунках соответствуют друг другу. На рис. 3, в и

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

Рис. 3. Результаты решения системы (13): (a) - множество решений в различные моменты времени; (б) -

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

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

Здесь наблюдается корреляция между графиками, увеличение уровня приводит к увели-

чению погрешности, что соответствует ранее сформулированному утверждению. Для оценки

погрешности в начальный момент времени случайным образом генерируется тестовое мно-

жество точек из области неопределённости параметров системы и строятся соответствующие

решения, с которыми выполняется сравнение.

Далее рассмотрим автономные системы ОДУ, для которых решение может быть найдено

аналитически. Первая система [11]:

(√

)

(

√

)

x^′ = a

y^′ = a -

t ∈ [0,4],

x(0) = x₀ ∈ [-1, 1], y(0) = y₀ ∈ [-1, 1], a ∈ [-1, 1].

(14)

Общее решение системы в (14) имеет следующее аналитическое представление:

(√

)(

)

⁽ at

⁽ at⁾⁾

x(x₀, y₀, a, t) = exp

at y₀ sin

+ x0 cos

(√

)(

)

⁽ at

⁽ at⁾⁾

y(x₀, y₀, a, t) = exp

at y₀ cos

- x0 sin

(15)

Проанализируем решение (15). Оно линейно зависит от начальных условий x₀ и y₀ и

нелинейно - от параметра a. Старшие производные отличны от нуля только для параметра a,

и, соответственно, коэффициенты (10) ненулевые только при l₁ = 0 и l₂ = 0. Следовательно,

уплотнение сетки будет происходить на четырёх подмножествах:

{-1} × {-1} × a,

{-1} × {1} × a,

{1} × {-1} × a,

{1} × {1} × a.

ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 57

№7

2021

984

МОРОЗОВ, РЕВИЗНИКОВ

На рис. 4, а и 4, б показаны множество решений системы (14) и адаптивная разреженная

сетка в различные моменты времени. Уплотнение сетки происходит только на рёбрах, которые

соответствуют параметру a. На рис. 4, в и 4, г представлены зависимости уровня сетки и

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

увеличению погрешности.

Рис. 4. Результаты решения системы (14): (a) - множество решений в различные моменты времени; (б) -

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

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

Введём в систему (14) ещё один интервальный параметр и получим вторую систему:

√

x^′ = a(bx +

1 - b²y), y^′ = a(-

1 - b²x + by), t ∈ [0,3],

x(0) = x₀ ∈ [-1, 1], y(0) = y₀ ∈ [-1, 1], a ∈ [0, 1], b ∈ [0, 1].

(16)

Общее решение системы в (16) аналитически записывается следующим образом:

√

x(x₀, y₀, a, b, t) = exp (abt)(y₀ sin(a

1 - b²t) + x0 cos(a

1 - b²t)),

√

y(x₀, y₀, a, b, t) = exp (abt)(y₀ cos(a

1 - b²t) - x0 sin(a

1 - b²t)).

В отличие от предыдущего примера здесь старшие смешанные производные не равны нулю

сразу по двум параметрам: a и b. Следовательно, уплотнение сетки будет происходить только

на подмножествах, соответствующих этим параметрам:

{-1} × {-1} × a × b,

{-1} × {1} × a × b,

{1} × {-1} × a × b,

{1} × {1} × a × b.

Построение решения выполняется при ϵ = 5 · 10-3. На рис. 5, а показано множество реше-

ний системы (16) в различные моменты времени. Решение системы в каждый момент времени

принадлежит пространству R⁴, поэтому изображение на фазовой плоскости является лишь

проекцией. На рис. 5, б показана адаптивная разреженная сетка, получающаяся в процессе ин-

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

параллелепипедом, поэтому, по аналогии с рис. 5, а, здесь показана проекция.

ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 57

№7

2021

АЛГОРИТМ АДАПТИВНОЙ ИНТЕРПОЛЯЦИИ НА РАЗРЕЖЕННЫХ СЕТКАХ

985

Рис. 5. Результаты решения системы (16): (a) - множество решений в различные моменты времени; (б) - узлы

адаптивной разреженной сетки.

Выполним сравнение вычислительных затрат алгоритма адаптивной интерполяции на ос-

нове разреженных сеток с классическим вариантом алгоритма. В работе [13] определён кри-

терий I, который характеризует вычислительные затраты и численно равен числу решённых

неинтервальных систем (2) при конкретных точечных значениях интервальных параметров.

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

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

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

вариант алгоритма.

Таблица. Результаты работы классического алгоритма адаптив-

ной интерполяции при p = 3 и алгоритма адаптивной интерполя-

ции на основе разреженных сеток

Алгоритм на основе

Система

Классический алгоритм

разреженных сеток

ОДУ

error

(13)

9.89 · 10-3

1 560

7.71 · 10-3

820

(14)

1.62 · 10-2

1 318

1.2 · 10-2

521

(16)

1.3 · 10-1

14 549

9.45 · 10-2

1 968

Заключение. В работе рассмотрены теоретические аспекты обобщения алгоритма адап-

тивной интерполяции на случай большого числа интервальных параметров с использованием

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

ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 57

№7

2021

986

МОРОЗОВ, РЕВИЗНИКОВ

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

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

которых строится интерполяционный полином в форме Лагранжа. При этом количество узлов

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

вопрос замены регулярных сеток на разреженные сетки. В процессе работы алгоритма проис-

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

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

глобальной погрешности и показано, что она линейно зависит от текущего уровня сетки. Про-

ведена апробация алгоритма адаптивной интерполяции на разреженных сетках на нескольких

системах ОДУ, демонстрирующая эффективность рассматриваемого подхода. Полученные ре-

зультаты согласуются с теоретическими оценками.

Работа выполнена при финансовой поддержке Министерства науки и высшего образования

Российской Федерации (проект 075-15-2020-799).

СПИСОК ЛИТЕРАТУРЫ

1. Мозохина А.С., Мухин С.И. Некоторые точные решения задачи о течении жидкости в сокращаю-

щемся эластичном сосуде // Мат. моделирование. 2019. Т. 31. № 3. С. 124-140.

2. Moore R. Interval Analysis. New Jersey, 1966.

3. Moore R.E., Kearfott R.B., Cloud M.J. Introduction to Interval Analysis. Philadelphia, 2009.

4. Шарый С.П. Конечномерный интервальный анализ. Новосибирск, 2019.

5. Eijgenraam P. The Solution of Initial Value Problems Using Interval Arithmetic. Math. Centre Tracts.

№ 144. Stichting Math. Centrum. Amsterdam, 1981.

6. Lohner R.J. Enclosing the solutions of ordinary initial and boundary value problems // Comput.

Arithmetic: Scientific Comput. and Program. Lang. 1987. P. 255-286.

7. Черноусько Ф.Л. Оценивание фазовых состояний динамических систем. Метод эллипсоидов. М.,

1988.

8. Makino K., Berz M. Models and their applications // Numerical Software Verification 2017: conference.

Heidelberg, Germany, July 22-23, 2017. Springer International Publishing AG, 2017. P. 3-13.

9. Nataraj P.S.V., Sondur S. The extrapolated Taylor Model // Reliable Computing. July, 2011. P. 251-278.

10. Рогалев А.Н. Гарантированные методы решения систем обыкновенных дифференциальных урав-

нений на основе преобразования символьных формул // Вычислит. технологии. 2003. Т. 8. № 5.

С. 102-116.

11. Морозов А.Ю., Ревизников Д.Л. Методы компьютерного моделирования динамических систем с

интервальными параметрами. М., 2019.

12. Sauer T., Xu Y. On multivariate Lagrange interpolation // Math. of Comput. 1995. V. 64. № 211.

P. 1147-1170.

13. Морозов А.Ю., Ревизников Д.Л. Алгоритм адаптивной интерполяции на основе kd-дерева для чис-

ленного интегрирования систем обыкновенных дифференциальных уравнений с интервальными

начальными условиями // Дифференц. уравнения. 2018. Т. 54. № 7. С. 963-974.

14. Морозов А.Ю., Ревизников Д.Л., Гидаспов В.Ю. Алгоритм адаптивной интерполяции на основе

kd-дерева для решения задач химической кинетики с интервальными параметрами // Мат. моде-

лирование. 2018. Т. 30. № 12. С. 129-144.

15. Morozov A.Yu., Reviznikov D.L. Modelling of dynamic systems with interval parameters on graphic

processors // Program. Ingeneria. 2019. V. 10. № 2. P. 69-76.

16. Khoirudin, Shun-Liang J. Gpu application in cuda memory // Adv. Comput.: An Int. J. 2015. V. 6. № 2.

P. 1-10.

17. Morozov A.Yu., Reviznikov D.L. Modeling of dynamic systems with interval parameters in the presence

of singularities // Rus. J. Nonlin. Dyn. 2020. V. 16. № 3. P. 479-490.

18. Морозов А.Ю., Журавлев А.А., Ревизников Д.Л. Анализ и оптимизация алгоритма адаптивной

интерполяции численного решения систем обыкновенных дифференциальных уравнений с интер-

вальными параметрами // Дифференц. уравнения. 2020. Т. 56. № 7. С. 960-974.

19. Morozov A.Yu., Zhuravlev A.A., Reviznikov D.L. Sparse grid adaptive interpolation in problems of

modeling dynamic systems with interval parameters // Mathematics. 2021. V. 9. № 4. P. 298.

20. Смоляк С.А. Квадратурные и интерполяционные формулы на тензорных произведениях некоторых

классов функций // Докл. АН СССР. 1963. Т. 148. № 5. С. 1042-1045.

ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 57

№7

2021

АЛГОРИТМ АДАПТИВНОЙ ИНТЕРПОЛЯЦИИ НА РАЗРЕЖЕННЫХ СЕТКАХ

987

21. Gerstner T., Griebel M. Sparse grids // Encyclopedia of Quantitative Finance / Ed. R. Cont. New York,

2010.

22. Garcke J. Sparse grids in a nutshell. Sparse grids and applications // Lect. Not. in Comput. Sci. and

Engin. V. 88. Berlin; Heidelberg, 2013. P. 57-80.

23. Brumm J., Scheidegger S. Using adaptive sparse grids to solve high-dimensional dynamic models

// Econometrica. 2017. V. 85. № 5. P. 1575-1612.

24. Bungatrz H-J., Griebel M. Sparse grids // Acta Numerica. 2004. V. 13. № 1. P. 147-269.

25. Judd K.L., Maliar L., Maliar S., Valero R. Smolyak method for solving dynamic economic models:

Lagrange interpolation, anisotropic grid and adaptive domain // J. of Economic Dyn. and Cont. 2014.

V. 44. P. 92-123.

26. Niesen J., Hall T. On the Global Error of Discretization Methods for Ordinary Differential Equations.

Cambridge, 2004.

27. Арнольд В.И. Обыкновенные дифференциальные уравнения. М., 2012.

Федеральный исследовательский центр

Поступила в редакцию 15.02.2021 г.

“Информатика и управление”

После доработки 15.02.2021 г.

Российской академии наук, г. Москва,

Принята к публикации 27.04.2021 г.

Московский авиационный институт

(национальный исследовательский университет)

ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 57

№7

2021