Вариационный подход к решению обратной задачи для интегральной динамической модели

Исследование выполнено за счет гранта Российского научного фонда № 2521-00743, ject/25-21-00743/

Рассмотрим интегральную динамическую модель:

/ t K (t

,s,x ⁽ s))ds = g^t),

0 <  s <  t <  T, g(0) = 0,

с ядром K(t, s, x(s)), имеющим разрывы I рода на множестве гладких кривых:

K i (t, s)G i (s, x(s)),   (t, s) G m i ,

K(t, s') = <

K n (t, s]G n (s, x(s)),   (t, s) G m _n .

Здесь m i = { t,s | a _{i- 1} (t) < s <  a i (t) } , a o (t) = 0, a _n (t) = t, i = 1,n, a i (t), g(t) G C 0 T ] . Компоненты K i (t, s) ядра K непрерывно дифференцируемы по переменной t при (t,s) G cl(m i ), где cl(m i ) — замыкание множества m i , а K _n (t,t) = 0. Функции G i (s, x(s)), i = 1,n удовлетворяют условию Липшица по второй переменной. Линии разрыва ядра удовлетворяют следующим условиям:

a i (0) = 0, 0 < a 1 (t) < a ₂ (t) <  • • • < a _n-1 (t) < t, t G (0, T ].

Пусть также функции a i (t),..., a _n-i (t) монотонно возрастают и

0 <  a 0 (0) < ••• <  а П- 1 (0) <  1.

Классическая постановка задачи для модели (1) - (2) заключается в определении x(t) при известных остальных компонентах. В этом случае мы имеем дело с уравнением Вольтерра I рода с кусочно-непрерывным ядром специального вида. Такие уравнения называются слаборегулярными, их теоретическому и численному исследованию посвящен ряд современных работ ([1–3], а также приведенные там ссылки).

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

Модели, описываемые интегральными уравнениями с неизвестными пределами интегрирования, берут свое начало в работах В. М. Глушкова [4] , их экономические приложения исследуются в работах N. Hritonenko и Yu. Yatsenko ( [5] , [6] ), ряд прямых и итерационных численных методов предложен в работах [7] – [12] .

Настоящая работа является продолжением работы [7] и посвящена построению вариационного метода решения обратной задачи для модели (1) - (2) в скалярном случае.

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

Рассмотрим задачу (1) - (2) при n = 2, состоящую в определении неизвестной функции разрыва a i (t) = a(t). Имеем функциональноинтегральное уравнение следующего вида:

^a ( t )

/ h l (t,s) о

t ds + У h2(t,s) ds = f (t), a(t)

0 6 t 6 T.

