Численное решение уравнений смешанного типа в неограниченной области на плоскости Текст научной статьи по специальности «Математика»

УДК 519.6

БОТ: 10.18698/2309-3684-2023-3-105124

Численное решение уравнений смешанного типа в неограниченной области на плоскости

© М.П. Галанин1, АР. Ухова2

1ИПМ им. М.В. Келдыша РАН, Москва, 125047, Россия 2МГТУ им. Н.Э. Баумана, Москва, 105005, Россия

Целью является построение и реализация алгоритма нахождения численного решения задачи для уравнений смешанного типа в неограниченной области. Рассматриваются задачи, в которых исследуемый процесс описывается в некоторой ограниченной области уравнением теплопроводности или волновым, а вне нее — уравнением Лапласа. Поставлены необходимые дополнительные условия в нуле, на бесконечности и условия сопряжения на границе внутренней области. Описан алгоритм нахождения численного решения задачи с волновым уравнением в ограниченной области в одномерном и двумерном случаях, задач с уравнением теплопроводности или волновым в двумерном случае. Разностные схемы построены интегро-интер-поляционным методом. Задача решается в ограниченной области. На ее границе поставлены нелокальные граничные условия так, что решение поставленной задачи в ограниченной области совпадает с проекцией на нее решения задачи в неограниченной области. При этом для решения введена искусственная граница в части области, в которой процесс описывается уравнением Лапласа. Построены итерационный алгоритм и алгоритм с нелокальным граничным условием. Представлены результаты вычислений для примеров в различных областях.

Ключевые слова: уравнения смешанного типа, неограниченная область, итерационный алгоритм, конечные разности

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

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

Решаемая задача для случая уравнения теплопроводности во внутренней ограниченной области имеет вид:

и( = А и, г е £); и\1=0=а(г), геД Аи = 0, ге№"\Д [и] = 0, ди

(1)

дп _ | и |< +ГО.

= 0, на дВ

Здесь знак квадратных скобок означает скачок величины. Необходимо разработать алгоритм решения задачи в ограниченной области Б2 так, чтобы решение задачи (1) и решение задачи в

ограниченной области Б2 совпадали в Вг.

Ниже рассмотрены следующие случаи: внутри области Б волновое уравнение при п = 1, 2 или уравнение теплопроводности при п = 2.

Иллюстрация задачи для ограниченной области В2: Б е В2 в виде

круга радиуса Я2 приведена на рис. 1. Именно такие варианты будут

рассмотрены в работе (в одномерном случае область — отрезок).

Рис. 1. Рассматриваемая область

Для численного решения задач в неограниченной области разработано множество методов, например, метод замены переменных, метод граничных интегральных уравнений, метод разностных потенциалов [1-4], метод введения бесконечных элементов совместно с конечными элементами [5, 6], использование квазиравномерных сеток [7], однако данные методы не всегда позволяют построить эффективный вычислительный алгоритм.

Для решения задачи во введенной области Б2 на ее границе необходимо поставить искусственные граничные условия. Различные способы задания таких условий описаны в [3, 4, 8, 9]. Искусственные граничные условия также можно задать с помощью основной интегральной формулы Грина. В [10] данный метод используется при решении внешней задачи для уравнения Лапласа в двумерном и трехмерном случаях, при этом для учета граничного условия также построен и исследован итерационный алгоритм. В [11, 12] алгоритм применен для решения задач смешанного типа с уравнением теплопроводности в одномерной ограниченной области и с волновым уравнением в двумерной ограниченной области соответственно, построен итерационный алгоритм для уточнения граничного условия.

Для этого в [10] рассмотрена дополнительная граница дД в области /)2 \ И, представляющая собой окружность радиуса Яг. Решение

внешней задачи для уравнения Лапласа для круга известно и задается интегралом Пуассона [13-15]. В двумерном случае он имеет вид:

1 2г г2—Я2

