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

1. ВВЕДЕНИЕ

Работа посвящена численному построению квазиизометрической параметризации двумерной криволинейной области. Алгоритм основан на теории квазиизометрических отображений [1]. Разработка алгоритма и расчетно-экспериментальные исследования велись под руководством С.К. Годунова и при его тесном участии в ИПМ им. М.В. Келдыша РАН в 1995–2004 гг.

Квазиизометрические параметризации естественно использовать для решения краевых задач. Расчетные сетки, построенные с помощью квазиизометрических отображений, обладают рядом достоинств [2]. Практическое использование таких сеток наталкивается на серьезные трудности, связанные с необходимостью численной минимизации обобщенного вариационного функционала Дирихле. В данной работе рассматриваются вычислительные проблемы минимизации дискретного аналога функционала [1], отвечающего квазиизометрическому отображению плоского криволинейного четырехугольника на квадрат. Заметим, что математические проблемы, связанные с единственностью и оптимальностью квазиизометрических параметризаций, до сих пор являются нерешенными [2].

Итак, рассмотрим односвязную область $\Omega $ на плоскости и будем считать ее четырехугольником с криволинейными сторонами, которые являются достаточно гладкими (см. [1]). Четырехугольник должен быть подчинен одному дополнительному, но не очень строгому, ограничению: минимальный угол ${{\varphi }_{{\min }}} = \min \{ {{\varphi }_{1}},{{\varphi }_{2}},{{\varphi }_{3}},{{\varphi }_{4}}\} $ среди четырех его углов должен удовлетворять неравенству ${{\varphi }_{{\min }}} \geqslant 0.5\,({{\varphi }_{1}} + {{\varphi }_{2}} + {{\varphi }_{3}} + {{\varphi }_{4}} - 2\pi )$.

Напомним, что отображение $u = u(s,t)$, $v = v(s,t)$ единичного квадрата ${{K}_{0}}$ в плоскости $(s,\,\,\,t)$ на криволинейный четырехугольник $\Omega $ является квазиизометрическим, если отношение расстояния между любыми (достаточно близкими) точками $({{s}_{1}},{{t}_{1}})$ и $({{s}_{2}},{{t}_{2}})$ к расстоянию между их образами $({{u}_{1}},{{v}_{1}})$ и $({{u}_{2}},{{v}_{2}})$ ограничено сверху и снизу:

Здесь и далее, говоря об отображениях области $\Omega ,$ мы имеем в виду взаимно-однозначное отображение $u(s,t)$, ${v}(s,t)$, частные производные которого ${{u}_{s}},\;{{u}_{t}},\;{{{v}}_{s}},\;{{{v}}_{t}}$ непрерывны по Гельдеру. Квазиизометричность означает ограниченность первых производных и ограниченность снизу якобиана ${{u}_{s}}{{{v}}_{y}} - {{u}_{t}}{{{v}}_{s}} \geqslant \delta > 0$. Квазиизометрическое отображение является квазиконформным, т.е. конформным в некоторой метрике. Конформное отображение может не быть квазиизометрическим. Отображение гладкого четырехугольника на другой является конформным и квазиизометрическим, если и только если соответствующие углы и конформные модули этих четырехугольников равны. Заметим, что наилучшим квазиизометрическим отображением при заданных ограничениях будет отображение с наименьшим значением величины ${{{{\sigma }_{2}}} \mathord{\left/ {\vphantom {{{{\sigma }_{2}}} {{{\sigma }_{1}}}}} \right. \kern-0em} {{{\sigma }_{1}}}}$. Численные процедуры оценки качества квазиизометрии проведены в [3].

Постановка задачи построения квазиизометрического отображения и отвечающей этому отображению расчетной сетки на основе минимизации вариационного функционала [1] приведена в разд. 2. Краткая общая схема решения этой задачи дана в разд. 3; схема включает расчет параметров метрики, управляющих граничных функций и координат узлов сетки. Основные этапы алгоритма детализированы в разд. 4–6. Вычислительные трудности минимизации обобщенного функционала Дирихле продемонстрированы в разд. 7.

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

В [1] сформулирована задача об отыскании квазиизометрического отображения криволинейного четырехугольника $\Omega \subset {{R}^{2}}$ на единичный квадрат ${{K}_{0}} = \{ (s,t):0 \leqslant s,\;t \leqslant 1\} $. С помощью такого отображения равномерная сетка в квадрате ${{K}_{0}}$ индуцирует криволинейную четырехугольную сетку в четырехугольнике $\Omega $.

