Журнал вычислительной математики и математической физики, 2020, T. 60, № 1, стр. 18-28

ВВЕДЕНИЕ

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

По своей сути эта методология представляет собой метод множителей Лагранжа, который применяется к дискретному варианту задачи оптимального управления. Рассматриваемая задача оптимального управления в этом случае формулируется в конечномерном пространстве. Целевой функционал и связи, наложенные на управления ($u \in {{R}^{r}}$) и фазовые переменные ($z \in {{R}^{n}}$), аппроксимируются, в результате чего целевому функционалу ставится в соответствие некоторая функция конечного числа переменных $W(z(u),u)$, а связям – система, состоящая из $n$ нелинейных алгебраических уравнений $\Phi (z,u) = {{0}_{n}}$ (${{0}_{n}}$ – это $n$-мерный нулевой вектор).

Уравнения связи $\Phi (z,u) = {{0}_{n}}$ всегда могут быть записаны в виде многошагового процесса ${{z}_{i}} = F(i,{{Z}_{i}},{{U}_{i}})$, $1 \leqslant i \leqslant n$, где ${{Z}_{i}}$ – набор векторов ${{z}_{j}}$, которые входят в правую часть последнего равенства, а ${{U}_{i}}$ – набор векторов ${{u}_{j}}$, которые входят в правую часть этого равенства. Предполагается, что $W(z,u)$ – непрерывно дифференцируемая скалярная функция, а функции $F(i,{{Z}_{i}},{{U}_{i}})$ дифференцируемы по всем компонентам векторов $z$ и $u$. БАД-методология предлагает канонические формулы, позволяющие точно вычислять градиент сложной функции $\Omega (u) = W(z(u),u)$:

Множители ${{p}_{i}} \in {{R}^{n}}$ определяются из следующей системы линейных алгебраических уравнений (дискретная сопряженная задача):

БАД-методология позволяет выбрать такую аппроксимацию сопряженной задачи, которая согласована с дискретизацией прямой задачи и с аппроксимацией целевого функционала. Эта методология значительно превосходит символьное дифференцирование и вычисление производных с помощью метода конечных разностей.

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

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

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

1. ПОСТАНОВКА ЗАДАЧИ

Рассмотрим ограниченную область $Q \subset {{R}^{n}}$ с кусочно-гладкой границей $S = \partial Q$. Эта область заполнена исследуемым веществом. Распределение температурного поля в каждый момент времени описывается следующей начально-краевой (смешанной) задачей:

Здесь $x = ({{x}_{1}},...,{{x}_{n}})$ – декартовы координаты; $t$ – время; $T(x,t)$ – температура материала в точке с координатами $x$ в момент времени $t$; $C(x)$ – объемная теплоемкость материала; $K(T)$ – коэффициент теплопроводности; ${{w}_{0}}(x)$ – заданная температура в начальный момент времени $t = 0$; ${{w}_{S}}(x,t)$ – заданная температура на границе объекта. Объемная теплоемкость вещества $C(x)$ рассматривается как известная функция координат.

Если зависимость коэффициента теплопроводности $K(T)$ от температуры $T$ известна, то, решив смешанную задачу (1.1)–(1.3), найдем распределение температуры $T(x,t)$ в $G = Q \times [0,\Theta ]$. Задачу (1.1)–(1.3) ниже будем называть прямой задачей.

Обратная коэффициентная задача сводится к следующей вариационной задаче: требуется найти такую зависимость коэффициента теплопроводности вещества $K(T)$ от температуры $T$, при которой температурное поле $T(x,t)$ и поток тепла $\left( { - K(T(x,t))\frac{{\partial T(x,t)}}{{\partial \,\bar {n}}}} \right)$ на границе объекта, полученные в результате решения прямой задачи (1.1)–(1.3), мало отличаются от температурного поля $Y(x,t)$ и потока тепла $P(x,t)$, полученных экспериментально. Мерой отклонения этих функций может служить величина:

Здесь $\varepsilon \geqslant 0$, $\beta (x,t) \geqslant 0$, $\mu (x,t) \geqslant 0$ – заданные весовые параметры; $Y(x,t)$ – заданное температурное поле; $P(x,t)$ – заданный тепловой поток на границе $S$ области $Q$; $\partial T{\text{/}}\partial{ \bar {n}}$ является производной от температуры по направлению внешней нормали к границе области; $[a,b]$ – интервал температур, где функция $K(T)$ будет восстановлена. Последнее слагаемое в функционале (1.4) используется для “сглаживания” решения вариационной задачи с умеренными возмущениями экспериментальных данных.

