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

УДК 519.6

О НОВЫХ ЧИСЛЕННЫХ МЕТОДАХ РЕШЕНИЯ КРАЕВЫХ ЗАДАЧ В ОБЛАСТЯХ СЛОЖНОЙ ФОРМЫ С ТАБЛИЧНО ЗАДАННЫМИ ИСХОДНЫМИ ДАННЫМИ

С.В. Трубников

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

Ключевые слова: численные методы; краевые задачи; дифференциальные уравнения

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

Исходные данные краевых задач обычно задаются формулами, но такой способ задания функций не всегда удобен. Нередко для проведения вычислительных экспериментов требуется задавать функции с определённой формой графика. Этого можно добиться, если задавать таблицу значений функции, а затем на её основе построить гладкую аппроксимацию задаваемой функции. Если для получения гладкой аппроксимации использовать кусочно-линейную интерполяцию и частичные суммы ряда Фурье по системе пробных функций метода Галёркина, то она оказывается особенно удобной для дальнейшего применения метода Галёркина. Впервые эта идея была опробована в 2004 году в книге [3] применительно к решению задачи Дирихле для уравнения Пуассона в прямоугольной области. Затем, в 2014 году она была развита в работе [4] на случай области сложной формы. Это обобщение основано на преобразовании координат, превращающем область, в которой решается дифференциальное уравнение в прямоугольник той или иной размерности, гибридной аппроксимации исходных данных, заданных таблично (кусочно-линейной интерполяции с последующим сглаживанием путём использования частичных сумм тригонометрического ряда Фурье), проекционной схеме решения аппроксимированной краевой задачи с теми же пробными и поверочными функциями, что используются при аппроксимации исходных данных. В результате применения описанного подхода получаются проекционные вычислительные схемы для решения стационарных краевых задач математической физики.

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

Пусть G = {( x, y): 0 < x < a, y1 ( x ) < y < y2 ( x )} - обобщённая криволинейная трапеция, y - граница обла-

сти G . Здесь й - заданная положительная постоянная, у1 (х ), у2 ренцируемые нужное количество раз функции на [а, ь ] (рис. 1). Введём замену переменной

(У - У (х))

( x )

(y (x) < y2 (x)) - заданные, непрерывно диффе-

y=yj(x)+^(y2 (x)-yj(x)) «л

(у2 (x)- У1 (x))'

Отображение (1) преобразует область G в прямоугольник р = х угольника Р. В прямоугольнике Р вводится прямоугольник ^ = х ^ )

0 < x < a ,0 < л

c1 < x <

d, <

(0 < с 1 < C 2 {(x, У ): ci

< a, 0 < d1 < d

(1)

11. Г - граница прямо-< d2 } . Здесь с1 , с2 , d вводится

<

Л - -2

< 1 ) - заданные положительные постоянные. А в области G y 1 (x )+ d 1 (y 2 (x )- y 1 (x ))< y < y1 (x) + d2 (y2 (x)- y1 (x))} , которая является прообра-

1 ' ° 2 ' 1 ' область

Q

зом прямоугольника R при отображении (1) (рис. 1).

1) 2) Рис. 1.

1) области G и Q; 2) прямоугольники P и R.

с

2

d

2