Зная функции h i (t, s), h ₂ (t, s') и f (t), поставим при известном значении параметра T задачу о воостановлении неизвестной функции a(t) в следующих дополнительных предположениях.

• 1.    Функция a(t) при 0 6 t 6 T является неубывающей. Если предположить, что неизвестная функция дифференцируема во всех точках t Е (0, T ), то тем самым должно выполняться условие:

  a ⁰ (t) >  0 при t Е (0,T ).

  
  
  
  

• 2.    При всех 0 6 t 6 T функция a(t) удовлетворяет двустроннему неравенству

0 6 a(t) 6 t.                             (5)

Здесь и далее без ограничения общности будем считать, что T = 1. В самом деле, в случае, если T = 1, можно выполнить в уравнении (3) замену независимой переменной t :

t0 = t/T, обеспечив, таким образом, справедливость неравенства 0 6 t0 6 1.

• 2    Описание метода

Предположим, что подынтегральные функции h^tt, s) и h 2 (t. s') уравнения (3) таковы, что найдутся функции H i (t. s) и H 2 (t. s), удовлетворяющие условиям:

dH i             дН2 .    ,

~Q^ = ^{h l (t,s)} ,         = h 2 ^(t,s) .

Тогда в результате применения к интегралам в левой части уравнения (3) формулы Ньютона — Лейбница можно переписать его в следующем эквивалентном виде:

H i (t.s)

s = a ( t )

+ ^H 2 ^(t . ^{s )} s =0

s=t s=a(t)

= f (t).

Заметим, что если для заданных функций h i (t,s) и h 2 (t,s) не находится функций H i (t,s) и H 2 (t,s), удовлетворяющих условиям (6) , то описываемый метод всё же может быть применён, если приближённо заменить в уравнении (3) функции h i (t, s) , h 2 (t, s) аппроксимирующими их функциями h _i (t. s), h ₂ (t. s), допускающими интегрирование по переменной s ; для построения таких аппроксимаций может быть, в частности, эффективно использован аппарат степенных рядов или рядов Фурье по подходящим ортогональным системам функций. В этом случае после нахождения в результате интегрирования таких функций H i (t.s), H ₂ (t,s), что H i = d h _i /ds, H ₂ = dh ₂ /ds может быть аналогичным образом выполнен переход к уравнению (7) .

Уравнение (7) в развёрнутом виде переписывается следующим образом:

H i (t. a(t)) - H i (t. 0) + H 2 (t. t) - H 2 (t. a(t)) = f (t).

Наконец, обозначив f (t) = f (t) + Hi(t. 0) - H2(t.t).

получим следующее уравнение, нелинейное относительно искомой функции a(t).

H i (t.a(t)) - H 2 (t.a(t)) = f (t).                     (8)

Уравнение (8) лежит в основе предлагаемого численного метода восстановления функции a(t).

Искомую функцию a(t) будем аппроксимировать функцией a(t), определяемой следующим образом:

m

a(t) = X C j t ^j .

j =i

Здесь

• •    m — параметр метода, натуральное число;

• •    C 1 , . . . , C _m — вещественнозначные коэффициенты, подлежащие определению.

Необходимо заметить, что, потребовав, чтобы аппроксимация a(t) удовлетворяла при 0 6 t 6 1 тем же дополнительным условиям неубывания функции a(t) и справедливости неравенства a(t) 6 t, заметим, что при t = 0 должно выполняться условие а(0) = 0. По этой причине в выражении (9) для функции a(t) отсуствует свободный член С о , который в силу а(0) = 0 должен быть равен нулю.

Используя аппроксимацию (9) , заменим уравнение (8) следующим приближенным уравнением:

H i (t,a(t)) - H (t,a(t)) = f(t).

Подставив в это уравнение представление (9) , получим:

m

m

H i  t, £ C j t ^j   - H t,£ C j t ^j  = f(t)

Введём в рассмотрение функцию Ф (t; C i ,..., C m ) в соответствии с формулой:

m

m

Ф(t; Ci,...,Cm)= Hi (t^Cjtj j=1

H 2  tJ^C j t j  - f (t)

j =1

Коэффициенты C j (j = 1, m) будем выбирать таким образом, чтобы стремилась к минимуму функция е (C i ,..., C m ), которую определим следующим образом.

Введём на сегменте t G [0,1] равномерную сетку из узлов t i = ih, где i = 1,N , h = 1/N и N — целое положительное число, являющееся параметром метода.

Коэффициенты C 1 , . . . , C m будем фиксировать таким образом, чтобы своего минимума достигала функция:

N е (Ci,..., Cm) = ]□ Ф(^; Ci,...,Cm)

N

H 1

m

m

t i CP^ t j   - H 2 tP^Ct j   - f(t i )

Таким образом, задача приближённого восстановления функции a(t) сведена к задаче минимизации нелинейной функции е (C 1 ,..., C m ). Эта задача может быть решена при помощи любого современного численного метода нелинейной оптимизации. При этом если функция е (C 1 ,..., C m ) имеет несколько локальных минимумов, то среди них необходимо выбрать наименьший.

Для корректного решения задачи минимизации функции в общем случае необходимо учесть наличие дополнительных условий 1, 2, наложенных на искомую функцию a(t).

• 1.    Условие a'(t) >  0 приводит к неравенству:

• 2.    Условие 0 6 a(t) 6 t порождает неравенство:

m

X jCt j ^-1 >  0.

j =1

Фиксируя в этом неравенстве последовательно t = t i , где i = 0,... ,N , получаем следующий набор дополнительных условий:

C 1 >  0, X jC j t j^{- 1} >  0, i = ТЖ          (13)

j =1

m

X C j t j 6 t.

j =1

Последовательно полагая в указанном неравенстве t = t i , где i = T,... ,N , приходим к следующей последовательности дополнительных условий:

m

X C j t j 6 t i .                       (14)

j =1

Тем самым приходим к задаче минимизации определенной формулой (12) функции е (C 1 ,..., C m ) при условиях (13) и (14) . Для завершения алгоритма найденные в результате решения задачи минимизации функции е ( C 1 , . . . , C _m ) значения коэффициентов C 1 , . . . , C _m необходимо подставить в формулу (9) для получения аппроксимации a(t) функции a(t).

Необходимо отметить, что предлагаемый метод в описанном виде имеет существенное ограничение применимости: хотя теоретически имеется возможность получать решение задачи с произвольной точностью за счёт увеличения степени m аппроксимирующего функцию a(t) полинома, на практике увеличение m сверх какого-то сравнительно небольшого значения приводит не к увеличению, а напротив, к падению точности и эффективности метода. Тем не менее имеется возможность достаточно простого обобщения описанного метода, позволяющего выполнять восстановление функции a(t) с произвольно заданной точностью; далее опишем это обобщение.

Сначала введём на сегменте t G [0,1] равномерную сетку из узлов t i = ih i , где i = 0, N i , N i — фиксированное целое положительное число и h i = 1/N i . Затем каждый из N i полученных сегментов [t i ,t i +i ] равной длины h i , где i = 0, N i — 1, в свою очередь, разбивается на N 2 промежутков [t i,j ,t i,j +i ]; здесь t i,j = t i + jh 2 , где j = 0,N 2 — 1, N — фиксированное целое положительное число и h 2 = h i /N 2 .

Вновь, как и ранее, зафиксировав не слишком большое целое положительное число m, являющееся параметром метода, будем аппроксимировать функцию a(t) с помощью функции a(t), определяемой следующим образом:

m

a(t) = ^C i,k t^k,если t i 6 t i +i ,               (15)

k =0

где i = 0, N i — 1, а коэффициенты C i^ , где i = 0, N i — 1 и k = 0,m, подлежат определению; заметим, что в силу условия а(0) = 0 можно сразу же зафиксировать С о , о = 0.

Подставив определенную таким образом функцию a(t) в уравнение (10) , введём в рассмотрение функцию:

Ф (t; C ) = | H (t, a(t)) — H 2 (t, a(t)) — f(t) i 2 .           (16)

Здесь C обозначает набор коэффициентов

C 0 , 1 , . . . C 0 ,m , C 1 , 0 , C 1 , 1 , . . . , C 1 ,m , . . . , C N 1 - 1 , 0 , C N 1 - 1 , 1 , . . . , C N 1 - 1 ,m .

Запишем определённую формулой (16) функцию Ф(t, C ) в точках t = t i,j , где i = 0, N i — 1 и j = 1,N 2 , и будем искать такой набор значений C , который доставлял бы глобальный минимум функции:

N 1 - 1 N 2

^г( С )= XX Ф(^; C ).             (17)

i =0 j =1

Минимизацию функции е( С ) будем проводить при следующих дополнительных условиях-ограничениях.

• 1.    Требование неубывания функции a(t) при 0 6 t 6 1 приводит к неравенствам:

• 2.    Требование a(t) 6 t при 0 6 t 6 1 означает справедливость неравенств:

• 3.    Наконец, потребовав непрерывности функции a(t), получим следующие дополнительные условия:

m

X kC i,k t k- 1 >  0 при i = 0,N i - 1, j = Ш.      (18)

k =1

m

X C i,k t kj 6 t при i = 0,N i - 1, j = 1,N₂.       (19)

k =0

mm

X C ik tw = X C i +i ,k t k +1 , 0 , i = 0,N 1 - 2.      (20)

k =0             k =0

Искомый набор коэффициентов C , полностью определяющий аппроксимацию a(t) искомой функции a(t), может быть получен в результате минимизации функции е ( С ) при условиях (18) - (20) . Подстановка найденного набора значений C в формулу (15) завершает решение задачи.

3 Численная иллюстрация

Для иллюстрации эффективности предлагаемого алгоритма выполним решение ряда модельных примеров. Первый модельный пример определяется следующими исходными данными:

h 1 = cos(t + s), h ₂ = sin(t — s), f (t) = 1 — sin(t) + sin (t + sin ³ (t) ) — cos (t — sin ³ (t) ) .

Точное решение задачи даётся функцией:

a(t) = sin ³ (t).

В соответствии с описанием метода зафиксируем T = 1, N = 100. Для решения задачи ограниченной нелинейной оптимизации использован метод SQP последовательного квадратичного программирования. Расчеты произведены в системе компьютерной математики Maple с количеством значащих цифр в вычислениях равным 16 (задан параметр Digits:=16 ).

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

m	ε min	δ
1	0.140	0.305
2	1.08 · 10 - 2	4.62 · 10 - 2
3	8.16 · 10 - 4	2.46 · 10 - 2
4	3.87 · 10 - 7	4.98 · 10 - 4
5	1.30 · 10 - 7	2.60 · 10 - 4
6	1.68 · 10 - 10	1.07 · 10 - 5
7	2.08 · 10 - 11	6.91 · 10 - 6

В этой таблице:

• •    m обозначает степень полинома в формуле (9) , используемого для аппроксимации искомой функции α(t);

• •    ε _mi_n обозначает минимальное значение определённой формулой функции (12) , соответствующее набору значений коэффициентов C 1 , . . . , C _m , который определяет приближенное решение α˜(t) поставленной задачи восстановления функции α(t);

• •    δ обозначает погрешность приближенного решения задачи, вычисляемую по формуле

δ = min k α(jτ) - α˜(jτ) k , τ = 1/M.

j =1 ,...,M

При численных расчётах значение M было зафиксировано равным M = 10 ³ .

В качестве примера на рисунке 1 приведён график точного и приближённого решения задачи при m = 3. При этом сплошной линией на графике показано точное решение задачи α(t), а пунктирной линией показано приближённое решение задачи α˜(t).

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

h 1 = (t + s) ² , h 2 = _{( t-s} ¹ ₎ 2 ₊₁ , f (t) = t 3 ⁽e ^t-1 + ^{1) 3} — t 3 — arctg (( e t^{- 1} — 1 ) t) .

Точное решение задачи 2 определяется функцией:

α(t) = te ^{t- 1} .

Рис. 1. Решение первого модельного примера при m = 3

При расчётах согласно описанию метода было зафиксировано T = 1. Значение N, как и ранее, было принято равным 100: N = 100.

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

m	ε min	δ
1	0.580	0.155
2	1.08 · 10 - 2	1.48 · 10 - 2
3	4.15 · 10 - 5	2.04 · 10 - 3
4	1.75 · 10 - 7	1.03 · 10 - 4
5	7.94 · 10 - 9	1.74 · 10 - 5
6	4.30 · 10 - 12	9.66 · 10 - 7

В качестве примера на рисунке 2 показан график точного и приближённого решения задачи при m = 2.

Рис. 2. Решение второго модельного примера при m = 2

Заключение

По представленным выше результатам расчетов можно судить об эффективности предложенного вариационного подхода к решению обратной задачи (3) . Даже при небольших степенях аппроксимирующих полиномов можно получить результат с приемлемой точностью. Дальнейшее увеличение степени m является нецелесообразным. Для более точных расчетов применяется аппроксимация локальными сплайнами с небольшими степенями компонентов. Нужно отметить, что предложенный подход дает хорошую альтернативу прямым методам дискретизации подобных нелинейных функционально-интегральных уравнений. Дело в том, что ввиду наличия неизвестных в пределах интегральных операторов возникает необходимость аппроксимации интегралов с неизвестными областями интегрирования. И тут приходится либо применять самые простейшие квадратуры (формулу прямоугольников), либо усугублять проблему ветвления решений. Дальнейшее развитие авторы видят в обобщении предложенного подхода на системы уравнений с произвольным количеством кривых разрыва, введенные в работе [7] .