Таким образом, задача идентификации функции $K(T)$ сводится к следующей задаче оптимального управления: определить оптимальное управление $K(T)$ и соответствующее решение $T(x,t)$ задачи (1.1)–(1.3), при котором целевой функционал (1.4) достигает минимального значения.

2. АНАЛИТИЧЕСКОЕ ВЫРАЖЕНИЕ ДЛЯ ГРАДИЕНТА ФУНКЦИОНАЛА

При определении градиента функционала (1.4) в непрерывном случае для простоты изложения потребуем, чтобы функции ${{w}_{0}}(x)$, ${{w}_{S}}(x,t)$, $C(x)$ (входные данные задачи), а также искомая функция $K(T)$ были такими, при которых решение задачи (1.1)–(1.3) принадлежит классу $T(x,t) \in {{С}^{2}}(G) \cap {{C}^{1}}(\overline G )$.

Уравнение (1.1) может быть записано в следующем дивергентном виде:

где $E(T) \in {{C}^{1}}\left( {[a,b]} \right)$ – удельная внутренняя энергия вещества.

Класс $\Gamma $ допустимых управляющих функций определяется множеством:

Пусть $K(z) \in \Gamma $ – некоторая допустимая функция, а $T(x,t)$ – соответствующее ей решение начально-краевой задачи (1.1)–(1.3). Рассмотрим произвольную функцию $\tilde {K}(z) \in \Gamma $, мало в норме пространства ${{C}^{1}}([a,b])$ отличающуюся от функции $K(z)$. С помощью соотношений (1.1)–(1.3) она порождает функцию $\tilde {T}(x,t)$. Обозначим через ${{\delta }_{z}}K(z) = \tilde {K}(z) - K(z)$, $z \in [a,b]$, вариацию управляющей функции $K(z)$, а через ${{\delta }_{{xt}}}T(x,t)$ – вариацию в точке $(x,t)$ фазовой переменной $T(x,t)$, т.е. ${{\delta }_{{xt}}}T(x,t) = \tilde {T}(x,t) - T(x,t)$, $(x,t) \in \bar {G}$.

Запишем функционал Лагранжа

где $p(x,t) \in {{С}^{2}}(G) \cap {{C}^{1}}(\bar {G})$ – произвольная функция. Используя общий метод множителей Лагранжа, первую вариацию исследуемого функционала, вызванную вариацией ${{\delta }_{z}}K(z)$ управляющей функции $K(z)$, можно записать в следующем виде: где $\Psi (x,t) = - K(T(x,t))\frac{{\partial T(x,t)}}{{\partial{ \bar {n}}}}$.

Если функцию $p(x,t)\,$ определить из следующих условий (2.1)–(2.3) (сопряженная задача):

то получим следующее выражение для первой вариации исследуемого функционала:

Преобразуем выражение (2.4). Разобьем область $G$ на ряд подобластей ${{G}_{i}}$ так, что $G = \bigcup\nolimits_i^{} {{{{\bar {G}}}_{i}}} $, причем ${{G}_{i}} \cap {{G}_{j}} = \emptyset $ при $i \ne j$. Также потребуем от разбиения, чтобы в каждой из подобластей ${{G}_{i}}$ не содержались две различные изотермические поверхности, отвечающие одной и той же температуре.

В каждой из подобластей ${{G}_{i}}$ введем новые независимые переменные $(\xi ,\tau ) = ({{\xi }_{1}},...,{{\xi }_{n}},\tau ) \in {{q}_{i}}$, так что ${{x}_{1}} = x_{{_{1}}}^{i}({{\xi }_{1}},...,{{\xi }_{n}},\tau )$, …, ${{x}_{n}} = x_{n}^{i}({{\xi }_{1}},...,{{\xi }_{n}},\tau )$, $t = {{t}^{{_{i}}}}({{\xi }_{1}},...,{{\xi }_{n}},\tau )$. Поверхности ${{\xi }_{n}} = {\text{const}}$ – изоповерхности поля температур. В качестве координатных поверхностей ${{\xi }_{1}} = {\text{const}}$, …, ${{\xi }_{{n - 1}}} = {\text{const}}$ и $\tau = {\text{const}}$ можно выбрать любые семейства поверхностей, сохраняющие в области ${{G}_{i}}$ невырожденность вводимого преобразования. В новых переменных справедливо соотношение