Рассматривается первая краевая задача для уравнения теплопроводности в криволинейной трапеции G. Ищется дважды непрерывно дифференцируемая в области (о \ у ) х (0, Т ] и непрерывная на о х [0, Т ] функция и (х, у, () , удовле-

творяющая двумерному уравнению теплопроводности

ди _ ди( х, у,/) (д2и (х, у,/) д2и( х, у,/)

д1

DAu:

дл

D

ч дх2 дУ2 у

(х, у )е о \ у , t е (0, Т ] ,

на границе у - однородному граничному условию первого рода

и 1(х,у)еу = 0, ^ е (Т ] ,

а в начальный момент времени t = 0 - начальному условию

и ( х, у ,0 )=^( х, у ) ; ( х, у ) е О .

= / (х, у,/)

(2)

(3)

(4)

Здесь Б и Т - заданные положительные постоянные. Функция ^ (х, у, /) задаётся следующим образом. Введём

функцию F (х ,-ц, /) такую, что Р

х.

(у - у1 (х)) г '(у2 (х)- у1 (х))'

= /(х, у, t) . В параллелепипеде R х [0, т ] вводится равно-

мерная прямоугольная сетка точек |( х,tp): х1 = с1 + Нх\, г]}. = d1 + НЦу, гр = Н(р, Н =(с2 - с1)/пгх , К =(d2 - d1 )/ , Н = ТП, / = 0,1,..., Д^, у = 0,1,..., , р = 0,1,..., п{ } . Натуральные числа П^, пЦ, ^

считаются исходными данными. Координатные линии сетки делят прямоугольник R на прямоугольники Ру = |(х, Ц ) • хг-1 < х < хг, Т]._х < Г) <Г) у |, которые, в свою очередь, разбиваются диагоналями на верхний и нижний

треугольники Т и ТУ (рис. 2). Узлам введенной сетки (хг,Т).) в прямоугольнике R соответствуют узлы криволинейной

сетки (хг, угу ) = (хг , у (хг) + Т/у (у2 (хг) _ у (хг))) в области Q. В качестве основных исходных данных, определяющих функцию / (х, у, t) , задаются её значения в этих узлах / (хг, у/у., tp ) = //у.р . Эти же значения будет иметь функция F (х t) в узлах прямоугольной сетки: F (хг ,Т)., tp ) = / (хг, у/у., tp ) = //у.р . Аппроксимация функции F (х t) строится в 2 этапа. В начале строится промежуточная аппроксимация р (х, Т]) функции F (х,Т], tp ) при фиксированном Р. В качестве аппроксимации Рр (х,Г]) выбрана линейная функция в каждом из треугольников ТЦ и Т, удовлетворяющая условиям интерполяции Рр (хг,Т).) = и обращающаяся в ноль вне прямоугольника R. Построить интерполяционную функцию Рр (х,ц) не сложно. Приведём формулы для её вычисления.

1)

2)

Рис. 2.

1) сетка в прямоугольнике R (П =4, п = 3 ); 2) треугольники Т и Ту •

Рр (х,л) =

Р»р (х,г,) , (х,г,) е ГЦ, / = 1,2,...,п{с,у = 1,2,...,п^,р = 0,1,...,п/

Р*р (х,П) , (х,П) е Ту, / = 1,2,...,п1,у = 1,2,...,п^,р = 0,1,..., Д

0, (х,ц)е Р \ Я.

Здесь

(х п)=а»]Р+р;]Р (х - хм )+у;]Р (] - ] -)

^¿р=(-/-1 ] р- ^ ]-1 р V И1.

аи = /

ЧР •/г-1 1 -1 Р' аг]р = ./1-1 1 -1 р,

=(/ур - /-1 ¿р)/И

(^ 1-1 р ■/i-l у-1 р VИ/

Р1р (x, 1) = а1р + Д^ (х - х-!) + ^ (л -111 -1)

Т']р = (/г] р / 1-1 р ")/И1 .

Коэффициенты , , , а^^, , в этих формулах подобраны так, что выполняются условия

интерполяции:

р (х-1, У/-1 ) = /г-1 /-1 р , Р"р (х-1, У/ ) = /г-1 / р , Рг"р (хг, У/) = /

Ри р (х-1, У/-1) = /г-1 /-1 р , Р7р (хг, У1-1 ) = /г/-1 р , Р'р (хг, У/ ) = /г

Окончательная аппроксимация ^^ (х, 7]] функции Р (х,], tp ) получается путём разложения промежуточной аппроксимации Рр (х,]) в ряд Фурье по известной ортогональной и полной в L2 ( Р) двумерной тригонометрической си

' 1 р:

' 1 р:

стеме

функций Фк (х,]) = sin

Ж

р (к)

х

V

а

^т (жq (к)]) . Окончательная аппроксимация функции Р (х, ], tp ) стро