Метрика в квадрате ${{K}_{0}}$ выбирается из некоторого семейства метрик ${{h}_{{ij}}} = {{h}_{{ij}}}(s,t,\vec {s},{{x}_{0}}(s),{{x}_{1}}(s),{{y}_{0}}(t),{{y}_{1}}(t)),$ ${{h}_{{ij}}} = {{h}_{{ij}}}(s,t,\vec {s},{{x}_{0}}(s),{{x}_{1}}(s),{{y}_{0}}(t),{{y}_{1}}(t)),$ $i,j = 1,2,$ определяемых вектором параметров $\vec {s} = ({{s}_{1}},{{s}_{2}},{{s}_{3}},{{s}_{4}},k)$ и управляющими граничными функциями ${{x}_{0}}(s),$ ${{x}_{1}}(s),$ ${{y}_{0}}(t),$ ${{y}_{1}}(t)$, определенными на сторонах квадрата ${{K}_{0}}$. Параметры метрики ${{s}_{1}},\;{{s}_{2}},\;{{s}_{3}},\;{{s}_{4}},\;k$ являются аналогами углов и конформного модуля криволинейного четырехугольника.

Функции $u = u(s,t),$ ${v} = {v}(s,t)$, осуществляющие отображение ${{K}_{0}}$ на $\Omega $, ищутся из условия минимума функционала

который минимизируется по трем группам переменных:

1) параметрам метрики $\vec {s} = ({{s}_{1}},{{s}_{2}},{{s}_{3}},{{s}_{4}},k)$;

2) управляющим граничным функциям ${{x}_{0}},\;{{x}_{1}},\;{{y}_{0}},\;{{y}_{1}}$;

3) функциям $u(s,t),\;v(s,t)$.

Функции ${{x}_{0}},{{x}_{1}},{{y}_{0}},{{y}_{1}}$ задают отображение сторон квадрата ${{K}_{0}}$ на вспомогательный единичный квадрат $K = \{ (x,y),0 \leqslant x,\;y \leqslant 1\} $ (см. фиг. 1).

Отображение нижней и верхней сторон квадрата ${{K}_{0}}$ определяется функциями ${{x}_{0}}(s)$, ${{x}_{1}}(s)$, отображения левой и правой сторон – соответственно функциями ${{y}_{0}}(t)$, ${{y}_{1}}(t)$. Предполагается, что выполнены условия согласования ${{x}_{0}}(0) = {{x}_{1}}(0) = {{y}_{0}}(0) = {{y}_{1}}(0) = 0$ и ${{x}_{0}}(1) = {{x}_{1}}(1) = {{y}_{0}}(1) = {{y}_{1}}(1) = 1$, управляющие граничные функции ${{x}_{0}},\,\,{{x}_{1}},\,\,{{y}_{0}},\,\,{{y}_{1}}$ являются достаточно гладкими, монотонно возрастающими и их первые производные ограничены снизу и сверху (см. [1]).

Если функции ${{x}_{0}},\;{{x}_{1}},\;{{y}_{0}},\;{{y}_{1}}$ известны, то отображение ${{K}_{0}}$ на $K$ задается явными формулами:

Геометрически это отображение можно представить как переход от равномерной сетки на ${{K}_{0}}$ к сетке на $K$, образованной прямолинейными отрезками, которые соединяют образы соответствующих граничных точек квадрата ${{K}_{0}}$, находящиеся на противоположных сторонах квадрата $K$ (см. фиг. 1).

Отображение квадрата $K$ на криволинейный четырехугольник $\Omega $ с углами ${{\varphi }_{1}},\;{{\varphi }_{2}},\;{{\varphi }_{3}},\;{{\varphi }_{4}}$ минимизирует функционал

где метрика ${{g}_{{ij}}}$ выбирается из пятипараметрического семейства метрик, заданных в квадрате $K$ следующими формулами: (здесь $\bar {x} = x - 0.5,$ $\bar {y} = y - 0.5$).