дифференцируя которое по переменным ${{\xi }_{1}},...,{{\xi }_{n}}$ и $\tau $, получаем систему линейных алгебраических уравнений относительно градиента температуры в ${{R}^{{n + 1}}}$где для краткости записи используются обозначения $\varepsilon = (\xi ,\tau ) = ({{\xi }_{1}},...,{{\xi }_{n}},\tau )$ и $x(\varepsilon ) = \left( {{{x}_{1}}(\varepsilon ),...,{{x}_{n}}(\varepsilon )} \right)$.

Введем в рассмотрение вектор $\bar {N}({{x}^{i}}(\varepsilon ),{{t}^{i}}(\varepsilon )) = \left( {{{N}_{1}},\,...,\,{{N}_{n}}} \right)$ с компонентами ${{N}_{k}}\left( {{{x}^{i}}(\varepsilon ),{{t}^{i}}(\varepsilon )} \right) = {{\Delta }^{k}}$, $k = 1,\,...,\,n$. Здесь ${{\Delta }^{k}}$ – детерминант матрицы, получаемой из основной матрицы системы уравнений (2.5) заменой ее $k$-го столбца $(n + 1)$-мерным столбцом ${{\left( {0,\,0,\,...\,,\,0,\,1,\,0} \right)}^{{\text{Т}}}}$ (столбцом свободных членов системы уравнений (2.5)).

Если обозначить через

якобиан преобразования в подобласти ${{G}_{i}}$, то, как следует из системы уравнений (2.5), имеем

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

Если обозначить интервал изменения температур в подобласти ${{G}_{i}}$ как $(\eta _{{Ini}}^{{(i)}},\eta _{{Fin}}^{{(i)}})$, а прообраз точек поверхности $\eta = {\text{const}} \in (\eta _{{Ini}}^{{(i)}},\eta _{{Fin}}^{{(i)}})$, содержащихся в ${{G}_{i}}$, как ${{\Omega }_{i}}$, то выражению для первой вариации функционала можно придать вид

где ${\rm M}({{\xi }_{n}}) = \sum\limits_i {\int\limits_{{{\Omega }_{i}}} {[{\text{sign}}\left( {{{J}_{i}}(\varepsilon )} \right)({{\nabla }_{x}}p({{x}^{i}}(\varepsilon ),{{t}^{i}}(\varepsilon )),\bar {N}({{x}^{i}}(\varepsilon ),{{t}^{i}}(\varepsilon )))]d{{\xi }_{1}}...d{{\xi }_{{n - 1}}}d\tau } } $.

Как хорошо известно, градиент целевого функционала (1.4) определяется в виде

3. ЧИСЛЕННОЕ РЕШЕНИЕ ЗАДАЧИ ОПТИМАЛЬНОГО УПРАВЛЕНИЯ

Сформулированная задача оптимального управления решалась численно. Отрезок $[a,b]$, на котором восстанавливается функция $K(T)$, $T \in [a,b]$, определялся как множество значений функций ${{w}_{0}}(x)$ и ${{w}_{S}}(x,t)$, т.е. границы отрезка $[a,b]$ задавались как минимальное и максимальное значения указанных функций. Этот отрезок разбивался точками ${{\tilde {T}}_{0}} = a$, ${{\tilde {T}}_{1}}$, ${{\tilde {T}}_{2}}$, …, ${{\tilde {T}}_{M}} = b$ на $M$ частей (части могут быть как равными, так и неравными). Каждой из точек ${{\tilde {T}}_{m}}$, $m = 0,...,M$, ставилось в соответствие число ${{k}_{m}} = K({{\tilde {T}}_{m}})$. Искомая функция $K(T)$ аппроксимировалась непрерывной кусочно-линейной функцией с узлами в точках $\left\{ {\,({{{\tilde {T}}}_{m}},\,\,{{k}_{m}})\,} \right\}_{{m = 0}}^{M}$, так что

Если в процессе решения задачи температура в точке выходила за границы отрезка $[\,a,b]$, то для определения функции $K(T)$ использовалась линейная экстраполяция.