ится в виде частичной суммы этого ряда:

N

К(х,]) = !Екр-Фк (х,]), р = 0,1,...,п/.

к=1

Здесь (р(к ), (к )) (к = 1,2, ...) - пары натуральных чисел, упорядоченные так, чтобы величины

( р(к2

ч а у

( q (к )) не убывали с ростом к, N - заданное натуральное число, которое мы будем считать исходным дан-

ным. Все функции фк ( х,]) обращаются в ноль на границе прямоугольника Р. Коэффициенты Фурье Р,р вычисляются

по стандартным формулам:

Р =

1 кр

( Рр ФФ ) (Фк Фк )

4ИРр (х,])ф (х]) а-1*

п/ п/ (

пх '1

■=11=1 ч т1

ЕЕ Л РиР (х,])Фк (х,7) dxd7+\\ Р'р (х,])фк (х,7) йхй] , к = 1,2,..., N, р = 0,l,..., п/.

Здесь круглыми скобками обозначено скалярное произведение пространства Ц>(Р:

а 1

(ф,^) = Цф( x]7dxd] = | dx |ф( x,])d].

Р 0 0

Важно отметить, что Р, р вычисляются аналитически.

к р

Аппроксимация (х, 7], г) функции р (хг) получается с помощью кусочно-линейной интерполяции коэф фициентов Фурье:

' Р., - Р,

И

* ч г у

Аппроксимация /М (х, У, г) для функции / (х, у ) получается на основе аппроксимации Р4 (х, 7], г) функции р (х, ], г) с помощью преобразования координат (1):

( у - у1(х))

N ( Р — Р 1 /

Р4 (х,],г) = Е Ркр-1 + (г - гр-1) -фк (х], г е [Гр-1,гр] , р = 1,2,...,П/.

(5)

для

N

N

к=1

/N (х, У, г ) = Р

Р - Р

77 , кр кр-

Р л Н--т-

кр-1 И/

х

V (У2(х)-У1(х))' у л /

1 (' - )

ф,

х

( у - у1 (х))

' (У2 (х)- у1 (х))

г е[гр 1, гр ] , р = 1,2,...,пТ.

4

а

т

Функция *( х, у) из начального условия (4) аппроксимируется аналогично. Введём функцию ¥( х,ц) такую

что

¥

(у _ у(х)) '(у2 (х)_ у1 (х))

= *(х, у). В прямоугольнике Р вводится равномерная прямоугольная сетка точек

г*

числа г1х и

|(х) : х = * ъ = И* = а/п* , И* = 1/<, г = 0,1,...,п*, у = 0,1,...,<} . Натуральные п* считаются исходными данными. Координатные линии этой сетки делят прямоугольник Р на прямоугольники Ру = |( х, ъ) • хг-1 < х < хг, ?]у_1 < ц < цу} , которые, в свою очередь, разбиваются диагоналями на верхний и нижний треугольники Т и Т?. Узлам введенной сетки (хг) в прямоугольнике Р соответствуют узлы криволинейной сетки (х, ) = (х , у1 (х) +(у2 (х ) _ у1 (х ))) в области Q. В качестве основных исходных данных, определяющих функцию * (х, у) , задаются её значения в этих узлах *(хг, у) = * . Эти же значения будет иметь функция ¥(х, Ц) в узлах прямоугольной сетки: ¥(хг, ) = *(хг,у у ) = у/^ . В качестве аппроксимации функции ¥(х,ц) строится функция ¥ (х, ц), которая является линейной в каждом из треугольников Т и Т?, удовлетворяет условиям интерполяции ¥( хг , ) = . Построить интерполяционную функцию ¥( х, ц) в прямоугольнике Р не сложно. Приведём формулы для её вычисления.

¥( х,ц) =

Здесь

¥ и (х,1) = ^и + р (х _ хм) + у (? _ 1] _), Щ =^ р = * _ */_!; )/н

уи = (*/_1 у _*/_1 ;_1 )/ К , ¥ ? (+ Р (х _ х_1 ) + У? ) , «у = *_1 ^