и(г,т) = — ( и(Я,ш) —--ёш, г > Я.

2я{ г2 - 2гЯ со8(^-ш) + Я

Вычислив интеграл Пуассона от значений искомой функции на дД, можно получить значения функции на границе области Д . Тогда исходная задача в части уравнения Лапласа примет вид:

Аи = 0, геД\Д

1 2 г Я2 - Я 2

и \дп =— ( и(Я,ш)—0---1-7ёш;

1д°2 2г( Я - 2ЗДсо8(р-ш) + Я (2)

[и] = 0;

ди

дп

= 0 на дБ.

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

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

Одномерное волновое уравнение. Плоский случай. В данном случае область Д — отрезок [0, Я ], Д — отрезок [0, Я2 ], а условие (2) в случае декартовых и цилиндрических координат становится равенством значений искомой функции в точках Я и Я . Рассмотрим задачу в декартовых координатах:

<

и„ = и„

0 < х < Я, г > 0;

и |,=0 = а(х), 0 < х < Я; щ |(=0 = Р(х), 0 < х < Я; и |х=0 = И>);

их = 0 Я < х < Я2;

и |х=Я2 = и |х=Я1 ; и |х=Я-0 = и |х=Я+0; их |х=Я-0 = их |х=Я+0 •

Построим разностную схему для решения задачи. Введем равномерную сетку в области 01 = {0 < х < Я ,0 < г < Т} с шагами Нг по пространству и Н по времени. Сетку построим так, что узел с номером N попадает в точку х = Я, N — в точку х = Я , а N — в точку х = Я2.

Схема с весами для волнового уравнения имеет вид:

(3)

У = У—

■'и ^ хх

Схема для уравнения Лапласа аналогична. Здесь

у(а) = а у + (1 - 2а) у + а у,

(4)

где у — сеточная функция, приближающая искомое решение и ; а — вес. Обозначения разностных производных и сеточных функций на разных временных слоях соответствуют [16, 17].

Данная схема является трехслойной, для нахождения значения ук+1 необходимо знать решения с двух предыдущих слоев. Значение на нулевом слое известно из начального условия и |г=0 =а(х), на первом слое можно вычислить решение из условия щ |г=0 = /3(х). Схема (3) имеет второй порядок аппроксимации по пространству и времени, поэтому аппроксимируем второе начальное условие тоже со вторым порядком. Для этого воспользуемся формулой Тейлора:

у = у0 + Н А х) + 0,5 • Н^а-х (х )•

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

у - = 2 у(а)- •

^ N ,гг ^ N ,хх

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

Для расчетов будем использовать симметричную схему, соответствующую а = 0,5 . Она является безусловно устойчивой.

Пример 1. Возьмем в качестве условий задачи

а( х) = соб(юс),

А х) = о,

щ(г) = соб(ж1 ). Точным решением является:

\соъ(яг) соб(^х), 0 < х < Я;

и(х, г) = ■

соб(^)соб(^Я), х > Я.

Размеры областей Я = 1, Я = 1,5, Я = 2 . Начальный шаг возьмем равным к = (кг; к) = (0,02; 0,0196078)

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

Таблица 1

Ошибки тестового решения для различных шагов

Шаг к к 2 к 4 к 8

Ошибка 1,67 •Ю-3 4,26 •Ю-4 1,08 •Ю-4 2,70 •Ю-5

Отношение — 3,92 3,96 3,98

Сравнение решений при уменьшении к показывает, что схема имеет второй порядок сходимости. Рис. 2 показывает решение при г е [0,1].

и 1,0

0,5

0,0

-0,5

.0,5-..

1,0

1,5

2,0 х

-1,0

Рис 2. Решение примера 1 в различные моменты времени:

— — г = 0;----г = 0,4;

----г = 0,6; • ---г = 1,0

Одномерное волновое уравнение. Цилиндрический случай. В

цилиндрических координатах задача в В2 принимает вид:

1 а

с

и„ =--

г аг

аи

г — V аг

и \г=0 = а(г), 0 <р< Я;

0 < г < Я, г > 0;

у

и

г ч=0

= Дг), 0 <р< Я;

1 а

г аг и

аи аг

= 0, Я < г < Я;

(5)

г=я и |г=Я, ;

иг |г=Я-0 иг |г=Я+0 •

Воспользуемся равномерной сеткой, построенной в плоском случае. Для данной задачи схема с весами для волнового уравнения имеет вид:

у гК (

г у(а) - г у(а) ) '+0,5 уг '-0,5 уг )•

Значение функции на первом слое найдем по формуле Тейлора:

у1 = у0 + К р( г) + 0,5 • К2^1- (г+0,5аг - г - 0,5а) •

Аппроксимируем условие непрерывности производной функции в точке г = Я так же, как и в плоском случае. Для этого проинтегрируем волновое уравнение на отрезке [Я - Нг /2, Я], а уравнение Лапласа —