Целевой  функционал  (1.4)  аппроксимировался  функцией $F({{k}_{0}},{{k}_{1}},...,{{k}_{M}})$   конечного числа переменных с помощью метода прямоугольников. Аппроксимация выражения $\left( { - K(T(x,t))\frac{{\partial T(x,t)}}{{\partial \,\bar {n}}}} \right)$, с помощью которого вычисляется поток тепла на границе объекта, представлена в работах [8] (одномерный случай) и [9] (двумерный случай).

Задачи оптимального управления, подобные этой, обычно решаются численно с помощью некоторого метода спуска, который требует знания градиента функционала (1.4), и при этом крайне важно определять точное значение этого градиента. Использовать аналитическое выражение (2.8) для градиента функционала при численном решении задачи практически невозможно в силу его сложности. Поэтому в предлагаемом алгоритме для вычисления градиента функционала использовалась эффективная методология быстрого автоматического дифференцирования. Вычисление градиента функционала с ее помощью требует затрат машинного времени, не превосходящих утроенных затрат на вычисление самой функции. Эффективность этой методологии при вычислении градиента функции обеспечивается использованием решения сопряженной задачи. Построение такой аппроксимации сопряженной задачи, которая согласована с выбранной аппроксимацией прямой задачи, и есть основное достоинство БАД-методологии. БАД-методология поставляет канонические формулы, вычисленный с помощью которых градиент целевой функции является точным для выбранной аппроксимации задачи оптимального управления.

Одним из основных элементов предлагаемого алгоритма решения обратной коэффициентной задачи является решение прямой задачи. Для численного решения прямой задачи вводилась временнáя и пространственная сетки (в общем случае, неравномерные). В каждом сеточном узле расчетной области $\bar {Q} \times [0,\Theta ]$ все функции определялись своими точечными значениями. Для аппроксимации уравнения теплопроводности использовалась некоторая подходящая конечно-разностная схема.

В работе [8] рассмотрена задача идентификации коэффициента теплопроводности в одномерном случае. Здесь обратная коэффициентная задача сводилась к следующей вариационной задаче: найти такую зависимость коэффициента теплопроводности от температуры, при которой температурное поле и поток тепла на левой границе образца (левый торец), полученные в результате решения прямой задачи, мало отличаются от данных, полученных экспериментально. В этой работе предполагалось, что коэффициент теплопроводности непрерывно зависит от температуры (отсутствуют скачки коэффициента в окрестности какой-нибудь точки).

При численном решении обратной задачи использовался разработанный авторами алгоритм. Было проведено большое число расчетов, которые продемонстрировали эффективность этого алгоритма. Тестовые примеры показали, что с помощью предложенного алгоритма коэффициент теплопроводности вещества воспроизводился с высокой точностью (относительная погрешность не превышала ${{10}^{{ - 10}}}$). Рассмотренные примеры также выявили некоторые особенности рассматриваемой обратной задачи и предложенного алгоритма.

Так, при использовании функционала “поле” ($\mu (x,t) \equiv 1,\quad \beta (x,t) \equiv 0,\quad \varepsilon = 0$) в некоторых случаях приходится сталкиваться с неединственностью решения. Будет решение обратной задачи единственно или нет – существенно зависит от заданного экспериментального поля $Y(x,t)$. В работе [8] приводятся необходимые условия на поле $Y(x,t)$, при которых возникает неединственность решения обратной задачи для одномерного уравнения теплопроводности.

Рекомендуется решать задачу идентификации несколько раз (3–4), выбирая в качестве начального приближения каждый раз новую функцию. Если при этом получается одно и то же решение задачи оптимального управления, то его можно брать в качестве решения задачи идентификации. Если же решение задачи оптимального управления будет зависеть от выбора начального приближения, то для решения задачи идентификации необходимо задать дополнительное условие, например, задавать точку, в которой искомый коэффициент теплопроводности известен. Важно отметить, что если для решения этой задачи использовать функционал “поток” с $\beta (x,t) > 0$, то всегда получалось единственное решение (даже при $\mu (x,t) \equiv 0$).

Предложенный алгоритм эффективно работает и при восстановлении разрывного коэффициента теплопроводности. Обратная задача для одномерного нестационарного уравнения теплопроводности в случае, когда коэффициент теплопроводности имел скачки I рода, рассмотрена в [10].