Р = (*'^_l _ */_: у_1 V, У =(* VК

Функция ¥(х,ц) используется в качестве промежуточной аппроксимации функции ¥(х,ц). Для получения окончательной аппроксимации функции ¥( х, ц) полученную промежуточную аппроксимацию ¥( х, ц) разложим в ряд Фурье по той же ортогональной и полной в ¿2 (Р) двумерной тригонометрической системе функций. В качестве окончательной гладкой аппроксимации функции ¥( х, ц) выберем Ж-ю частичную сумму этого ряда Фурье

N

¥Ж (х,?) = £¥, -фк (хц). (6)

k=1

Коэффициенты Фурье ¥ к вычисляются по стандартным формулам:

(¥,ф ) = 4 ' 4

¥",(х,ц) , (х,ц)еТ;", г = 1,2,...,п*,у = 1,2,...,п*, ¥у (х,ц) , (х,фТ/,'' = 1,2,...,п*,у = 1,2,...,п*.

¥ =

т k

4 ^ 5

(Ф ,Ф)

а

{{¥(х,^ (хц)йхйЛ=—££ ((х,ц)йЫц+Л^(х,ц)ф, (хц)

аЬ г=1 ,'=1 Т„

V Т

Важно отметить, что коэффициенты ¥ к также вычисляются аналитически.

Аппроксимация (х, у) для функции *( х, у ) получается на основе аппроксимации ¥Ж (х, Т]) для функции ¥(х,ц):

* (х, у ) = ¥Ж

(у _ у(х)) ' (у2 (х)_ у (х))

N

k=1

х

(у _ у1 (х)) ' (у2 (х)_ у1 (х))

Далее вместо краевой задачи в криволинейной трапеции для уравнения теплопроводности (2), (3), (4) решается аппроксимирующая её аналогичная краевая задача

ди

N

дг

DAu

N

дг

( ^2 N

и

дuN (х, у,г) ^ ди (х, у, г ) + ди (х, у,г)

дх2

+-

дУ2

=Г (х, У,г),

(х, у)еО\у, г е( 0,т], (7)

иЛ е = 0, г е[0,т], (8)

1( х, у)еу ^ *

и (х, У,0)= Ц (х, У), (х, У )еО. (9)

Обозначим иЫ (х, у, г) точное решение краевой задачи (7), (8), (8). Сделаем в этой задаче преобразование коор-

(у - У1(х))

динат (1). Введём функцию и (х,], г) такую, что и (х, у, г) = и

' (У2 (х)- У1 (х))'

(10)

В результате замены (1) краевая задача (5), (6), (7) превращается в следующую краевую задачу в прямоугольнике Р для линейного дифференциального уравнения 2 порядка с переменными коэффициентами:

(дгг^ д. д2UN

L (UN ):

и

дг

- D

ди

дх2

■ + 2

дЩо (х)]+Д (х))+-] - (а2 (х)1+2а (х) Д (х) ]+Д2 (х)+Л2 (х)) +

+-] („2 (х)]+Д (х)) = РN (х,],г), (х,]е Р \ Г, г е( 0,т],

д] у

= 0, г е[0,Т].

и"

I (х,])еГ

UN (х,],0) = ^ (х,]), (х,])е Р.

Здесь

' 'Iх

(11)

(12) (13)

„ , Д(х)= У1'(х))

(У2 (х)-У1 (х))

(У2 (х)-У1 (х))

Л( х ) =

а(х) -(У2№( х)-УГ( х)) + 2( у2( х)- У1(х)^ Д()

„ ()=(У2 (х)-У1 (х)) (У2 (х)-У1 (х))2, А (х) = -

(У2 (х)- У1 (х)) '

УГ ( х ) 2 у1 ( х )( у2( х )- у1( х ))

(у2 (хУ1 (х))+ (у2 (х)- У1 (х))

Краевая задача (11), (12), (13) решается проекционным методом, приводящим её решение к задаче Коши для системы обыкновенных дифференциальных уравнений. Будем искать приближённое решение этой задачи в виде суммы

N

иЦ (х] ) = Ег, (г )ф (х,]).

