Журнал вычислительной математики и математической физики, 2023, T. 63, № 11, стр. 1839-1848

1. ВВЕДЕНИЕ

При численном решении дифференциальных уравнений наиболее широко используются конечно-разностный и спектральные методы (см. [1, 2]). Конечно-разностный метод прост в программной реализации и экономичен (см. [3, 4]). Спектральный метод более удобен в задачах частотного зондирования Земли (см. [5]), дефектоскопии и др.

Важным показателем качества численного алгоритма является точность получаемого решения уравнений (см. [6]). Повышению точности решений посвящено большое число работ (см. [7]–[9]). Существуют различные способы повышения точности решения, например, построение компактных схем (см. [8]), использование разностных схем более высокого порядка аппроксимации (см. [9]) или построение разностных схем, минимизирующих погрешность дисперсионного соотношения (dispersion-relation-preserving) (см. [10]). К последнему типу схем относятся так называемые оптимальные разностные схемы.

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

В [11] была предложена оптимальная разностная схема для решения волнового уравнения в спектральной области. В разностное уравнение 2-го порядка аппроксимации для заданной гармоники Фурье вводятся три дополнительных параметра. Значения этих параметров определяются минимизацией погрешности численного решения на точном аналитическом решении. Алгоритм рассматривается при равных пространственных шагах разностной сетки. Обобщение для неравных шагов было предложено в [12] введением средних значений в пространственные производные. В этом случае оптимизация проводилась по четырем параметрам.

В [13, 14] такая методика построения оптимальных разностных схем была применена при решении уравнений с использованием спектрального преобразования Лагерра.

В ряде задач спектральный метод, основанный на преобразовании Лагерра, по эффективности в несколько раз превосходит метод Фурье. Эффективность обусловлена видом уравнений для гармоник Лагерра (см. [13–15]). В частности, система уравнений для гармоник Лагерра всегда содержит только действительные переменные.

В [13, 14] была рассмотрена оптимальная разностная схема 2 -го порядка аппроксимации для решения уравнений Максвелла и для решения волнового уравнения. В разностной схеме для гармоник использовались оптимальные параметры. Решение уравнений получалось обратным преобразованием Лагерра.

Такая методика получения оптимальных параметров легко распространяется на разностные схемы 4 -го порядка аппроксимации (см. [15]).

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

В разностное уравнение 2-го порядка аппроксимации для гармоник Лагерра вводятся дополнительные параметры. Оптимальные значения этих параметров определяются минимизацией погрешности численного решения на точном аналитическом решении уравнений для гармоник.

Применяя обратное преобразование Лагерра к уравнениям для гармоник, получают дифференциально-разностные уравнения. Эти уравнения являются разностными по пространству и дифференциальными по времени. Эти уравнения решаются с использованием конечно-разностных схем 2 -го порядка аппроксимации и по времени, и по пространству.

Предлагаемые конечно-разностные оптимальные схемы дают более точное решение в сравнении с неоптимальными конечно-разностными схемами 2 -го порядка аппроксимации.

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

Приводятся результаты для одномерных и двумерных уравнений Максвелла.

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

Система уравнений Максвелла в 2-мерном случае состоит из двух независимых систем уравнений (см. [1, 16]). Будем рассматривать одну (TE-мода) из них:

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

Здесь ${\mathbf{H}} = ({{H}_{x}},{{H}_{y}},{{H}_{z}})$ – напряженность магнитного поля, ${\mathbf{E}} = ({{E}_{x}},{{E}_{y}},{{E}_{z}})$ – напряженность электрического поля, ${\mathbf{J}} = ({{J}_{x}},{{J}_{y}},{{J}_{z}})$ – ток внешнего источника, являющийся функцией времени. В эту систему уравнений входит только одна компонента вектора тока ${{J}_{y}}$, две другие компоненты тока ${{J}_{x}}$, ${{J}_{z}}$ входят в другую независимую систему уравнений.

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

Здесь введены эффективные значения диэлектрической проницаемости ε и проводимости $\sigma $:

где ${{\varepsilon }_{s}}$ – статическое значение диэлектрической проницаемости, ${{\tau }_{E}}$ – время релаксации напряженности электрического поля, ${{\tau }_{D}}$ – время релаксации электрической индукции, ${{\sigma }_{s}}$ – статическое значение проводимости, ${{\tau }_{\sigma }}$ – время релаксации проводимости (см. [1, 2]).

В частном случае малых ${{\tau }_{E}}$, ${{\tau }_{D}}$, ${{\tau }_{\sigma }}$и ${{{{\tau }_{E}}} \mathord{\left/ {\vphantom {{{{\tau }_{E}}} {{{\tau }_{D}}}}} \right. \kern-0em} {{{\tau }_{D}}}} \simeq 1$ эта система уравнений принимает стандартный вид уравнений Максвелла (см. [16]).

Проведем преобразование Лагерра (см. [17]) по времени уравнений (1):

где $l_{n}^{\alpha }(ht)$ – ортогональная функция Лагерра (см. [17]) степени $n$, $\alpha $ – целая константа, $h$ – параметр преобразования Лагерра, который далее положим равным единице: $h = 1$.

В результате получим уравнения для $n$-й гармоники Лагерра ${{E}_{{ny}}}$, ${{H}_{{nx}}}$, ${{H}_{{nz}}}$ :

Здесь опущен индекс $y$ у гармоники Лагерра электрического поля ${{E}_{{ky}}}$.

Рассмотрим разностную аппроксимацию этих уравнений.

3. АППРОКСИМАЦИЯ УРАВНЕНИЙ

Определим электрическое поле ${{E}_{{ny}}}$ и параметры $\varepsilon $, $\sigma $, $\mu $, ${{J}_{{ny}}}$ в целых $i,j$ узлах разностной сетки. Магнитное поле ${{H}_{{nx}}}$, ${{H}_{{nz}}}$ определим в промежуточных узлах $i + {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0em} 2}$, $j + {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0em} 2}$ разностной сетки. Первые производные заменим конечными разностями второго порядка аппроксимации (см. [3, 4]).

Запишем уравнения (4) в разностном виде, используя средние значения гармоник Лагерра:

В первом уравнении у компонент магнитного поля опущен номер гармоники. В правой части этого уравнения гармоники поля ${{E}_{k}}$ заменены средними по пяти точкам (см. [13–15]):

где $c$, $d$, $g$ – весовые множители.

В [13] рассмотрен случай, когда накладывалось ограничение на коэффициенты

После замены $с = (1 - 2d - 2g)$ в (6) остаются два параметра.

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

В разностных производных по $z$ использованы средние значения (см. [12]) для поля вида

и в разностных производных по $x$ использованы средние значения вида

Разностные уравнения (5), содержащие дополнительные параметры, аппроксимируют уравнения (4) со вторым порядком точности.

Подберем введенные параметры $\alpha $, $\beta $, $c$, $d$, $g$ таким образом, чтобы такая аппроксимация была оптимальной для точности решения уравнений.

4. ВЫБОР ОПТИМАЛЬНЫХ ПАРАМЕТРОВ

В однородной среде уравнения (4) для нулевой гармоники поля $E = {{E}_{{0y}}}$ можно представить в виде уравнения, которое без учета источников принимают простой вид

Это уравнение на разностной сетке можно записать, применяя средние значения поля, приведенные в предыдущем разделе:

Уравнение (10) имеет точное решение

Подставим это решение в разностное уравнение (11). После простых преобразований получим уравнение вида ${{V}^{2}}(\theta ,k) = 1$, где

и $r = {{\Delta x} \mathord{\left/ {\vphantom {{\Delta x} {\Delta z}}} \right. \kern-0em} {\Delta z}}$, $k = {{k}_{0}}\Delta x$

Будем искать параметры $\alpha $, $\beta $, $c$, $d$, $g$, требуя максимально точного выполнения равенства $V(\theta ,k) = 1$ в пределах допустимых значений $\theta ,\,k$. Для этого будем минимизировать по параметрам $\alpha $, $\beta $, $c$, $d$, $g$ функционал