Исследованию сформулированной обратной коэффициентной задачи в двумерном случае посвящены работы [9], [11]–[14]. Как и предполагалось, двумерный случай в сравнении с одномерным потребовал больших затрат машинных ресурсов (как по памяти, так и по времени). Здесь также обнаружилось появление неединственности решения в случае, когда искомая функция восстанавливается только по заданному температурному полю (функционал – “поле”). Необходимые условия на поле $Y(x,t)$, при которых возникает неединственность решения обратной задачи для двумерного уравнения теплопроводности, приводятся в [12]. Вместе с тем в двумерном случае проявились новые особенности обратной коэффициентной задачи и алгоритма ее решения.

Во-первых, проведенные исследования позволили заключить, что если в целевом функционале используется тепловой поток не на всей границе области ($\beta (x,t) \equiv 0$ на части границы), то коэффициент теплопроводности может идентифицироваться не полностью (см. [13]). На основе полученных результатов сформулирована рекомендация учитывать в целевом функционале экспериментальный тепловой поток на всей границе области.

Похожая ситуация наблюдается и в случае использования функционала “поле”, когда температурное поле задано только на некотором подмножестве рассматриваемой области (см. [12]). Объяснить это явление можно тем, что из-за выбора весовой функции целевой функционал не зависит от значений коэффициента теплопроводности для всех значений интервала температур $[a,b]$.

Во-вторых, отметим еще одну особенность обратной коэффициентной задачи, с которой пришлось сталкиваться при решении задачи в многомерном случае. Здесь градиент целевого функционала распределен по отрезку $[a,b]$ неравномерно: значения его модуля вблизи концов отрезка близки к нулю и в разы (а то и на порядки) меньше значений в остальной части отрезка. Указанное отличие возрастает с ростом числа $M$ разбиений отрезка $[\,a,b]$ на подотрезки. Объяснить отмеченное явление можно тем, что градиент функционала вычисляется интегрированием некоторой функции по множеству точек $(x,t) \in {{Q}_{m}} \in Q \times (0,\Theta )$, $m = 1,...,M$, в котором температура $T$ принадлежит подотрезку $[{{\tilde {T}}_{{m - 1}}},{{\tilde {T}}_{m}}]$, $m = 1,...,M$. Максимальное и минимальное значения температуры часто достигаются лишь в одной точке, поэтому мера множеств ${{Q}_{1}}$ и ${{Q}_{M}}$ близка к нулю и заметно меньше мер остальных множеств ${{Q}_{m}}$. Неравномерное распределение модуля градиента сильно затрудняет получение численного решения задачи.

Эффективно устранить описанную трудность при решении задачи позволил подход, основанный на последовательном увеличении числа $M$ разбиений отрезка $[\,a,b]$. Начинать процесс желательно с $M = 1$. После получения оптимального решения его следует использовать в качестве начального приближения для варианта с $M = 2$. Получив оптимальное решение для $M = 2$, использовать его в качестве начального приближения для $M = 4$, и т.д. Описанный подход активно использовался при решении рассматриваемой обратной задачи. Отметим, что в одномерном случае это явление также присутствует, но гораздо слабее, так что на него можно было не обращать внимание.

В-третьих, возник вопрос, связанный с аппроксимацией прямой задачи. При рассмотрении обратной задачи в одномерном случае для аппроксимации уравнения теплопроводности использовалась двухслойная неявная схема с весами. В силу нелинейности получающейся системы алгебраических уравнений для ее решения строился некоторый итерационный процесс (см. [8]). Для аппроксимации сопряженной задачи БАД-методология в этом случае также “строит” неявную схему с весами. Однако получающаяся при этом система алгебраических уравнений является линейной, и для ее разрешения применялся эффективный метод прогонки.

В двумерном случае ситуация усложняется. При аппроксимации прямой задачи неявной схемой с весами аппроксимация сопряженной задачи будет представлять собой линейную систему алгебраических уравнений, для разрешения которой (в силу двумерности по пространству) также требовалось строить итерационный процесс. В связи с этим возникла идея использовать для аппроксимации прямой задачи какую-либо экономичную разностную схему (см. [15]). В работах [9], [12], [13] для этой цели была выбрана схема переменных направлений. В этом случае для решения сопряженной задачи не требуется итерационного процесса, линейная система алгебраических уравнений расщепляется на ряд подсистем, каждая из которых решается методом прогонки [15].