При этом функции $u(x,y),$ $v(x,y)$ должны обеспечить отображение границы $\partial K$ квадрата $K$ на границу $\partial \,\Omega $ криволинейного четырехугольника $\Omega $, а параметры метрики ${{s}_{1}},\;{{s}_{2}},\;{{s}_{3}},\;{{s}_{4}},\;k$ удовлетворять нелинейным соотношениям:

На углы ${{\varphi }_{1}},\;{{\varphi }_{2}},\;{{\varphi }_{3}},\;{{\varphi }_{4}}$ четырехугольника $\Omega $ накладывается ограничение $0 < {{\varphi }_{j}} < \pi ,$ $2{{\varphi }_{j}} > {{\varphi }_{1}} + {{\varphi }_{2}} + {{\varphi }_{3}} + {{\varphi }_{4}} - 2\pi ,$ $j = 1,2,3,4$. Функционал (2.3) минимизируется по $u,\,v$ как функциям аргументов $x,\,\,y$ и параметрам метрики $\vec {s} = ({{s}_{1}},\,\,{{s}_{2}},\,\,{{s}_{3}},\,\,{{s}_{4}},\,\,k)$. Искомая метрика ${{g}_{{ij}}}$ представляет естественную метрику на поверхности постоянной кривизны (плоскости, сфере или плоскости Лобачевского), а геодезические в $K$ являются отрезками прямых; общая форма таких метрик описана в [4].

Поставим задачу отыскания функций $u$, ${v}$ как функций $s,\;t$. Для удобства обозначим $\varphi = {{x}_{1}} - {{x}_{0}},$ $\psi = {{y}_{1}} - {{y}_{0}},$ $w = 1 - \varphi \psi .$ Тогда

Матрица Якоби и производные этого отображения находятся по формулам: После замены переменных (2.6), (2.7) функционал (2.3) преобразуется в функционал (2.1), где или в матричной форме В (2.8) функции ${{g}_{{ij}}} = {{g}_{{ij}}}(s,t)$ определены как функции $s,\,\,t$ в соответствии с (2.4), (2.6).

Вариационная задача (2.1) требует для замыкания краевых условий на сторонах квадрата ${{K}_{0}}$. На каждой из сторон ${{K}_{0}}$ может быть задано условие Дирихле либо естественное краевое условие. Условие Дирихле означает, что на соответствующей границе области $G$ граничные точки фиксированы, а соответствующая этой границе управляющая функция (из ${{x}_{0}},\;{{x}_{1}},\;{{y}_{0}},\;{{y}_{1}}$) неизвестна. В случае естественного краевого условия граничные точки свободны, а соответствующая управляющая функция известна.

Минимальное значение функционала $\Phi = S$, где $S$ – площадь четырехугольника $\Omega $, достигается на квазиизометрическом отображении $u,v$ квадрата ${{K}_{0}}$ на $\Omega $ [1].

3. ОБЩАЯ СХЕМА АЛГОРИТМА

Задача поиска минимума функционала в формулировке (2.1) является трудной. Основные трудности стали ясны в процессе разработки алгоритма и численных экспериментов с областями различной формы. Построенная вычислительная схема является итерационной; на каждой итерации решаются следующие задачи:

1. Расчет параметров метрики $\vec {s} = ({{s}_{1}},{{s}_{2}},{{s}_{3}},{{s}_{4}},k)$; функции $u,\;v$ и ${{x}_{0}},\;{{x}_{1}},\;{{y}_{0}},\;{{y}_{1}}$ известны.

2. Расчет управляющих функций ${{x}_{0}},\;{{x}_{1}},\;{{y}_{0}},\;{{y}_{1}}$ (или эквивалентная этому шагу процедура движения точек по границам четырехугольника $\Omega $, см. разд. 6); функции $u,\;v$ и параметры метрики $\vec {s} = ({{s}_{1}},{{s}_{2}},{{s}_{3}},{{s}_{4}},k)$ известны.

3. Расчет функций $u,\;v$; считаются известными функции ${{x}_{0}},\;{{x}_{1}},\;{{y}_{0}},\;{{y}_{1}}$ и значения $u,\;v$ на границах квадрата ${{K}_{0}}$.

Приведенная схема представляет собой итерационный процесс координатного спуска; в качестве направлений спуска выступают поочередно параметры метрики, управляющие граничные функции и отображение $u(s,t),$ $v(s,t)$. Каждая из задач 1–3 является достаточно сложной и подробно описывается в следующих разделах 4–6. Для численного решения этих задач в квадрате ${{K}_{0}}$ вводится прямоугольная сетка, на которой исходный функционал заменяется разностным. Разностная дискретизация достаточно стандартна.