Пределы интегрирования по углу $\theta = [0,{\pi \mathord{\left/ {\vphantom {\pi 2}} \right. \kern-0em} 2}]$. Пределы интегрирования по второй переменной от $k = 0$ до $k = K$. Величина $k$ определяет отношение шага разностной сетки $\Delta x$ к характерному размеру $1{\text{/}}{{k}_{0}}$ изменения решения. Поэтому брать величину верхнего предела интегрирования $K$ значительно больше единицы не имеет смысла по причине очевидной потери точности.

Если в (13) под интегралом вместо ${{(1 - V(\theta ,k))}^{2}}$ использовать ${{\left( {1 - {{V}^{2}}(\theta ,k)} \right)}^{2}}$, то оптимальные значения параметров в этих двух случаях совпадают с точностью до 4–5 знаков после запятой.

Из вида разностной схемы (5) следует, что использование оптимальных параметров не повышает порядок аппроксимации уравнений. При выполнении условия (7) это – схема 2 -го порядка аппроксимации по шагу разностной схемы. При невыполнении условия (7) нельзя указать порядок аппроксимации по шагу разностной схемы, но точность аппроксимации дифференциальных уравнений этой разностной схемой выше, чем точность аппроксимации неоптимальной разностной схемой. Таким образом, обе схемы аппроксимируют дифференциальные уравнения независимо от условия (7). Порядок аппроксимации по шагу разностной схемы и точность аппроксимации разностной схемы – это разные величины.

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

Это следует из сравнения значений величины ${{(1 - V(\theta ,k))}^{2}}$ при оптимальных значениях параметров и при неоптимальных значениях параметров, т.е. при $\alpha = 1$, $\beta = 1$, $c = 1$, $d = 0$, $g = 0$ в области допустимых $\theta ,k$. Ввиду простоты вида $V(\theta ,k)$это легко сделать. Это и подтверждают приводимые ниже решения уравнений с использованием оптимальных разностных схем.

Разностная схема без оптимальных параметров является 5-точечной, а разностная схема с оптимальными параметрами является 9-точечной. Увеличение числа точек ведет к увеличению числа вычислительных операций.

Если положить $\alpha = 1$, $\beta = 1$ и искать минимум функционала по трем параметрам $c$, $d$, $g$, то такая оптимальная схема остается 5-точечной, как и обычная схема 2 -го порядка. Такая простая модернизация обычной схемы 2 -го порядка не меняет структуры матрицы системы разностных уравнений (5).

Опустив производные по $x$ в системе (1), получим одномерные уравнения Максвелла. Тогда в системе (5) нет усреднения по пространству (8)–(9) и остаются только два оптимальных параметра $c$, $d$. Положив в (12) $\alpha = 1$, $\beta = 1$, $\theta = 0$, $\Delta x = \Delta z$, $g = 0$ и в функционале (13) оставив интегрирование только по $k$, получим формулу для ${{V}^{2}}(\theta ,k)$ в одномерном случае.

Производные функционала (13) по параметрам $\alpha $, $\beta $, $c$, $d$, $g$ вычисляются аналитически. Поэтому используя известные методы минимизации легко получить оптимальные значения $\alpha $, $\beta $, $c$, $d$, $g$. В данном случае для этого использовался метод Ньютона. Отметим, что на всех стадиях расчетов параметров $\alpha $, $\beta $, $c$, $d$, $g$ функция ${{V}^{{}}}(\theta ,k)$ (как и ${{V}^{2}}(\theta ,k)$) остается положительной и близкой к единице.

В табл. 1 приведены некоторые оптимальные значения параметров и функционала $I(r,K)$ в зависимости от соотношения шагов разностной сетки $r = \Delta x{\text{/}}\Delta z$ и верхнего предела $K$ в функционале (13) для трех разных случаев .

Первые четыре строчки таблицы получены при минимизации по пяти параметрам $\alpha $, $\beta $, $c$, $d$, $g$ и при $r = {{\Delta x} \mathord{\left/ {\vphantom {{\Delta x} {\Delta z}}} \right. \kern-0em} {\Delta z}} = 1.5$.

