Градиентная идентификация эволюционных сеточных задач Текст научной статьи по специальности «Математика»

Труды Петрозаводского государственного университета

Серия “Математика” Выпуск 18, 2011

УДК 517.954

С. В. Маничева, И. А. Чернов

ГРАДИЕНТНАЯ ИДЕНТИФИКАЦИЯ ЭВОЛЮЦИОННЫХ СЕТОЧНЫХ ЗАДАЧ

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

§ 1. Введение

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

© С. В. Маничева, И. А. Чернов, 2011

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

Введем в прямоугольнике [0,1] х [0,1] равномерную сетку Б с шагами т и Н. Ее узлы (п, г) — это точки (Ьп, х*) = (пт, гк) при 0 < п < N,

0 < г < I, Nт =1, 1к =1. Пусть Б — подмножество внутренних узлов (0 < п < N, 0 < г < I). Значение сеточной функции и в узле (п, г) обозначим игп. Совокупность значений игп при фиксированном п обозначим ип, а объединение ип при п < т — через ип<т. Рассмотрим сеточную задачу (систему алгебраических уравнений)

/п(ит<п,|) = 0, (п, г) е Б, (1)

где функции /т гладкие, а вектор £ е Дм — это вектор параметров. Предполагаем, что эта система разрешима при любом £ е Дм.

Подчеркнем «независимость от будущего» — независимость функций /п от икт при т > п, которую предполагаем; при трактовке индекса п как времени указанное свойство и означает эволюционность задачи. В самом деле, при п = 0 функции /О(й0, £) при заданном £ зависят только от иО, поэтому, решая систему /0(ио, £) = 0, определяем иО при 0 < г < I. Далее решаем систему /|(мт<1,£) = 0 относительно и[, с учетом известных £ и и0. Продолжая далее, определяем последовательно все игп. Итак, эволюционная задача размерности ^ +1) х (I +1) сводится к N + 1 задаче размерности I +1 без каких-либо дополнительных предположений.

К таким задачам сводятся, в частности, разностные схемы для параболических краевых задач [4]. В этом случае функции / являются сеточными аппроксимациями уравнения в частных производных и краевых условий. Зависимость / от п и г связана с зависимостью коэффициентов исходной задачи от времени и точки пространства, а зависимость от предыстории — со внутренними состояниями различной физической природы [5]. Отметим, что перенумерация узлов позволяет записать в виде (1) также и разностные схемы с несколькими пространственными переменными [6].

Система (1) позволяет определить и1п для всех (п, г) е Б N при заданном £; для вычисления £ требуется дополнительная информация. Предположим, что вектор £ должен минимизировать функционал

(2)

Функцию Г предполагаем гладкой и ограниченной снизу на пространстве RN+1+2+м. Если речь идет об измерениях, то фактически функция Г зависит не от всех значений игп, а только от некоторых, например, только от иП — значений на поверхности.

Имеем обычную задачу условной оптимизации большой размерности с системой ограничений особой структуры. Поскольку аналитическое решение системы (1) невозможно или очень громоздко, то применим метод множителей Лагранжа, перейдя к задаче безусловной оптимизации функции

Множитель Л либо равен нулю, либо может быть положен равным единице. Случай Л = 0 исключается требованием линейной независимости градиентов Уи/п — векторов частных производных по перемен-

полагаем Л = 1.

Частные производные функции Лагранжа Ь можно разделить на три группы: по переменным игп, по множителям Лагранжа Хгп и по параметрам — компонентам вектора £. Условие обнуления производных второй группы равносильно системе (1). Найдем и приравняем нулю производные второй группы:

Отметим, что нижний индекс суммирования равен п в силу эволю-ционности. Также слагаемое, содержащее производные от Г, равно нулю в большинстве узлов, если функционал зависит только от некоторых значений; это позволяет существенно сократить вычисления в прикладных задачах.

Наконец, производные по параметрам имеют вид

§ 3. Метод множителей Лагранжа

N I

Ь(ит^, е, А П) = ЛГ + £ ]Г лп/п.

(3)

п=0 і=0

ным игп — на решении задачи (1) при любых £ Є Rм. В дальнейшем

N I

(4)

Это — компоненты градиента функции Г в подпространстве параметров Дм в любой точке £ € Дм.