4. РАСЧЕТ ПАРАМЕТРОВ МЕТРИКИ

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

В данном разделе предполагается, что величины ${{G}_{{ij}}}$ известны; они вычисляются по следующим формулам (точнее по их разностным аналогам):

Заметим, что ${{u}_{x}}{{{v}}_{y}} - {{u}_{y}}{{{v}}_{x}} = {{\left( {\det J} \right)}^{{ - 1}}}({{u}_{s}}{{{v}}_{t}} - {{u}_{t}}{{{v}}_{s}})$ и

Используя (4.1)–(4.3), мы запишем вместо (2.3) эквивалентный функционал:

Для определения нового итерационного приближения параметров метрики заменим этот функционал квадратичным функционалом относительно вариаций параметров $\delta \vec {s} = (\delta {{s}_{1}},\delta {{s}_{2}},\delta {{s}_{3}},\delta {{s}_{4}},\delta k)$. В каждой точке $(x,y)$ представим ${{g}_{{ij}}}$ и ${{G}_{{ij}}}$ в виде

Тогда подынтегральное выражение в (4.4) примет вид $f = \operatorname{ch} (p - P) - \cos (\varphi - \Phi ),$ или Функция $f(p - P,\varphi - \Phi )$ является неотрицательной и выпуклой, если

На начальной стадии итерационного процесса это условие может нарушаться, поэтому для его обеспечения в каждой точке $(u,{v})$ проведем регуляризацию, положив

Всюду далее будем использовать старое обозначение $\Phi $ вместо $\hat {\Phi }$. Заменим функционал (4.4) приближенным с квадратичным интегрантом

Пусть Из (4.5) находим

Учитывая, что аргументы $(p,\varphi )$ функции $f$ являются дифференцируемыми соответствующее число раз функциями переменных ${{g}_{{11}}},\;{{g}_{{22}}},\;{{g}_{{12}}}$, приходим к следующему представлению для второго дифференциала:

где а ${{D}_{1}}$ и ${{D}_{2}}$ есть матрицы вторых производных $p$ и $\varphi $ как функций неизвестных переменных ${{g}_{{11}}},\;{{g}_{{22}}},\;{{g}_{{12}}}$: Здесь функция $z$ определена по формуле а значения ее аргументов следующие: Связь между $\vec {\delta } = (\delta p,\delta \varphi )$ и $\delta \vec {g} = (\delta {{g}_{{11}}},\delta {{g}_{{22}}},\delta {{g}_{{12}}})$ имеет вид $\vec {\delta } = R\delta \vec {g}$, $R\,$ – $\,2 \times 3$-матрица: где $d = \sqrt {{{g}_{{11}}}{{g}_{{22}}} - g_{{12}}^{2}} $.

Введем матрицы $D = {{f}_{p}}{{D}_{1}} + {{f}_{\varphi }}{{D}_{2}},$ $H = \operatorname{diag} \{ {{f}_{{pp}}},{{f}_{{\varphi \varphi }}}\} $ и перепишем (4.7) в виде

а затем в виде квадратичной формы где $\operatorname{grad} f = R{\kern 1pt} {\kern 1pt} {\text{*}}\vec {G},$ $C = R{\kern 1pt} {\text{*}}HR + D.$

Связь между $\delta \vec {g}$ и вариациями параметров метрики $\delta \vec {s} = (\delta {{s}_{1}},\delta {{s}_{2}},\delta {{s}_{3}},\delta {{s}_{4}},\delta k)$ в каждой точке $(x,y)$ квадрата $K$ имеет вид

где матрица $T$ определяется с помощью (2.4): где $\bar {x} = x - 0.5,$ $\bar {y} = y - 0.5$.

Переменные ${{g}_{{11}}},\;{{g}_{{22}}},\;{{g}_{{12}}}$ являются линейными функциями параметров метрики $({{s}_{1}},{{s}_{2}},{{s}_{3}},{{s}_{4}},k)$ (в силу (2.4)), поэтому пользуясь инвариантностью формы второго дифференциала в этом частном случае после подстановки (4.10) в (4.9), мы получаем в каждой точке $(x,y)$:

Записывая выражение (4.11) для каждой ячейки сетки в $K$ и суммируя по всем ячейкам, мы получаем квадратичную относительно $\delta \vec {s}$ аппроксимацию функционала (2.3):

Если мы находимся в достаточно малой окрестности точки минимума, то формально задача минимизации функционала (4.12) может быть легко решена в виде $\delta \vec {s} = - Q_{0}^{{ - 1}}{{G}_{0}}$. Обозначим старые значения параметров метрики через ${{\vec {s}}^{{(n)}}}$, тогда новoе итерационное приближение есть ${{\vec {s}}^{{(n + 1)}}} = {{\vec {s}}^{{(n)}}} + \delta \vec {s}$. Однако параметризация $({{s}_{1}},{{s}_{2}},{{s}_{3}},{{s}_{4}},k)$ метрики ${{g}_{{ij}}}$ является не очень удачной для численного алгоритма. Встречаются ситуации, когда малым возмущениям параметров соответствуют большие возмущения функций ${{g}_{{ij}}}$ (обычно экстремумы возмущения метрики находятся в некоторых особых точках – углах области, точке минимума функции $\theta = {{g}_{{11}}}{{g}_{{22}}} - g_{{12}}^{2}$). Это может приводить к вырождению метрики $(\theta \leqslant 0)$, либо к постановке сильных ограничений на шаге вдоль ньютоновского направления. Кроме того, на параметры метрики наложены ограничения (2.5), которые при минимизации могут не обеспечиваться. Поэтому мы введем новую параметризацию $\vec {\xi } = ({{\xi }_{1}},{{\xi }_{2}},{{\xi }_{3}},{{\xi }_{4}},{{\xi }_{5}})$, где

$\delta ({{\bar {\omega }}_{l}})$ – вариация относительного отклонения ${{\bar {\omega }}_{l}}$ $l$-го угла квадрата $K$ (в метрике ${{g}_{{ij}}}$) от прямого: ${{\bar {\omega }}_{l}} = \,\,\,\,\,{{\left( {{{\omega }_{l}} - 0.5\pi } \right)} \mathord{\left/ {\vphantom {{\left( {{{\omega }_{l}} - 0.5\pi } \right)} {\left( {0.5\pi } \right)}}} \right. \kern-0em} {\left( {0.5\pi } \right)}},$ причем согласно (2.5) ${{\bar {\omega }}_{l}}$ выражается через значения функций ${{g}_{{ij}}}$ углах области: Пятый параметр ${{\xi }_{5}}$ имеет смысл вариации “локального конформного модуля”: где $(x{\kern 1pt} {\text{*}},y{\kern 1pt} {\text{*}})$ – точка минимума функции $\theta = {{g}_{{11}}}{{g}_{{22}}} - g_{{12}}^{2}$. Преобразуем вариации (4.13) и (4.14) к виду где $d = \sqrt {{{g}_{{11}}}{{g}_{{22}}} - g_{{12}}^{2}} .$

Используя (4.10), запишем вариацию $\delta \vec {g} = (\delta {{g}_{{11}}},\delta {{g}_{{22}}},\delta {{g}_{{12}}})$ в четырех углах квадрата $(l = \overline {1,4} )$ и в точке $(x{\kern 1pt} {\text{*}},y{\kern 1pt} {\text{*}})$ через $\delta \vec {s}$ (со своей для каждой точки матрицей $T$), и с помощью (4.15) получим $\vec {\xi } = Q\,\delta \vec {s}$с некоторой $(5 \times 5)$-матрицей $Q$. Обозначим $A = {{Q}^{{ - 1}}}$, тогда $\delta \vec {s} = A\vec {\xi }.$ В результате квадратичная форма (4.12) преобразуется к виду

Обозначая

получаем

Квадратичная форма (4.16) симметрична и минимизируется при условиях

где ${{\bar {\varphi }}_{l}} = ({{\varphi }_{l}} - 0.5\pi ){\text{/}}0.5\pi $ – относительное отклонение $l$-го угла области от прямого. Условия (4.17) обеспечивают выполнение условий (2.5) (в линейном приближении).