Следующие четыре строчки получены при минимизации по трем параметрам $c$, $d$, $g$ также при $r = \Delta x{\text{/}}\Delta z = 1.5$, но в отсутствие усреднения в пространственных производных, т.е. при $\alpha = 1$, $\beta = 1$.

Последние четыре строчки таблицы соответствуют одномерному случаю. Здесь минимизация проводилась по двум параметрам $\,c,\,d$.

Малые значения $K < 0.5$ соответствуют слабо меняющимся решениям. В табл. 1 приведены значения параметров при $0.5 < K < 2.$ В этой области оптимальное значение функционала $I(r,K)$ растет с ростом $K$ в первом и третьем случаях, во втором случае значение функционала $I(r,K)$ падает с ростом $K$.

Наименьшее значение $I(r,K)$ достигается в первом случае при минимизации по пяти параметрам. Несколько больше значения $I(r,K)$ в третьем случае, который относится к одномерным уравнениям. Во втором случае, при минимизации по трем параметрам $\,c,\,d,\,g$, величина функционала $I(r,K)$ более чем на порядок превышает соответствующие величины других двух случаев.

Из значений $I(r,K)$ следует, что применение оптимальных параметров наиболее эффективно в случае оптимизации по пяти параметрам $\alpha $, $\beta $, $c$, $d$, $g$, и наименее эффективно во втором случае при минимизации по трем параметрам. Одномерный случай занимает промежуточное положение между этими двумя случаями.

Если искать точку минимума по двум параметрам $\alpha $, $\beta $ при $c = 1$, $d = 0$, $g = 0$, т.е. учитывая только усреднение пространственных производных, то независимо от значений $r = \Delta x{\text{/}}\Delta z$ и $K$ минимальное значение функционала (13) находится вблизи одного значения $I(r,K) \approx 0.6$. Таким образом, в этом случае оптимальных значений параметров $\alpha $, $\beta $, не существует.

Если в системе уравнений (5) провести обратное преобразование Лагерра, то с учетом оптимальных параметров уравнения системы примут вид

Это система дифференциально-разностных уравнений. Они разностные по пространству и дифференциальные по времени.

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

Из табл. 1 видно, что параметры $d$, $g$ значительно меньше параметра c. Для первого случая (первые четыре строчки табл. 1) разница в 8–9 раз. Для двух последних случаев разница в значениях превышает 10 раз.

Поэтому в первом уравнении системы (14) производные по времени с параметрами $d$, $g$ являются малой поправкой. Учитывая это, можно использовать итерационный алгоритм численного решения системы дифференциально-разностных уравнений Максвелла (14).

Возможны (см. [15]) различные варианты формулы (6) усреднения $\left\langle {{{E}_{k}}} \right\rangle $. Таких вариантов не менее десятка. В общем случае среднее $\left\langle {{{E}_{k}}} \right\rangle $ можно представить в виде

Если параметры ${{s}_{m}}$ значительно меньше параметра $c$, то при итерационном методе учета малых поправок в сумме (15) можно использовать большое количество пространственных точек и большее число дополнительных параметров ${{s}_{m}}$.

В работе применялся метод простой итерации, где на первом шаге решалась система уравнений при $d = 0$, $g = 0$. На следующих шагах полученные значения электрического поля ${{E}_{{i,j}}}$ использовались для аппроксимации производных по времени при параметрах $\,d,\,g$.

Такая простая схема учета оптимальных параметров делает простой программную реализацию модифицированной разностной схемы (14) . Для этого в программу с неоптимальной схемой достаточно добавить всего около двух десятков строк (операторов), что составляет менее 5% размера всей программы. Кроме того, итерации не увеличивают значительно время счета задачи.

5. РЕЗУЛЬТАТЫ ЧИСЛЕННЫХ РАСЧЕТОВ

Решение системы уравнений Максвелла (1), полученное с использованием разностной схемы 2 -го порядка аппроксимации с оптимальными параметрами (14), будем сравнивать с решением этой системы уравнений высокоточным (см. [18]) конечно-разностным методом (4-го порядка аппроксимации по времени и 4-го порядка аппроксимации по пространству).

Источник тока брался в виде

где ${{f}_{0}}$ – несущая частота источника, ${{t}_{0}}$ – момент центра импульса источника, ${{z}_{s}}$ – точка расположения источника.