В [14] проведен численный анализ устойчивости решений, полученных в [8]–[13]. Было показано, что возмущение восстановленного коэффициента теплопроводности имеет тот же порядок, что и вызвавшее его малое возмущение в экспериментальных данных.

4. О РЕШЕНИИ ОБРАТНОЙ КОЭФФИЦИЕНТНОЙ ЗАДАЧИ В ТРЕХМЕРНОЙ ПОСТАНОВКЕ

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

Пусть образец исследуемого вещества имеет форму параллелепипеда длины $R$, ширины $V$ и высоты $W$. Для численного решения прямой задачи вводятся временнáя и пространственная сетки (в общем случае, неравномерные).

Временнáя сетка задается набором узловых значений $\{ {{t}^{j}}\} _{{j = 0}}^{J}$, ${{t}^{0}} = 0$, ${{t}^{J}} = \Theta $. Шаги ${{\tau }^{j}}$ этой сетки определяются соотношениями ${{\tau }^{j}} = {{t}^{{j + 1}}} - {{t}^{j}}$, $j = \overline {0,J - 1} $.

Вводятся также две пространственные сетки: основная и вспомогательная. Для построения основной пространственной сетки на отрезке $[0,R]$ выбирается система “опорных” точек $\left\{ {{{x}_{n}}} \right\}_{{n = 0}}^{N}$ так, что ${{x}_{0}} = 0$, ${{x}_{N}} = R$ и ${{x}_{n}} < {{x}_{{n + 1}}}$ для всех $0 \leqslant n < N$. При этом $h_{n}^{x}$ – расстояние между опорными точками ${{x}_{n}}$ и ${{x}_{{n + 1}}}$, т.е. $h_{n}^{x} = {{x}_{{n + 1}}} - {{x}_{n}}$, $n = \overline {0,N - 1} $. Аналогично на отрезках $[0,V]$ и $[0,W]$ выбирается система “опорных” точек $\left\{ {{{y}_{i}}} \right\}_{{i = 0}}^{I}$ и $\left\{ {{{z}_{l}}} \right\}_{{l = 0}}^{L}$ так, что ${{y}_{0}} = 0$, ${{y}_{I}} = V$, ${{y}_{i}} < {{y}_{{i + 1}}}$ для всех $0 \leqslant i < I$ и ${{z}_{0}} = 0$, ${{z}_{L}} = W$, ${{z}_{l}} < {{z}_{{l + 1}}}$ для всех $0 \leqslant l < L$ соответственно. При этом $h_{i}^{y} = {{y}_{{i + 1}}} - {{y}_{i}}$, $i = \overline {0,I - 1} $, $h_{l}^{z} = {{z}_{{l + 1}}} - {{z}_{l}}$, $l = \overline {0,L - 1} $. Точки основной сетки – совокупность точек $\left\{ {{{x}_{n}},{{y}_{i}},{{z}_{l}}} \right\}$, где $n = \overline {0,N} $, $i = \overline {0,I} $ и $l = \overline {0,L} $.

Вспомогательная сетка $\left\{ {{{{\tilde {x}}}_{n}},{{{\tilde {y}}}_{i}},{{{\tilde {z}}}_{l}}} \right\}$, $n = \overline {0,N + 1} $, $i = \overline {0,I + 1} $, $l = \overline {0,L + 1} $ строится похожим образом. Отличие состоит лишь в выборе опорных точек:

Линии $x = {{\tilde {x}}_{n}}$, $y = {{\tilde {y}}_{i}}$ и $z = {{\tilde {z}}_{l}}$ делят область $Q$ на так называемые расчетные ячейки. Расчетной ячейке будем приписывать индексы $\left( {n,\;i,\;l} \right)$, если вершина этой ячейки, наиближайшая к началу координат, совпадает с узловой точкой $\left( {{{{\tilde {x}}}_{n}},\;{{{\tilde {y}}}_{i}},\;{{{\tilde {z}}}_{l}}} \right)$, а ее объем будем обозначать через ${{V}_{{nil}}}$, т.е.

Обозначим также среднюю температуру в ячейке с индексами $\left( {n,i,l} \right)$ через ${{T}_{{nil}}}\left( t \right)$, а через $T_{{nil}}^{j}$ – ее значение в момент времени ${{t}^{j}}$.