Точка $(x{\kern 1pt} {\text{*}},y{\kern 1pt} {\text{*}}) = {\text{argmin}}\theta $, в которой задается параметр ${{\xi }_{5}}$, выбрана после серии численных экспериментов по изучению характера зависимости возмущений от положения точки $(x{\kern 1pt} {\text{*}},y{\kern 1pt} {\text{*}})$. Отметим, что направление в точку минимума функционала (4.16) уже не является ньютоновским, так как связь между переменными $({{s}_{1}},{{s}_{2}},{{s}_{3}},{{s}_{4}},k)$ и $\vec {\xi }$ является нелинейной.

Полная процедура поиска параметров метрики $\vec {s}$ является итерационной и включает итерации по нелинейности (переход от точки $\vec {s}$ к новой точке $\vec {s} + \delta \vec {s}$) и итерации при решении задачи одномерной минимизации – выборе шага спуска вдоль “ньютоновского” направления). В результате такой процедуры условия (2.5) выполняются практически точно.

5. РАСЧЕТ УПРАВЛЯЮЩИХ ГРАНИЧНЫХ ФУНКЦИЙ

Вычислим на прямоугольной сетке, заданной в квадрате ${{K}_{0}}$, сеточные аналоги величин

Для расчета управляющих граничных функций вместо функционала (2.1) рассмотрим ему эквивалентный (с переобозначением на $F$)

Функционал минимизируется по функциям ${{x}_{0}},\;{{x}_{1}},\;{{y}_{0}},\;{{y}_{1}}$, которые являются функциональными аргументами функций ${{h}_{{ij}}}$ (см.формулы (2.5)–(2.7)). Расчет пар функций ${{x}_{0}},\;{{x}_{1}}$ и ${{y}_{0}},\;{{y}_{1}}$ будем производить независимо; алгоритм расчета запишем только для пары функций ${{x}_{0}},\,\,{{x}_{1}}$. В каждой ячейке прямоугольной сетки в квадрате ${{K}_{0}}$ будем аппроксимировать подынтегральное выражение (5.1) суммой трех первых членов отрезка ряда Тейлора относительно четырех неизвестных величин ${{x}_{0}}(s),$ $x_{0}^{'}(s),$ ${{x}_{1}}(s),$ $x_{1}^{'}(s)$ в точке $s$. Сначала мы преобразуем функционал (5.1) к более удобной форме, также как и в разд. 4, с помощью введения новых переменных – неизвестных $p,\;\varphi $ и заданных $P,\;\Phi $:

Величины $p,\varphi $ являются искомыми; от множителей $\rho ,\,\,R$ функционал (5.1) не зависит. В результате подынтегральное выражение в (5.1) преобразуется к виду

Функция $f(p - P,\varphi - \Phi )$ неотрицательна и выпукла, если ${\text{|}}\varphi - \Phi {\kern 1pt} {\text{|}} \leqslant 0.5\pi .$ Для обеспечения этого условия мы проводим регуляризацию поля $(P,\Phi )$ по формулам

Аналогичная регуляризация проводилась в разд. 4, но роль $P$ и $\Phi $ там играли другие величины.

В каждой ячейке $(s,t)$ мы аппроксимируем функцию $f(p - P,\varphi - \Phi )$ тремя членами ряда Тейлора с учетом того, что аргументы $p$ и $\varphi $ являются сложными функциями переменных ${{x}_{0}},\;x_{0}^{'},\;{{x}_{1}},\;x_{1}^{'}$. Мы используем следующую суперпозицию преобразований:

При отыскании функций ${{y}_{0}},\;{{y}_{1}}$ суперпозиция имеет аналогичный вид.

В переменных $(p,\varphi )$ функционал (5.1) приближенно заменяется квадратичным с интегрантом

где а 2 × 2-матрица Гессе является диагональной: $H = \operatorname{diag} \{ \operatorname{ch} (p - P),\cos (\varphi - \Phi )\} $.

Пусть $z = ({{z}_{1}},{{z}_{2}},{{z}_{3}},{{z}_{4}}) \equiv ({{x}_{0}},x_{0}^{'},{{x}_{1}},x_{1}^{'}).$ Найдем в каждой ячейке $(s,t)$ – сетки все частные производные до второго порядка для каждого преобразования ${{Q}_{1}},\,\,{{Q}_{2}},\,\,{{Q}_{3}}$, затем с помощью стандартной достаточно громоздкой вычислительной процедуры найдем вектор градиента ${\text{grad}}\;f$ и матрицу Гессе $A$ функции $f$ как функции $z$.