на отрезке [Я, Я + Нг /2]. После преобразований получим:

2 ( _ (

у~и = ^"(г+ 0,5у( ' - г-0,5уг

,(а)

(а)

Далее в примерах радиус Я внутренней области Б выбран так, чтобы обеспечить выполнение решением условия непрерывности нормальной производной на границе аБ .

В нуле ставится условие ограниченности. Для аппроксимации этого условия проинтегрируем волновое уравнение (4) на отрезке [0, К /2] и после преобразований получим уравнение схемы. Оно имеет вид:

4

у 0,й 1. у0

(а)

г

Рассмотрим далее симметричную схему при а = 0,5 . Схема имеет аппроксимацию О(h2 + hr2) .

Пример 2. Решим задачу с начальными условиями

a{r) = Jо{яг ), р{г ) = nJo(^r ).

Ее точное решение:

\(cos(nt) + sin(^t)) J0 (^r), 0 < r < R; [(cos(^t) + sin(^t)) J0 (^R), r > R.

Здесь и далее J (z) — функция Бесселя к -го порядка. Возьмем R = 1,21967, R = 1,75, R = 2,43934. Начальный шаг h = (hr; ht) = (0,0243934; 0,0196078) В табл. 2 приведены ошибки решения при уменьшении шага h .

Таблица 2

Ошибки тестового решения для различных шагов

u(r, t) = •

Шаг h h 2 h 4 h 8

Ошибка 3,49 •Ю-3 8,72 •Ю-4 2,18 •Ю-4 5,48 •Ю-5

Отношение — 4,00 3,99 3,99

Сравнение решений примера при уменьшении шага к показывает, что схема имеет второй порядок сходимости. На рис. 3 представлено решение задачи при t е [0,1].

1,5

1,0

0,5

0,0

-0,5

0,5 V 1,0-------1,5------2,0-------2,5 r

-1,0

Рис 3. Решение примера 2 в различные моменты времени: — — t = 0;----t = 0,3;----t = 0,7; • ---t = 1,0

u

Двумерная задача с нелокальным граничным условием. Волновое уравнение. В случае, когда в области 1>еМ2 имеет место волновое уравнение, исходная задача выглядит следующим образом:

1 а

и" =

г аг

аи аг

1 а и

Г2 аф2"

(г, ф) е Б, г > 0;

и |г=0 =а(г,ф), (г,ф) е Б; = р(г,ф), (г, ф) е Б;

иг |г=0

1 д ( диЛ 1 д2и . . . „2 . ^ г дг\ аг

г2 аф2

е аБ; = 0, (г,ф) е аБ;

(6)

[и ] = 0, (г,ф) е аБ;

аи _ап

и <

Рассмотрим область Б2 — круг радиуса . Введем равномерную сетку в 0 = {Д / 2<г <^,0< ф <2я",0<г <Т} с шагом кг по координате г, шагом Ь по углу ф и шагом кг по времени г. При этом первая точка по г сдвинута на половину шага от нуля. Узел под номером N попадает на границу аБ, N — на аБ1 и N2 — на аБ2.

Как и ранее, построим разностную схему для волнового уравнения интегро-интерполяционным методом [16, 17]. Получим следующую разностную схему с весами для волнового уравнения (5):

гЬ

Схема является трехслойной, для нахождения значения у необходимо знать значения с двух предыдущих слоев. Значение на первом слое найдем с помощью формулы Тейлора:

У = У + ЬгР + 0,5 • к]

I

гк

1 / л 1

— г0,5«г - гг-0,5« ) + " «ф

Аналогично волновому построим схему для уравнения Лапласа:

(7)

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

k V R2 - R2)

yN+m

m-0,Np R - 2R1R COS(^; - ) + R

k+1

Nv-1 2 V

^ s Nm

-1 R2 -2RiR2cos(vJ -Pm) + Ri2

Пусть область D — круг радиуса R . Решение задачи должно быть ограничено в нуль по радиусу. Проинтегрируем первое уравнение (5)

по ячейке [0, r05 ] х [р_0 5, р+0 5 ] и получим следующую схему:

V - — V(ff) + — V(а)

Jtt и и2 S PP

hr hr

На границе области D поставлены условия сопряжения. Для аппроксимации условия равенства потоков рассмотрим ячейки

[rN-0,5 , rN Х [Pj-0,5 , Pj +0,5 ] и [rN > rN+0,5 ] Х [Pj-0,5' Pj+0,5 ] . В первой ячейке

проинтегрируем волновое уравнение, а во второй — уравнение Лапласа. Из данных равенств выразим поток в точке rN, приравняем и получим аппроксимацию условия сопряжения:

V = — О v(e7)-г ) + —v(a)

rh +0'5 г r2 w'

Далее использована симметричная схема с а - 0,5 . Построенный метод имеет порядок аппроксимации 0(h2 + h2 + h2).

Для решения системы линейных алгебраических уравнений схемы использовался решатель BiCGStab — стабилизированный метод бисо-пряженных градиентов для задач с разреженными матрицами из библиотеки Eigen [18].

Примеры применения представлены в разделе, посвященном итерационному алгоритму.

Двумерная задача с нелокальным граничным условием. Уравнение теплопроводности. Рассмотрим задачу (1) в данном случае. Рассмотрим область D2 такую же, как и для (5). Точно так же введем сетку в области Qj.

Для построения схемы уравнения теплопроводности воспользуемся интегро-интерполяционным методом. В итоге разностная схема с весами для уравнения теплопроводности выглядит следующим образом:

m

Аппроксимация уравнения Лапласа аналогична полученной выше. Интеграл Пуассона приближенно вычисляется по формуле трапеций. Ошибка аппроксимации схемы для уравнения теплопроводности при любом значении веса а — О(Д + кг2 + к2) [16, 17].

Рассмотрим примеры. Пусть внутренняя область Б — круг радиуса Я . Решение задачи должно быть ограничено в нуль по радиусу. Для аппроксимации данного условия необходимо проинтегрировать

уравнение теплопроводности по ячейке _0, г05 ]х[ф-_0 5ф+о5] . Получим выражение, аналогичное волновому уравнению:

2 (а) . 4 (а)

у =-уУ ' + —- уУ-'.

■Уг к к2 фф

При расчетах использована неявная схема, соответствующая весу а = 1. Система линейных уравнений схемы решена с использованием решателя BiCGStab [18].

Пример 3. Начальное условие:

а(г, ф) = (г) зт(ф).

Точное решение задачи таково:

ехр(-г)(г) в1п(ф), 0 < г < Я;

и(г,ф, 0) =

, ч 1,24846 . . . ехр(-г)-Б1п(ф), г > Я.

Возьмем Я = 2,40483 , Я = 2Я , Я2 = 4Я .

Начальные шаги сетки к = 0,1 и к = (Нг; \) = (0,229031; 0,785538). В табл. 3 приведены ошибки решения для разных шагов к и к .

Таблица 3

Ошибки тестового решения для различных шагов

<

г

Шаг к, к к 2, Ы 4 к/ 4, 16 к/ 8, 64

Ошибка 0,018200 0,004600 0,001160 0,000280

Отношение — 3,85 3,87 3,88

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

Итерационный алгоритм для интегрального граничного условия. Построим итерационный алгоритм для учета интегрального граничного условия (2) [12]. Значение и* на границе области Б2 на каждой новой итерации будет вычисляться через уже известные значения

и5 1 на предыдущей итерации, используемые в интеграле Пуассона. Получим модифицированную задачу (2) в области £)2 \ £).