(14)

к=1

Для определения неизвестных функций

гк

(г)

используем проекционное соотношение:

(- (и% )-РN , фт )= 0, т = 1,2,..., N . (15)

Здесь круглыми скобками обозначено скалярное произведение пространства (Р) . Заметим, что при любом значении N функция иN (х, ]) удовлетворяет граничному условию (12), поскольку этому условию удовлетворяют все пробные функции Фк ( х,]) . Подставляя в (15) выражения (14), (5) и, учитывая условия ортогональности тригонометрических функций, а также линейность дифференциального оператора -(и) (- („и + Ди2 ) = а- (и) + (и2)) ), получим систему обыкновенных дифференциальных уравнений относительно неизвестных функций Г, ( г) вида:

dr (г) N Р _ Р т Щ = Е А . г (г) +Р _1т

, ¿_1Лтк к\) т р-1 ^

к=1

т = 1,2,...,N .

к/

- (г - гр-1) при г е [гр_1, гр ] , р = 1,2,...,п/,

(16)

Здесь величины Ат, в конечном счёте представляются в виде определённых интегралов по [ 0, а], которые можно вычислить с заданной точностью 5 , например, методом Симпсона. Положительное число 5 будем считать исходным данным. Подставляя в условия (13) выражения (14) и (6), получим начальные условия для системы (16)

1

2

(18)

Гт (0)=¥т , т = 1,2,..., N . (17)

Полученная задача Коши (16), (17) может быть решена, например, методом Рунге-Кутта. Приближённое решение исходной задачи (2), (3), (4) получим с помощью преобразования координат (1):

„Ж- (х¥, ) = ^ х (у _ у1 (х)) ^ = У>(!) ф {х (у _ у1 (х)) Л

и (хуt)=и Vх' (у(х)_л(х))• г(0 фVх'(у(х)_д (х))

Метод Рунге-Кутта позволяет получить приближённое решение задачи Коши (16), (17) с заданной точностью р (считается исходным данным) на сетке = Н • 5 : Н = Т/, 5 = 0,1,..., } . Значение N подбирается автоматически с помощью правила Рунге исходя из заданной точности р решения задачи Коши. Поэтому окончательно мы по-

N (х у т\ Т

N

лучим только значения приближенного решения в узлах этой сетки и^^ (х, у, ^) (ts = 5 • -Т—, 5 = 0,1,..., ). Сходи-

мость полученного приближённого решения к точному решению U ( X, y, t) краевой задачи (2), (3), (4) при

, , nf, nj, vi^, N ^ да, р ^ 0 , S ^ 0 будет иметь место, если исходная задача устойчива, а исходные данные

(функции) достаточно гладкие.

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

We propose a new approach to the development of numerical methods for solving boundary problems in domains of complex shape with given data table. Keywords: numerical methods; boundary problems; differential equations

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

1. Самарский А. А. Теория разностных схем. Учебное пособие. - М.: Наука, 1983.

2. Тихонов А. Н. Самарский А. А. Уравнения математической физики. - М.: Наука, 1977.

3. Трубников С. В. Компьютерное моделирование. Учебное пособие для студентов вузов. - Брянск: Изд-во БГУ 2004. - 336 С.

4.Трубников С.В. Проекционный метод решения задачи Дирихле для уравнения Пуассона с таблично заданной правой частью в области сложной формы. Электронный ресурс. // Материалы VI международной научно-практической интернет конференции "Инновационные технологии обучения физико-математическим дисциплинам" (Мозырь, 25-28 марта 2014г.). URL: http://ic14.belarusforum.net.

5.Флетчер К. Численные методы на основе метода Галеркина. - М.: Мир, 1988.

Об авторе

Трубников С.В. - кандидат физико-математических наук, доцент, заведующий кафедрой информатики и прикладной математики Брянского государственного университета имени академика И.Г. Петровского, serg.trubnickov@yandex.ru