В результате в каждой точке $(s,t)$ имеем аппроксимацию функционала квадратичной положительно-определенной симметричной формой; суммируя по $t$ при каждом фиксированном $s$, получаем квадратичную аппроксимацию

где ${\text{grad}}\,f$ и матрица $A$ получены суммированием по $t$ (по вертикальным столбцам ячеек) их локальных аппроксимаций.

Искомые функции ${{x}_{0}}(s)$ и ${{x}_{1}}(s)$ (их сеточные аналоги) определены в узлах сетки

В (5.2) вектор вариаций $\delta z = (\delta {{x}_{0}},\delta x_{0}^{'},\delta {{x}_{1}},\delta x_{1}^{'})$ определен в центре каждого интервала $(i,i + 1)$ сетки с помощью простейших формул:

Это преобразование линейно, поэтому квадратичная форма (5.2) инвариантна относительно него. Условия минимума квадратичной формы относительно неизвестных функций ${{\vec {\delta }}_{i}} = (\delta {{x}_{0}}(i),\;\delta {{x}_{1}}(i),\;i = 1,...,N + 1)$ имеет вид трехточечного векторного уравнения:

с однородными краевыми условиями: ${{\vec {\delta }}_{1}} = \vec {0},$ ${{\vec {\delta }}_{{N + 1}}} = \vec {0},$ где ${{\hat {A}}_{i}},\,\,{{\hat {B}}_{i}},\,\,{{\hat {C}}_{i}}$ – $2 \times 2$-матрицы. Это уравнение легко решается матричной прогонкой. После расчета $\vec {\delta }$ находим новые итерационные приближения ${{x}_{0}}(s)$ и ${{x}_{1}}(s)$ и алгоритм повторяется.

Приведенный алгоритм является методом Ньютона минимизации функционала (5.1) относительно пары функций ${{x}_{0}},\;{{x}_{1}}$ (и ${{y}_{0}},\;{{y}_{1}}$). На каждом шаге полного итерационного процесса обычно делается несколько (2–3) итераций для расчета управляющих функций.

6. РАСЧЕТ ДВИЖЕНИЯ ТОЧЕК ПО ГРАНИЦАМ ОБЛАСТИ

Пусть положение граничных точек неизвестно, например, на нижней границе четырехугольника $\Omega $ граничные точки сетки не фиксированы. Мы найдем положение этих точек из условия минимума функционала (5.1), т.е. рассмотрим случай естественного краевого условия, в такой постановке управляющая функция ${{x}_{0}}(s)$ считается известной, например, можно задать ${{x}_{0}}(s) = s$.

Пусть $l(s)$ есть функция, задающая положение точек на нижней границе; $s$ есть координата точки на границе квадрата ${{K}_{0}}$, $l$ есть координата точки на границе четырехугольника $\Omega $; в качестве $l$ удобно взять расстояние от начальной точки границы, нормированное на длину границы; тогда $0 \leqslant s \leqslant 1$, $0 \leqslant l \leqslant 1$.

Пусть ${{l}^{{(\nu )}}}$ есть предыдущее итерационное приближение. Найдем малую вариацию функции ${{x}_{0}}(s)$ из условия минимума функционала (5.1) с помощью алгоритма разд. 5 (функция ${{x}_{0}}(s)$ известна):