Точность решения оценивалась по величине относительной погрешности решения $D$, которая определялась выражением

где ${{E}_{y}} = E$ – решение уравнения (14), полученное с использованием оптимальных параметров, ${{E}_{{0y}}}$ – решение уравнений (1), полученное с использованием конечно-разностной схемы 4 -го порядка аппроксимации. На графиках, приведенных ниже, решение этой конечно-разностной схемой показано сплошной линией.

На фиг. 1 показано прохождение электромагнитной волны от точечного источника (${{E}_{y}}$ – компоненты поля) через слой, расположенный в однородной среде. Точками показано расположение слоя в среде. Показанные решения получены после пробега импульсом расстояния примерно в 12–15 длин волн.

Сплошная линия соответствует решению уравнений (1) разностной схемой 4 -го порядка аппроксимации.

Штрихованная линия соответствует решению обычной, неоптимальной разностной схемой 2 -го порядка. Указана величина погрешности $D1$ для этого решения, вычисленная по формуле (16).

На этом же рисунке приведено решение оптимальной разностной схемой 2 -го порядка аппроксимации (14) с пятью оптимальными параметрами $\alpha $, $\beta $, $c$, $d$, $g$, полученными при $\Delta x = 0.05$, $r = 1.5$, $K = 2.0$ и приведенными в табл. 1.

Это решение графически совпадает с решением, полученным разностной схемой 4 -го порядка аппроксимации. Для этого решения также приведена величина погрешности $D2$, вычисленная по формуле (16).

Видно, графически и численно, что оптимальная разностная схема с пятью оптимальными параметрами дает более точное решение. Погрешность $D2$ примерно на порядок меньше погрешности $D1$ неоптимальной схемы 2 -го порядка.

Для сравнения на этом же рисунке приведена погрешность решения $D3$ разностной схемой 2 ‑го порядка аппроксимации (14) с тремя оптимальными параметрами $c$, $d$, $g$, полученными при $\alpha = 1$, $\beta = 1$, $\Delta x = 0.05$, $r = 1.5$, $K = 2.0$ и приведенными в табл. 1.

Погрешность $D3$ в 2–3 раза меньше погрешности $D1$ неоптимальной схемы 2 -го порядка.

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

На фиг. 2 приведено решение одномерных уравнений Максвелла. Показано прохождение электромагнитной волны от точечного источника (${{E}_{y}}$ – компоненты поля) через слой, расположенный в однородной среде. Точками показано расположение слоя в среде. Показанные решения получены после пробега импульсом расстояния примерно в 10–12 длин волн.

В разностной схеме (14) в одномерном случае опущены разностные производные по $x$, параметры $g = 0$, $\alpha = 1$, $\beta = 1$.

Здесь сплошная линия соответствует решению одномерных уравнений Максвелла разностной схемой 4 -го порядка аппроксимации.

Штрихованная линия соответствует решению обычной, неоптимальной разностной схемой 2 -го порядка. Указана величина погрешности $D1$ для этого решения, вычисленная по формуле (16).

На этом же рисунке приведено решение оптимальной разностной схемой 2 -го порядка аппроксимации (14) с двумя оптимальными параметрами $c,\,d$, полученными при $\Delta z = 0.035$, $K = 1.5$ и приведенными в табл. 1.

Это решение графически совпадает с решением, полученным разностной схемой 4 -го порядка аппроксимации. Для этого решения также приведена величина погрешности $D2$, вычисленная по формуле (16).

Оптимальная разностная схема с двумя оптимальными параметрами дает более точное решение. Погрешность $D2$ более чем на порядок меньше погрешности $D1$ неоптимальной схемы 2‑го порядка.

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

Применение разностных схем с оптимальными параметрами позволяет повысить точность решения уравнений Максвелла в сравнении с обычными схемами 2 -го порядка аппроксимации. Это верно для одномерных и двумерных уравнений. Использование преобразования Лагерра позволяет получить оптимальные параметры и построить дифференциально-разностные уравнения. Эти уравнения можно решать, используя метод итерации по малым оптимальным параметрам.