§ 4. Декомпозиция

Пусть точка £ € Дм задана. Выше уже описана декомпозиция системы уравнений (1), сводящейся к последовательному решению подсистем для фиксированного п в ходе возрастания п. Это позволяет эффективно вычислить значения переменных «п. Приравнивая нулю выражения (4), получаем систему линейных уравнений относительно Ап размерности (Ж +1) х (I + 1), также допускающую декомпозицию. Эта система невырождена в силу предположения о линейной независимости векторов Уи/. Параметры £& заданы, а значения «п уже найдены; подсистема (4) при п = N содержит только неизвестные А^ при

0 < г < N и поэтому может быть решена независимо. Но тогда подсистема (4) при п = N — 1 содержит только неизвестные А^-1 и уже вычисленные А^. Продолжая далее, последовательно находим Ап при фиксированном п в порядке его убывания. Выражения (5) позволяют явно найти градиент У Г функции Г в пространстве параметров.

Таким образом, учет структуры системы ограничений позволяет эффективно вычислить градиент, который может быть использован, например, для методов спуска или тяжелого шарика. Метод тяжелого шарика (движение в пространстве параметров вдоль траектории системы ^£ + 7^ + УГ = 0, где точкой обозначено дифференцирование по параметру траектории, а коэффициенты «трения» 7 и «инерции» ц подбираются) выглядит предпочтительнее градиентного спуска из-за «инерции», которая позволяет избежать остановки в критических точках, не являющихся минимумами, а при удачном выборе коэффициента «трения» 7 — и миновать «ложные минимумы». Разумеется, система дифференциальных уравнений, описывающая движение тяжелого шарика, решается численно.

§ 5. Частные случаи

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

Пусть сетка состоит из шести точек (п, г), п = 0,1, г = 0,1, 2.

Функции /п заданы следующим образом:

Го = «О — 1 /0 = «1 — 1 /1 = «1— (2 + £1 )и1 + «1 + Clu0,

/2 = м1 — «1 + £1£2(«1)2.

Такой выбор связан с аппроксимацией уравнения диффузии на сетке: при единичных шагах сетки /11 является разностной схемой для уравнения диффузии ((£1)-1 — это коэффициент диффузии), /О задают постоянное начальное распределение, /0 — граничное условие I рода, а /2 — нелинейное граничное условие «объемной десорбции» [1].

Решим систему уравнений (1). Очевидно, «1 = «0 = 1 для всех г = 0,1, 2, «1 является положительным корнем квадратного уравнения

(2 + £1)£1£2 («1)2 + (1 + £1 )и1 — (1 + £1) = 0

(при £& > 0, что подразумевается, единственный вещественный положительный корень существует), а и1 выражается через и2:

1 1 + £1 + «1

«1 =-----------1.

1 2 + £1

Функционал зададим формулой

Г = (£2<“2»2 — 9 Г

Это среднеквадратичное отклонение «квадратичного потока» С!^!)2 от модельного значения, соответствующего значению ад2 = 1/3, которое получается при £1 = 1, £2 =4.

Решим систему (4). Поскольку «[] и ад2 входят только в /° ив /О соответственно, то А° = А2 = 0. Величина ад1 линейно входит в / и //, так что получаем уравнение А° + £1 А1 = 0. Аналогично получаются уравнения А° + А1 =0 и (2 + £1)А1 + А2 =0. Из последнего уравнения А1 выразим через А2. Наконец, «1 нелинейно входит в формулы для функционала Г и функций // и /2, что приводит к уравнению

^и2 ^Ыи2)2 — 9^ + А1 + А1(1 + ) = °.

После подстановки A11 получаем

2_ (2 + Сі) (ЭЫ^О2 - 4)

A] = —-

9 1+ Сі + 2(2 + a^u]-

Теперь найдем производные по параметрам по формулам (5):

Ц = Л2С2(«;)2 - A]u1, Ц = 2(u;)2 (dCu?)2 - 9) + AQ(uQ)2.

Очевидно, эти выражения обращаются в ноль в точке £ = (1, 4), так как при этом u ] = 1/3. Нетрудно убедиться, что F = О на кривой

С2 = 36 (і - 4^Сі)-2 , 4^С 1< 9.

V 1 + Сі ) 1 + Сі 1

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