Для нахождения решения итерационным способом на каждом временном слое воспользуемся схемами, построенными ранее.

Волновое уравнение. Пример 4. Решим задачу с использованием итерационного алгоритма. Возьмем начальные условия

а(г ,ф) = 3о(г ); Р(г,ф) = 0.

Точным решением задачи является функция

Гсо8(03 (г), 0 < г < Я;

ы(г,ф, г) = •

со8(Г)3 (Я), г > Я.

Возьмем Я = 3,83171, Я = 2Я, Я2 = 4Я .

Начальные шаги сетки к = (кг; \; \) = (0,364924; 0,785398; 0,1). В табл. 4 приведены ошибки решения при уменьшении шага к .

Таблица 4

Ошибки тестового решения для различных шагов

Шаг к к 2 к 4 к 8

Ошибка 0,015100 0,003190 0,000838 0,000214

Отношение — 4,74 3,81 3,91

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

Исследуем влияние положения границ на сходимость. Для этого зафиксируем Я = 3,83171 и Я = 8Я, а радиус Я1 будем менять.

В табл. 5 приведено количество итераций, необходимое для получения решения при различных значениях радиуса Я .

Таблица 5

Количество итераций при различных положениях границы 5Д

Я\1 Я 1,5 2,0 2,5 3,0 3,5 4,0 4,5 5,0

Количество 9 12 15 18 22 27 33 45

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

значений решения с текущей и предыдущей итераций меньше 10-5). Полученные результаты подтверждают оценки [12].