Тогда $\delta l = {{l}^{{(\nu + 1)}}} - {{l}^{{(\nu )}}} = l_{{{{x}_{0}}}}^{'}\delta {{x}_{0}}.$ Отсюда для каждой точки с координатой $s$ имеем:

Тем самым новое положение ${{l}^{{(\nu + 1)}}}$ граничных точек найдено; после этого выполняется этап 3 общей схемы счета – находятся функции $u(s,t)$, $v(s,t)$ с помощью решения эллиптических уравнений с краевыми условиями Дирихле, отвечающими новому положению точек ${{l}^{{(\nu + 1)}}}$. Все производные ${{u}_{s}},\;{{u}_{t}},\;{{v}_{s}},\;{{v}_{t}},\;{{x}_{s}},\;{{y}_{s}},\;{{x}_{t}},\;{{y}_{t}}$ и т.д. вычисляются как разделенные центральные разности в точке $(s,t)$ по значениям с предыдущей итерации. В приведенных формулах вариация $\delta {{x}_{0}} = x_{0}^{{new}} - {{x}_{0}}$ известна; вариации $\delta {{x}_{1}},\;\delta {{y}_{0}},\;\delta {{y}_{1}}$ вычисляются так же, как и вариация $\delta {{x}_{0}}$, если точки на соответствующей границе подвижны; в противном случае вариации полагаются равными нулю.

7. РАСЧЕТ ФУНКЦИЙ, ЗАДАЮЩИХ ОТОБРАЖЕНИЕ

Построение квазиизометрической системы координат в отдельной области основано на решении вариационной задачи (2.1). При известных управляющих функциях и параметрах метрики функции $u = u(s,t),$ $v = v(s,t)$ можно найти как решение вариационных уравнений Эйлера-Лагранжа:

Здесь $A,\;B,\;C$ – коэффициенты уравнений, заданные как функции переменных $s,t$:

Дискретную задачу получим с помощью дифференцирования по функциям $u,\;{v}$ разностного аналога функционала (2.1), используя стандартные формы для разностного представления в каждой ячейке выражений $u_{s}^{2} + v_{s}^{2},$ ${{u}_{s}}{{u}_{t}} + {{v}_{s}}{{v}_{t}},$ $u_{t}^{2} + v_{t}^{2}$. Дифференцирование дает два разностных уравнения:

которые аппроксимируют уравнения (7.1). Будем считать, что разностные уравнения записаны так, что ${{L}_{h}}$ – самосопряженный положительно-определенный оператор в пространстве сеточных функций, обращающихся в нуль на границе области.

Для решения линейных эллиптических уравнений (7.3) используем итерационный многосеточный метод [5]. Будем пользоваться вариантом метода, в котором разностная аппроксимация записывается только на исходной сетке; разностные уравнения на вспомогательных сетках формируются из исходных с помощью простых формул, фактически их комбинированием. Это обеспечивает достаточную свободу и удобства по заданию вспомогательных сеток – они образуются из исходной прямоугольной сетки удалением произвольных сеточных линий. В качестве сглаживателя в многосеточном методе используется явно-итерационная чебышёвская схема ЛИ-М [6], а для решения уравнений на самой грубой сетке используется стандартный чебышёвский метод.

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

8. ЧИСЛЕННЫЕ ЭКСПЕРИМЕНТЫ

Продемонстрируем трудности минимизации функционала (2.1) на простом примере. Рассмотрим область, показанную на фиг. 2 с построенной в ней квазиизометрической сеткой. Для наглядности сетка показана грубая, но расчеты по приведенному алгоритму проводились на подробных сетках, обеспечивающих установление метрических параметров. При минимизации функционала получены следующие значения метрических параметров: $s_{1}^{*} = s_{4}^{*} = 0,$ $s_{2}^{*} = s_{3}^{*} = - 0.893;$ $k{\text{*}} = 0.551$.

Покажем, как зависит функционал (2.1) от параметров метрики. Зафиксируем найденные остальные функциональные параметры (сетку и управляющие функции). Изучим зависимость функционала (2.1) в окрестности точки минимума как функции параметров ${{s}_{2}}$ и $k$. Напомним, что согласно (2.5) параметр метрики ${{s}_{2}}$ соответствует правому нижнему углу области (в данном случае тупому), а параметр $k$ выполняет роль “конформного модуля”. Пусть $f({{s}_{2}},k) = \Phi - S$, где $S$– площадь области. Тогда точное значение функции $f$ в точке минимума $f(s_{2}^{*},k{\text{*}}) = 0$. Профили $f({{s}_{2}},k{\text{*}})$ и $f(s_{2}^{*},k)$ показаны на фиг. 3а и 3б соответственно. Видно, что в окрестности точки минимума метрика вырождается (в точках $({{s}_{2}},k)$, где метрика не существует, мы полагаем $f = 10$). Эти рисунки иллюстрируют вычислительные трудности, с которыми мы сталкиваемся в процессе разработки алгоритма для решения проблемы минимизации обобщенного функционала Дирихле.

Численные эксперименты с областями различной формы проведены в достаточно большом количестве. С.К. Годунов высказывал предложение о создании атласа областей с указанием оптимальных метрических параметров, отвечающих наилучшему квазиизометрическому отображению. Эта задача до сих пор является актуальной.

9. ЗАКЛЮЧЕНИЕ

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