Теперь рассмотрим частный случай (1): двухслойную разностную схему для эволюционного уравнения, разрешенного относительно производной по времени. Требование «двуслойности» означает, что для аппроксимации производных используются значения сеточного решения только на двух слоях, то есть un и un_ і при фиксированном n и произвольных i. Итак, функции /n теперь имеют вид

/i i i i i / —* —* \

n = un - vnun- 1 - gn(un, un- 1 ),

где vn — числовой коэффициент (он может понадобиться, например, для задания краевых условий — в этих случаях он равен нулю). Система (4) имеет вид

= dF + Ai vi Ai v (Aj dgn + Aj dgn+i А = О dun = dun + An- vn+iAn+Q - z. dun + An+1^urJ =О.

Обратим внимание на «обмен знаками»: вместо un - vnun_і в исходной системе (1) имеем An - vn+]An+i в сопряженной системе (4). Этот феномен хорошо известен в теории сопряженных задач [7]. Если, например, vn = 1 при (n, i) Є D, а #n = КТ1 - 2un + u^1 )^h-2 при

(п, і) Є Б, то в Б задана неявная двухслойная разностная схема для уравнения теплопроводности/диффузии:

■^n+l un ________ un 2un + u

i+i

n

Т Л.2

Пусть также функционал ^ не зависит от при (п, і) Є Б (то есть от значений во внутренних узлах). Тогда система (4) примет вид

i i i 1 i i+1 An+1 — An = An — 2 An + An

т h2 '

(n, i) Є D, n < N.

Эта система также является двухслойной разностной схемой, причем также неявной, но в «обратном времени», то есть при расчете в ходе убывания n. Если формально положить Л]у +1 = 0, то указанные уравнения окажутся применимы и для n = N.

Отметим, что если функция / не зависит от параметров £ (что типично для начальных условий или граничных условий I рода), то соответствующий ЛП не входит в (5) и может не вычисляться.

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

В этом случае систему (1) запишем так:

/n G(um<n? sn? О ° (n? 0 ^ D?

sn = sn-1 + Гп (sn- Ъ un-1 ? 0 = °

Достаточно вычислить производные

d/m = dGm + dGL дГп+i / = дсп + дсп dr„

dun dsm dun ’ d£fc d£fc ds„ d£fc '

Все выкладки выше применимы и в этом частном случае зависимости от «предыстории».

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

Resume

We construct the gradient of finite-dimensional functional in the space of large dimension with constraints. The system of constraints has special structure: it is a nonlinear evolutionary grid boundary-value problem.

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

[1] Заика Ю. В., Родченкова Н. И. Моделирование высокотемпературного пика ТДС-спектра дегидрирования // Математическое моделирование. 2006. Т. 18. № 4. С. 100-112.

[2] Чернов И. А. Математическая модель экзотермичного формирования гидрида // Математическое моделирование. 2010. Т. 22. № 1. С. 3-16.

[3] Maykut G. A., Untersteiner N. Some Results from a Time-Dependent Thermodynamic Model of Sea Ice // Journal of Geophysical Research. 1971. V. 76. No 6. P. 1550-1575.

[4] Chernov I. A. Maximum Principle for the Parabolic Boundary-Value Problems with Nonlinear Boundary Conditions and Inner State // Mathematical Modelling / Christopher R. Brennan. Mathematics Research Development series. 2011. P. 203-246.

[5] Чернов И. А. Сходимость сеточно-интерполяционных аппроксимаций решения квазилинейной параболической краевой задачи на отрезке // Труды ПетрГУ. Серия Математика. 2010. Вып. 17. С. 26-37.

[6] Rodchenkova N. I., Zaika Yu. V. Numerical modelling of hydrogen desorption from cylindrical surface // International Journal of Hydrogen Energy. Elsevier. 2011. V. 36. P. 1239-1247.

[7] Марчук Г. И., Агошков В. И., Шутяев В. П. Сопряженные уравнения и методы возмущений в нелинейных задачах математической физики. М.: Физматлит, 1993. 224 с.

ИПМИ КарНЦ РАН, 185910, Петрозаводск, ул. Пушкинская, 11;

КГПА, 185910, Петрозаводск, ул. Пушкинская, 17 E-mail: [email protected], [email protected]