Пример 5. Рассмотрим задачу с начальными условиями

а(г,ф) = 0; Р(г ,ф) = Зх(г )81п(ф). Точное решение задачи таково:

81п(г)^ (г) 81п(ф), 0 < г < Я;

и (г ,ф, 0) =

. , Л 1,24846 . . . 81п(г)-81п(ф), г > Я.

Исследуем количество итераций в зависимости от границы Б . Зафиксируем Я = 2,40483 и Я2 = 8Я, а радиус Я1 будем менять.

В табл. 6 приведено количество итераций, необходимое для получения решения при различных значениях радиуса Я .

Таблица 6

Количество итераций при различных положениях границы 5Д

Я,/ Я 1,5 2,0 2,5 3,0 3,5 4,0 4,5 5,0 5,5 6,0

Количество 5 5 6 7 8 9 10 12 16 31

Из табл. 6 (аналогично табл. 5 примера 4) видно, что чем ближе дополнительная граница к границе внутренней области, тем меньше итераций требуется для получения решения с заданной. Результаты также подтверждают теоретические оценки [12].

Пример 6. Рассмотрим задачу, в которой внутренняя область Б имеет форму полукруга радиуса Я (рис. 4).

Условия сопряжения при г = Я, -л/ 2 < ф < л/ 2 не меняются. На левой границе области Б путем интегрирования уравнений по ячейке

-0,5, г+0,5 ]х

ф

ф

МБ-0,5'^ МБ +0,5

ТБ

где Мф — номер узла по ф, попадающего на левую границу Б , получим аппроксимацию условия сопряжения:

У =

гк

[ г+0,у(а) - г-0

5 У

(а)

) + 4 у (а).

! гг2 ^ фф

В качестве начальных условий возьмем:

а (г ,ф) = 0; /3(г,ф) = г 81п(ф).

г

При расчетах Я = 1, Я = 2Я, Я2 = 4Я.

Начальный шаг сетки к = 0,05 и к = (кг;к ) = (0,0852381; 0,785388).

В данном случае точное решение не известно. Поэтому для проверки алгоритма в качестве тестового решения возьмем решение данной задачи без использования итераций со следующими изменениями. Вместо условия (2) на внешней границе поставим условие равенства нулю производной по радиусу, которое часто используется для имитации условия на бесконечности. Радиус внешней области Я увеличим, чтобы уменьшить влияние данного граничного условия на решение в исходной области Б2. Данный метод описан, например, в [19]. Таким

образом, необходимо решить задачу в области Б2' радиуса Я' = дЯ, д > 1 — параметр, Я — радиус исходной области Б2, с условием равенства нулю производной вместо нелокального условия. Его аппроксимируем также интегро-интерполяционным методом.

Решим данную задачу при различных значениях д , чтобы посмотреть, как размер области Б2' влияет на решение в исходной области Б2. В табл. 7 представлены результаты вычисления относительной ошибки совпадения решения задачи в исходной области с решением задачи в большей области при различных положениях внешней границы области Б2'. Вид решения представлен на рис. 5.

Для вычисления ошибки использована равномерная норма. Относительная ошибка решений мала и уменьшается с увеличением д, Это означает, что проекции решения данной задачи при различных значениях д на область исходной задачи близки по значению.

Рис. 4. Область — полукруг для примера 6

Таблица 7

Сравнение решений вспомогательной задачи при различных значениях q в области Д решения исходной задачи для примера 6

ч 7 9 11 13

Относительная ошибка 0,004820 0,002020 0,001040 0,000702

Теперь сравним решения исходной задачи с использованием итерационного алгоритма для уточнения решения на внешней границе с

решением тестовой задачи при различных размерах области Д'. В табл. 8 приведены результаты сравнения.

Таблица 8

Сравнение решений вспомогательной задачи при различных значениях q в области Д решения исходной задачи для примера 6

ч 5 7 11 13

Ошибка 2,27 -10-3 1,17-10-3 9,00-10-4 7,18-10-4