При исследовании обратной коэффициентной задачи для трехмерного уравнения теплопроводности для аппроксимации начально-краевой задачи предлагается использовать один из вариантов схемы переменных направлений (схемы расщепления по направлению $x$, $y$ и $z$). В этом случае, как показал опыт решения двумерной обратной задачи, для решения сопряженных уравнений, полученных с помощью БАД-методологии на основе дискретной прямой задачи, не потребуется построения итерационного процесса.

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

Наряду с основными значениями искомой сеточной функции $T_{{nil}}^{j}$ введем еще два промежуточных значения $T_{{nil}}^{{j + 1/3}}$ и $T_{{nil}}^{{j + 2/3}}$ функции температуры, которые можно формально рассматривать как значения температуры в момент времени ${{t}^{j}} + {{\tau }^{j}}{\text{/}}3$ и ${{t}^{j}} + 2{{\tau }^{j}}{\text{/}}3$ (см. [15]). Для краткости написания рассматриваемых ниже схем введем следующие обозначения:

Простейшая схема расщепления для трехмерного уравнения теплопроводности – это локально-одномерная схема (см. [15]). Она имеет следующий вид:

$x$ – направление

$y$ – направление

$z$ – направление

Как отмечено в [15], при дополнительных предположениях относительно ограниченности вторых производных, указанную выше схему можно рассматривать в двух вариантах:

1) $p = j + 1{\text{/}}3$, $q = j + 2{\text{/}}3$, $r = j + 1$;

2) $p = j$, $q = j + 1{\text{/}}3$, $r = j + 2{\text{/}}3$.

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

Одна из наиболее известных и распространенных схем переменных направлений для трехмерного уравнения теплопроводности – это разностная схема Дугласа–Рекфорда (см. [16]). Ее называют еще схемой стабилизирующей поправки (см. [17]). Она имеет следующий вид:

$x$ – направление

$y$ – направление

$z$ – направление

Если использовать $p = j$ (что часто и делают), то на каждом этапе получаются линейные системы уравнений с трехдиагональной матрицей.

Следующая рекомендуемая схема для решения начально-краевой задачи (1.1)–(1.3) – это схема Писмена–Рекфорда:

$x$ – направление

$y$ – направление

$z$ – направление

Системы нелинейных алгебраических уравнений, которые получаются на первом этапе дискретизации прямой задачи, решаются относительно неизвестных значений $T_{{nil}}^{{j + 1/3}}$. Эта система уравнений (4.1) разбивается на ряд подсистем. Каждая подсистема характеризуется фиксированными значениями индексов $i$ и $l$ и решается независимо от других подсистем. Эти подсистемы решаются итерационно с привлечением метода прогонки. Коэффициенты подсистем определяются по известным значениям функции $T_{{nil}}^{{j + 1/3}}$ на данной итерации. С помощью метода прогонки определяется значение функции $T_{{nil}}^{{j + 1/3}}$ на последующей итерации. Итерации прекращаются при достижении требуемой точности.

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

Насколько нам известно, до настоящего времени не было представлено доказательство численной устойчивости и сходимости указанных выше схем для нелинейных задач. В то же время было замечено, что численные эксперименты с методом переменных направлений свидетельствуют в пользу успешности таких применений, хотя строгого обоснования этому до сих пор не было построено. Этот факт отмечается как в работах американских ученых, например, в работах Дж. Дугласа, так и в работах А.А. Самарского и его учеников. Авторы этих схем рекомендуют использовать несколько разностных схем при решении трехмерных уравнений теплопроводности с зависящими от температуры коэффициентами. При этом необходимо проводить исследования по выбору расчетных сеток. Заметим, что в работе [18] показано, что для трехмерной задачи схема Писмена–Рекфорда является условно устойчивой даже в случае, когда коэффициенты уравнения не зависят от температуры.

Выбор какой-нибудь из указанных схем для решения поставленной обратной коэффициентной задачи можно считать обоснованным лишь в случае, когда результаты решения задачи с использованием нескольких схем практически совпадают. Этот критерий успешно использовался авторами настоящей статьи при выполнении цикла работ, в которых исследовалась задача управления процессом кристаллизации металла в литейном деле (см. [4]–[7]). Для дискретизации задачи была выбрана схема Писмена–Рекфорда. Она сравнивалась с локально-одномерной схемой и с двухслойной неявной схемой с весами (см. [19]).

ЗАКЛЮЧЕНИЕ

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