Из таблицы видно, что чем больше область решения тестовой задачи, тем меньше ошибка при сравнении с решением исходной задачи.

-4

0,2

0,0

-0,2

4

Рис. 5. Решение примера 6 при г = 1

Уравнение теплопроводности. Пример 7. Рассмотрим задачу с начальными условиями:

а(г,р) = 30(г); Р(г,р) =

Точным решением задачи является функция

и(г, р, г) =

[ехр(-г)3 (г), 0 < г < Я; I ехр(-г)3 (Я), г > Я.

4

Возьмем Я = 3,83171, Я = 2Я, Я = 4Я.

Начальные шаги сетки к = 0,1 и к = (кг; к ) = (0,364824; 0,785388).

В табл. 9 приведены ошибки решения для разных шагов к и к . Само решение представлено на рис. 6.

Таблица 9

Ошибки тестового решения для различных шагов

Шаг к, к и/ 2, Ы 4 И/4, к/16 к/8, к/64

Ошибка 0,022000 0,005780 0,001470 0,000368

Отношение — 3,81 3,84 3,88

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

Исследуем зависимость числа итераций от положения границы 5Б1. Будем менять положение 5Б1, а остальные границы зафиксируем. В табл. 10 приведено количество итераций, необходимое для получения решения при различных значениях радиуса Я .

Таблица 10

Количество итераций при различных положениях границы д01

Я/ Я 1,5 2,0 2,5 3,0 3,5 4,0 4,5 5,0 5,5

Количество 11 14 17 21 25 30 37 51 99

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

Полученные результаты подтверждают теоретические оценки [12]. Однако количество итераций в случае табл. 10 заметно больше таковых из табл. 5, 6. Объясняется это различием решаемых уравнений (теплопроводности и волнового соответственно).

Пример 8. Возьмем начальное условие и точное решение из примера 3. Зафиксируем Я = 2,40483, Я = 2Я и Я = 4Я . Начальные шаги сетки к = 0,1 и к = (к;к) = (0,229031; 0,785398).

В табл. 11 приведены ошибки решения для разных шагов к и к .

Таблица 11

Ошибки тестового решения для различных шагов

Шаг к, к к/ 2, Ы 4 к/4, к,/16 к/ 8, 64

Ошибка 0,016800 0,004280 0,001080 0,000271

Отношение — 3,92 3,96 3,99

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

Приведенные данные также показывают, что ошибка при решении с использованием итерационного алгоритма меньше, чем в безытерационном случае.

Проведем исследование влияния положения дополнительной границы на количество итераций. Будем менять радиус Я.

В табл. 12 приведено количество итераций, необходимое для получения решения при различных значениях радиуса Я .

Таблица 12

Количество итераций при различных положениях границы дД

Я,/ Я 1,5 2,0 2,5 3,0 3,5 4,0 4,5 5,0 5,5 6,0

Количество 5 6 6 7 8 9 11 13 17 25

Чем ближе дополнительная граница к границе внутренней области, тем меньше итераций требуется для получения решения.

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

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

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

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

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

Рассмотрен случай, когда внутренняя область D — полукруг радиуса R . Численные решения задачи в такой области при использовании итерационного алгоритма и с нелокальным граничным условием близки.

Исследование выполнено за счёт Российского научного фонда (грант 22-21-00260).

ЛИТЕРАТУРА

[1] Koleva M.N. Numerical solution of the heat equation in unbounded domains using quasi-uniform grids. Lecture Notes in Computer Science, 2006, vol. 3743, pp. 509-517.

[2] Рябенький В.С. Метод разностных потенциалов для некоторых задач механики сплошных сред. Москва, Наука, 1987, 391 с.

[3] Брушлинский К.В., Рябенький В.С., Тузова Н.Б. Перенос граничного условия через вакуум в осесимметричных задачах. Журнал вычислительной математики и математической физики, 1992, т. 32, № 12, с. 1929-1939.

[4] Брушлинский К.В. Математические и вычислительные задачи магнитной гидродинамики. Москва, БИНОМ. Лаборатория знаний, 2009, 200 с.

[5] Bettess P. Infinite Elements. Paris, Penshaw Press., 1992, 264 p.

[6] Zienkiewicz O.C., Emson C., Bettess P. A novel boundary infinite element. International Journal for Numerical Methods in Engineering, 1983, vol. 83, no. 3, pp. 393-404.

[7] Калиткин Н.Н., Алъшин А.Б., Алъшина Е.А., Рогов Б.В. Вычисления на квазиравномерных сетках. Москва, Физматлит, 2005, 223 с.

[8] Галанин М.П., Низкая Т.В. Разработка и применение численного метода линейных эллиптических уравнений в неограниченной области. Препринты ИПМим. М.В. Келдыша, 2005, № 2, с. 1-29.

[9] Tsynkov S.V. Numerical solution of problems on unbounded domains. A review. Applied Numerical Mathematics, 1998, vol. 27, iss. 4, pp. 465-532.

[10] Галанин М.П., Сорокин Д.Л. О решении внешних краевых задач для уравнения Лапласа. Дифференциальные уравнения, 2020, т. 56, № 7, с. 918-926.

[11] Галанин М.П., Сорокин Д.Л., Ухова А.Р. Методы численного решения дифференциального уравнения смешанного типа в неограниченной области. Математическое моделирование и численные методы, 2021, № 1, с. 91-109.

[12] Галанин М.П., Сорокин Д.Л., Ухова А.Р. О решении уравнения смешанного типа в неограниченной области. Дифференциальные уравнения, 2022, т. 58, № 7, с. 921-929.

[13] Тихонов А.Н., Самарский А.А. Уравнения математической физики. Москва, Наука, 1972, 735 с.

[14] Свешников А.Г., Боголюбов А.Н., Кравцов В.В. Лекции по математической физике. Москва, Изд-во МГУ, 1993, 352 с.

[15] Мартинсон Л.К., Малов Ю.И. Дифференциальные уравнения математической физики. Москва, Изд-во МГТУ им. Н.Э. Баумана, 1996, 228 с.

[16] Галанин М.П., Савенков Е.Б. Методы численного анализа математических моделей. Москва, Изд-во МГТУ им. Н.Э. Баумана, 2010, 591 с.

[17] Самарский А.А. Введение в теорию разностных схем. Москва, Наука, 1971, 552 с.

[18] Eigen is a C+ + template library for linear algebra: matrices, vectors, numerical solvers, and related algorithms [Электронный ресурс]. URL: https://eigen.tuxfamily.org (дата обращения: 31.03.2022)

[19] Галанин М.П., Сорокин Д.Л. Разработка и применение численных методов решения задач в неограниченной области на основе третьей формулы Грина. Препринты ИПМим. М.В. Келдыша, 2018, № 246, с. 1-24.

Статья поступила в редакцию 27.05.2023

Ссылку на эту статью просим оформлять следующим образом: Галанин М.П., Ухова А.Р. Численное решение уравнений смешанного типа в неограниченной области на плоскости. Математическое моделирование и численные методы, 2023, № 3, с. 105-124.

Галанин Михаил Павлович — д-р физ.-мат. наук, главный научный сотрудник, и.о. заведующего отделом ИПМ им. М.В. Келдыша РАН, профессор кафедры «Прикладная математика МГТУ им. Н.Э. Баумана. e-mail: [email protected]

Ухова Анна Романовна — студентка кафедры «Прикладная математика» МГТУ им. Н.Э. Баумана. e-mail: [email protected]

Numerical solution of equations of mixed type in unlimited region on a plane

© M P. Galanin1, A.R. Ukhova2

1Keldysh Institute of Applied Mathematics, Moscow, 125047, Russia 2Bauman Moscow State Technical University, Moscow, 105005, Russia

The purpose of the work is to build and implement an algorithm for finding a numerical solution to a problem for mixed-type equations in an unlimited region. In this case problems are considered in which the process under study is described in some limited area by the thermal conductivity equation or wave equation, and outside it by the Laplace equation. The necessary additional conditions at zero, at infinity and the conditions for conjunction at the border of the inner region are set. There is described an algorithm for finding a numerical solution to a problem with a wave equation in a limited region in one-dimensional and two-dimensional cases, problems with a thermal conductivity equation or a wave equation in a two-dimensional case. Difference schemes are built by the integro-

interpolation method. The task is solved in a limited area. Nonlocal boundary conditions are set on its border so the solution of task in limited area coincides with projection of problem in unlimited area. In this case, an artificial boundary is introducedfor the solution in the part of the region in which the process is described by the Laplace equation. An iterative algorithm and an algorithm with a non-local boundary condition are built. The results of calculations for examples in various fields are presented.

Keywords: mixed type equations, unlimited domain, iterative algorithm, finite differences

REFERENCES

[1] Koleva M.N. Numerical solution of the heat equation in unbounded domains using quasi-uniform grids. Lecture Notes in Computer Science, 2006, vol. 3743, pp. 509-517.

[2] Ryabenky V.S. Metod raznostnyh potencialov dlya nekotoryh zadach mekhaniki sploshnyh sred [The method of difference potentials for some problems of continuum mechanics]. Moscow, Nauka Publ., 1987, 391 p.

[3] Brushlinskii, K.V., Ryaben'kii, V.S., Tuzova, N.B. The transfer of boundary conditions across a vacuum in axisymmetric problems. Computational Mathematics and Mathematical Physics, 1992, vol. 32, no. 12, pp. 1757-1767.

[4] Brushlinsky K.V. Matematicheskie i vychislitel'nye zadachi magnitnoj gidro-dinamiki [Mathematical and computational problems of magnetic hydrodynamics]. Moscow, BINOM. Laboratory of Knowledge Publ., 2009, 200 p.

[5] Bettess P. Infinite Elements. Paris, Penshaw Press., 1992, 264 p.

[6] Zienkiewicz O.C., Emson C., Bettess P. A novel boundary infinite element. International Journal for Numerical Methods in Engineering, 1983, vol. 83, no. 3, pp. 393-404.

[7] Kalitkin N.N., Aleshin A.B., Alshina E.A., Rogov B.V. Vychisleniya na kvaziravnomernyh setkah [Calculations on quasi-uniform grids]. Moscow, Fiz-matlit Publ., 2005, 223 p.

[8] Galanin M.P., Nizkaya T.V. Development and application of a numerical method for solution of linear elliptic equations in unbounded region. Keldysh Institute Preprints, 2005, no. 2, pp. 1-29.

[9] Tsynkov S.V. Numerical solution of problems on unbounded domains. A review. Applied Numerical Mathematics, 1998, vol. 27, iss. 4, pp. 465-532.

[10] Galanin M.P., Sorokin D.L. Solving exterior boundary value problems for the Laplace equation. Differential Equations, 2020, vol. 56, no. 7, pp. 890-899.

[11] Galanin M.P., Sorokin D.L., Ukhova A.R. Methods for numerical solution of a mixed type differential equation in an unbounded domain. Mathematical Modeling and Computational Methods, 2021, no. 1, с. 91-109.

[12] Galanin MP, Sorokin DL, Ukhova AR. On solving the equation of mixed type in an unlimited region. Differential equations, 2022, vol. 58, no. 7, pp. 921-929.

[13] Tikhonov A.N., Samarsky A.A. Uravneniya matematicheskoj fiziki [Equations of mathematical physics]. Moscow, Nauka Publ., 1972, 735 p.

[14] Sveshnikov A.G., Bogolyubov A.N., Kravtsov V.V. Lekcii po matematicheskoj fizike [Lectures on mathematical physics]. Moscow, MSU Publ., 1993, 352 p.

[15] Martinson L.K., Malov Yu.I. Differencial'nye uravneniya matematicheskoj fiziki [Differential equations of mathematical physics]. Moscow, BMSTU Publ., 1996, 228 p.

[16] Galanin M.P., Savenkov E.B. Metody chislennogo analiza matematicheskih mod-elej [Methods of numerical analysis of mathematical models]. Moscow, BMSTU Publ., 2010, 591 p.

[17] Samarskiy A.A. Vvedenie v teoriyu raznostnyh skhem [Introduction to the theory of difference schemes]. Moscow, Nauka Publ., 1971, 552 p.

[18] Eigen is a C+ + template library for linear algebra: matrices, vectors, numerical solvers, and related algorithms [Electronic resource]. URL: https://eigen.tuxfam-ily.org (accessed: 31.03.2022)

[19] Galanin M.P., Sorokin D.L. Development and application of numerical methods for solving tasks in unlimited regions based on the third green formula. Keldysh Institute Preprints, 2018, no. 246, pp. 1-24.

Galanin M.P., Dr. Sc. (Phys. — Math.), Chief Researcher, Acting Head of the Department of the Keldysh Institute of Applied Mathematics, Professor of Department of Applied Mathematics, Bauman Moscow State Technical University. e-mail: [email protected]

Ukhova A.R., student of Department of Applied Mathematics, Bauman Moscow State Technical University. e-mail: [email protected]