Смешанная задача для параболического уравнения

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

Аннотация научной статьи по математике, автор научной работы — Сазонов Анатолий Юрьевич

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

Похожие темы научных работ по математике , автор научной работы — Сазонов Анатолий Юрьевич

ON CLASSICAL SOLUTION OF MIXED PROBLEM FOR PARABOLIC EQUATION CONTAINING BESSEL OPERATOR ON SOME SPACE VARIABLES

The work derives the sufficient conditions on the area boundary, the operator coefficients, the right-hand side, and the initial function, under which the Fourier series is a classical solution of the mixed problem for a second order parabolic equation containing the Bessel operator on some space variables.

Текст научной работы на тему «О классическом решении смешанной задачи для параболического уравнения, содержащего оператор Бесселя по части пространственных переменных»

О КЛАССИЧЕСКОМ РЕШЕНИИ СМЕШАННОЙ ЗАДАЧИ ДЛЯ ПАРАБОЛИЧЕСКОГО УРАВНЕНИЯ, СОДЕРЖАЩЕГО ОПЕРАТОР БЕССЕЛЯ ПО ЧАСТИ ПРОСТРАНСТВЕННЫХ ПЕРЕМЕННЫХ

Клкгиитс слоаа: оператор Бссссля; смотанная задача; параболическое уравнение. Устанавливаются лостаточттые условия тта границу области, коэффициенты оператора, правую часть и начальную функцию, при которых ряд Фурье представляет классическое рентеттт-те сметанной задачи для параболического уравнения второго порядка, содержащего оператор Ьесселя по нескольким пространственным переменным.

П 1 произвольная ограниченная область, расположенная и Е» 1и прилегающая к гиперплоскостям рі 0 . ут 0. Обозначим через Iм/. часть границы области И1, лежащей на гиперплоскостях у\ 0 . ут 0, Г 1 і и — — • ЦИ ? I’1 замыкание оставшейся

где 6 некоторое положительное число.

Классической задаче (1) (3) (т 0, /.,/ самосопряженный эллиптический опера гор) посвящено болыное количество работ, напиная с самого первого результата (п I, т 0)

И.А. Стеклова [2| и до наиболее iio.iim.ix к настоящему времени результатом (т 0) В.А. Ильина [3]. В ноклассической задаче; случай т I рассмотрен в работе [1].

Общее решение задачи (I) (3) прсдставимо рядом Фурье

и котором 1’р(х) ор-гоиормироншг 11ые собственные функции, а Аг> соотнетствующие соб-

ственные значения краеной задачи:

Ь,.>г I V = 0 в области 121. Нг = 0. —-‘ <>!)>■

Через фр н /р(1) обозначены коэффициенты Фурье ршпожошш функций У (/Л’1 • — — •, /Л/г у с)у ‘ ( /Л/1/Л/„, )■

Через |) замыкание множества С] (П+ и Г°) по норме

«II//’_ ( \ч,\2ук(1х \ X / \1Уг’ 1/,/ ,1\2!

Положим я“ (12+) =Ь2.*.(П+). Замыкание подмножества 1 х и ран11]>1х пулю вблизи Г1, сходящихся но норме пространства | (П 1 ) к функции р. Для функции ^ имеем

Согласно нераиенстиу Коши Вуняконского.

У Ш — 0. Переходя к пределу при 6—>0 в равенстве (8) и неравенствах (9). (10), получим формулу (7). Ра пенсию (6) устанавливается аналогично.

«I е м м а 2. Для любой функции 2а,2

Раскрывая скобки и учитывая (12) и (13). получим

отку,ча вытекает неравенство Весса;1я (11).

. I е м м а 3. Для любой функции у €//* | (^ 1 ) « обладающей, обобщенными ‘производными второго порядка принадлежащих классу Ьч,к(И 1) справедливо неравенство вида:

И частности, утверждается сходимость числового ряда, стоящего в левой части (14). Дока з а т о л ь с, т в о. Для любой функции *(х), удовлетворяющей условиям леммы 3 справедливо тождество (12). Производя в этом тождестве интегрирование по частям с помощью леммы 1. в которой полагаем р = гр, получим

I і’р1;уггук(іх Ар І рі’Рукііх Хрірр

‘Здесь (/^з/г)р обозначает коэффициен т Фурье функции І^ір. Записывая для функции Ьу>ір пераиепстио Весселя н учитыиая раиеистио (15). получим перапеиспю (14).

. I (; м м а \. Пусть коэффициенты оператора Ьу> имеют в замкнутой обла-

сти 121 непрерывные производные до порядка я. коэффициенты Ы(х) и г(:г) до порядка ,ч — 1, « — любое целое положительное число. Пусть функция, є ІІ%.+\ (^+) удовлетворяет следующему требовал шю:

ір. Ьцгг. Ь’у,р> принадлежат пространству Н *.+(£2 ).

Тогда для функции р(х) справедливы неравенства вида: для четного в

для нечетного я

Х>^+, , начальная функция р<х), и правая

часть уравнения /(а\ I) удовлетворяют следующим требованиям:

1) коэффициенты а^(х) и с(а?) удовлетворяют условию И-эллиптичности (4) и условию с(.г) р. Ь , 1 р принадлежат про-

страпству II | (І21);

I» і І і І |2., . І «■ -ъ 111

•1) / Є НІ | 2 ( Не можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

и;> которой следует сходимость

для почти всех г € [0, Г] . Следоиателыю. существует интеграл

Из теоремы о предельном переходе под знаком интеграла Лебега вытекает сходимость второго ряда в (18).

Лемма -г) в классической формулировке.

. I е м м а 0. Пусть коэффициенты оператора /-у, начальная функция у(.т), и правая часть уравнения, /(хЛ) удовлетворяют следующим требованиям:

1) коэффициенты (121 );

1п 0, /-/;(/’у?||Ч 0 /1 i U ■ ■ ‘v г

1) f(x, t) непрерывно дифференцируема до порядка н | 1; производная, порядка ,ч I 2 принадлежит классу lj2,k

Т о о р е м а 1. Если f(x,t) = 0 коэффициенты оператора. Ьу> и (функция. (р(х) удовлетворяют условиям леммы о. то ряд (5) сходится равномерно во всем, замкнутом цилиндре Qt; о. ряды, полученные однократным почленным дифференцированием ряда (5) по I. и двукратным дифференцированием вида О2 /dXjiJXj, Иш сходятся равномерно в любой стропа внутренней подобласти Q[ CQ j. При этом сумма ряда (5) определяет классическое решение задами (I)—(3).

Д о к а з а г е л 1> с т к о. И условиях теоремы 1 ряд (5) имеет вид

Применяя к (19) неравенство Копти Вуняковского

и используя сходимость числового ряда ^ ррАр

установленную 1> лемме 5. и равно-

___ ОС I 71 —т — к I |

мерную в И 1 сходимость ряда г’р(а:)Ар 2 . установленную в |4|, получаем равно-

мерную сходимость ряда (19) в замкнутом цилиндре (^. Обозначим ряды, полученные однократным дифференцированием ряда (Г>) по /, дифференцированием ряда (Г>) вида ()2/дх<дху 1 е, г - произвольное положите. плюс чис.;ю. Докажем сходимость этих рядов при единственном условии 1;-2.к№ 1 )• 11ользуясь неравенством

СХОДИТСЯ ДЛЯ ПОЧТИ псех / € |0, / |, И час тности. ДЛЯ некоторого /■() € |0, 7’| Применяя к (21) неравенство Копти Вуняковского

п учитывая равномерную сходимость в 121 первого из рядов в правой части (22), получим

сходимость почти всюду в [0.7’| ря,д,а (21). в частности, для /о € |0,7’|. Из леммы 4, примененной к /*(а\ £)- предельным переходом под знаком интеграла Лебега следует сходимость

где 1’р(т) коэффициент Фурье функции //.(;г,/). Как и вытпе, применяя неравенство К01 ВуПЯКОВСКОГО, покажем равномерную СХОДИ МОСТ], в ().р ряда

X ‘>(•’•) I Тр(т) Не можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

Применяя теорему Фубинп и переходя к пределу при ] —*•(), полупим равенство

/ / /,И*)^(*) ), справедливы утвероюдения, теоремы. I.

Д о к а з а і е л і> с т в о. І Іокажем, ч то ряд

В силу равномерной сходимости второго ряда в (21) в произвольной подобласти 12С 121 и леммы 5 вытекает равномерная сходимость в С 0 замена »(.т,/) п(хЛ)еп1′ приводит к уравнению

в котором коэффициент с(х) —а I с(х) в силу непрерывности г(.г) на П+ отрицат(‘лсн при достаточно большой постоянной а > 0.

Автор выражает благодарность Т.Д. Воробьевой и М.В. Ворзовой за внимание к работе и обсуждения.

1. Куприянов И.А. О краевых задачах для уравнений в частных производных с дифференциальным оператором Бесселя . / ДА11 СССР. 1964. Т. 158. .4* 2. С. 275-278.

2. Стеклоо В. А. Остютшыо задачи математической физики. ТТг. 1922. Т I.

‘Л.Ильин- В.А. О разрешимости смешанных задач для гиперболических и параболических уравнений ,7 УМП. I960. ГГ. 15. Выи. 2. С. 97-151.

4. Сазонов А.Ю. О классическом решении смешанной задачи для сингулярного параболического уравнения / Дифферепциальпые уравпепия. 1990. Т. 26. Л»» 8.

5.Киприм/шо И.А. Оиигуляриыс эллиптические краевые задачи. М.: Наука, 1997.

(>. Салонов Л. ТО, Фомичева ТО. Г. О разрешимости ехгептаттттой задачи для гиперболического уравпетптя, содержащего оператор Бесселя по части пространственных переменных Междупар. конф. «Дифферепци-альпые уравнения. Функциональные пространства. Теория приближений», посвященная 100-летию со дня рождения (’..1. (Соболева. 11овосибирск: Октябрь, 2008.

БЛАГОДАРНОСТИ: Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (проекты .V» 1 ••1-01-00877, № М-01-97501).

Поступила в редакцию 21 ноября 2013 г.

OX CLASSICAL SOLI. ТЮХ OF MIXED PRO В L НМ FOR. PARABOLIC EQUATION COX’TAINIXC

bessel operator ox so mi: space variables

The work derives the sullicient conditions on the area boundary, the operator cocHieicnts, the right,-hand side, and the initial function, under which the Fourier series is a classical solution of the mixed problem for a second order parabolic equation containing (he Bessel operator on some space variables.

Key words’. Bessel operator: mixed problem; parabolic equation.

Решение систем линейных алгебраических уравнений

Главная > Решение

7.1 Метод сеток для решения смешанной задачи для уравнения параболического типа (уравнения теплопроводности)

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

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

, k =const>0.

Задано начальное условие

и заданы краевые условия первого рода

Требуется найти функцию u (x,t) , удовлетворяющую в области D (0 x a , 0 t T) условиям (7.5) и (7.6). Физически это можно представить как стержень, на концах которого поддерживается требуемый температурный режим, заданный условиями (7.6).

Рисунок 10 – Неявная схема

При проведении замены получим , т.е. k =1. Задача решается методом сеток : строим в области D равномерную сетку с шагом h по оси x и шагом  по t (см. рисунок 10).

Приближенное значение искомой функции в точке — обозначим через . Тогда ; ; i =0,1. n ; ;

j =0,1. m ; .
Заменим производные разностными отношениями

;

.

В результате получим неявную двухслойную схему с погрешностью O (  +h 2 )

.

Используя подстановку , выразим из этой схемы u i,j-1

,

где: u 0, j =  1 ( t j ) ; u n , j =  2 ( t j ) .

Получаем разностную схему, которой аппроксимируем уравнение (7.4). Эта схема (7.7) неявная, и выглядит так, как показано на рисунке 10. При построении схемы (7.7) получается система линейных уравнений с трехдиагональной матрицой. Решив ее любым способом (в частности, методом прогонки), получаем значения функции на определенных временных слоях. Так, на нулевом временном слое используем начальное условие U i,0 =f ( x i ), т.к. j =0. Эта неявная схема более устойчива для любых значений параметра >0.

Есть и явная схема (рисунок 11), но она устойчива только при , т.е. при .

Рисунок 11 — Явная схема

7.2 Решение задачи Дирихле для уравнения Лапласа методом сеток

Рассмотрим уравнение Лапласа

.

Уравнение (7.8) описывает распространение электромагнитных волн(полей). Будем рассматривать уравнение Лапласа в прямоугольной области с краевыми условиями

; ; ; ,

где -заданные функции. Заметим, что чаще всего область бывает не прямоугольной.

Введем обозначения u ij = u ( x i , y j ). Накладываем на прямоугольную область сетку ; i =0,1,…, n ; ; j =0,1,…, m . Тогда , .

Частные производные аппроксимируем по формулам

и заменим уравнение Лапласа конечно-разностным уравнением

Рисунок 12 – Схема “крест”

,

где: i =1,…, n -1, j =1. m -1 (т.е. для внутренних узлов).

Погрешность замены дифференциального уравнения разностным составляет величину О(). Выразим u i , j при h =l, и заменим систему

Получаем систему (7.10) линейных алгебраических уравнений, которые можно решить любым итерационным методом (Зейделя, простых итераций и т.д.). При этом построении системы использовалась схема типа “крест”(рисунок 12). Строим последовательность итераций по методу Гаусса-Зейделя

,

где s -текущая итерация.

Условие окончания итерационного процесса

.

Условие (7.11) ненадежно и на практике используют другой критерий

где .

Схема “крест “- явная устойчивая схема ( малое изменение входных данных ведет к малому изменению выходных данных).

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

Рассмотрим уравнение колебания однородной и ограниченной струны.

Задача состоит в отыскании функции u ( x , t ) при t >0, удовлетворяющей уравнению гиперболического типа

,

где: 0 x a ; 0 t

и краевым условиям

Заменим С на с t и получим уравнение

и в дальнейшем будем считать С =1.

Для построения разностной схемы решение задачи (7.12)-(7.14) построим в области сетку ; i = 0,1,…, n ; ; ; j =0,1,…, m ;  m = T .

Аппроксимируем (7.12) разностными производными второго порядка точности относительно шага

.

Полагая  =  / h перепишем (7.15), выразив U i , j +1. Таким образом получим трехслойную разностную схему

,

где: i =1,…, n ; j =1,…, m . Задаем нулевые граничные условия  1 ( t ) =0,  2 ( t ) =0. Тогда в (7.16) , для всех j .

Схема (7.16) называется трехслойной, т.к. она связывает значения искомой функции на трех временных слоях j -1, j , j +1.

Численное решение задачи состоит в вычислении приближенных значений решения u ( x , t ) в узлах при i =1,…, n ; j = 1,…, m . Алгоритм решения основан на том, что решение на каждом следующем слое ( j = 2,3. n ) можно получить пересчетом решений с двух предыдущих слоев ( j = 0,1. n — 1) по формуле (7.16). При j =0 решение известно из начального условия . Для вычисления решения на первом слое ( j = 1) положим

,

тогда , i = 1,2,…, n . Теперь для вычисления решений на следующих слоях можно использовать формулу (7.16).

Описанная схема аппроксимирует задачу (7.12)-(7.14) с точностью O (  + h ). Невысокий порядок аппроксимации по  объясняется грубостью аппроксимации по формуле (7.17).

Схема будет устойчивой, если выполнено условие .

Лабораторная работа № 1

Решение систем линейных алгебраических уравнений

Входные параметры: n—целое положительное число, равное порядку n системы; а — массив из n х n действительных чисел, содержащий матрицу коэффициентов системы (а(1) = а 11 , а(2) = a 12… а(n) = а n 1 , а(n + 1) = а 12 , . а(n х n) = а nn ); b — массив из n действительных чисел, содержащий столбец свободных членов системы (b(1) = b 1 , b(2)=b 2, …b(n)=b n ) .

Выходные параметры: b—массив из n действительных чисел (он же входной); при выходе из программы содержит решение системы b(l) = x 1 , b(2) = x 2 , … b(n) = х n ; error—признак правильности решения (код ошибки): если ks = 0, то в массиве b содержится решение системы, если error= 1, исходная система не имеет единственного решения (определитель системы равен нулю).

Перед обращением к подпрограмме SIMQ необходимо:

1) описать массивы а и b. Если система содержит n уравнений, то массив а должен содержать n 2 элементов, а массив b – n элементов;

2) присвоить значение параметру n, который равен числу
уравнений системы;

3) присвоить элементам массивов а и b значения коэффициентов системы следующим образом: a(l) = a 11 , а(2) = а 21 , а(3) = а 31 ,…а(n) = а n1 а(n+1) = а 12 , а(n+2) = а 22 … а(n x n) = а nn . b(1) = b 1 , b(2)=b 2, …b(n)=b n

4) проверить соответствие фактических параметров по типу и порядку следования формальным параметрам подпрограммы SIMQ. Параметры а и b — величины вещественного типа, n и error — целого типа.

Задание. Используя программу SIMQ, решить заданную систему трех линейных уравнений. Схема алгоритма приведена на рисунке 13.

Порядок выполнения лабораторной работы:

1. Составить головную программу, содержащую обращение к SIMQ и печать результатов;

2. Произвести вычисления на ЭВМ.

Пример. Решить систему уравнений

Рисунок 13 – Схема алгоритма метода Гаусса

PROCEDURE SIMQ(Nn:Integer;Var Aa:TMatr;Var Bb:TVector;Var Ks:Integer);

Var Max,U,V : Real; I,J,K1,L : Integer;

For I:=1 To Nn Do Aa[i,Nn+1]:=Bb[i];

For I:=1 To Nn Do

For L:=I+1 To Nn Do If (Abs(Aa[l,i])>Max) Then

For J:=I To Nn+1 Do

Begin U:=Aa[i,j]; Aa[i,j]:=Aa[k1,j]; Aa[k1,j]:=U;

For J:=I To Nn+1 Do Aa[i,j]:=Aa[i,j]/V;

V:=Aa[l,i]; For J:=I+1 To Nn+1 Do Aa[l,j]:=Aa[l,j]-Aa[i,j]*V;

For I:=Nn-1 Downto 1 Do

For J:=I+1 To Nn Do Bb[i]:=Bb[i]-Aa[i,j]*Bb[j];

Вычисления по программе привели к следующим результатам:

X(1)= .100000E+01 Х(2)= .200000Е+01 Х(3)= .З00000Е + 01

признак выхода 0

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

Метод квадратных корней Холецкого

Входные параметры: n—целое положительное число, равное порядку n системы; а — массив из n х n действительных чисел, содержащий матрицу коэффициентов системы (а(1) = а 11 , а(2) = a 12… а(n) = а n 1 , а(n + 1) = а 12 , . а(n х n) = а nn ); b — массив из n действительных чисел, содержащий столбец свободных членов системы (b(1) = b 1 , b(2)=b 2, …b(n)=b n ) .

Выходные параметры: b—массив из n действительных чисел (он же входной); при выходе из программы содержит решение системы b(l) = x 1 , b(2) = x 2 , … b(n) = х n ; p—количество операций.

Схема алгоритма приведена на рисунке 14.

Пример. Решить систему уравнений

Procedure Holets(n:integer;a:TMatr;b:TVector;var x:TVector;var p:integer);

Var i,j,k:integer; a11:real;

For i:=1 To n Do Begin

If i<>1 Then Begin

If a[i,i]=0 Then Begin

p:=0; error:=2; MessageDlg(‘. ‘,mtError,[mbOk],0);

For j:=1 To i Do Begin

For k:=1 To j-1 Do Begin

For i:=1 To n Do Begin

For j:=1 To i-1 Do b[i]:=b[i]-a[i,j]*b[j];

If a[i,i]=0 Then Begin

p:=0; error:=2; MessageDlg(‘. ‘,mtError,[mbOk],0);

For i:=n DownTo 1 Do Begin

For j:=n DownTo i+1 Do b[i]:=b[i]-a[i,j]*b[j];

Вычисления по программе привели к следующим результатам:

X(1)= .100000E+01 Х(2)= .200000Е+01 Х(3)= .З00000Е + 01

Рисунок 14 — Схема алгоритма метода Холецкого

Тема лабораторной работы №1 для контроля знаний проиллюстрирована контрольно – обучающей программой.

Смешанные задачи для нестационарных уравнений

Одним из базовых уравнений математической физики является одномерное параболическое уравнение второго порядка. В прямоугольнике $$ \bar_T = \bar <\Omega>\times [0, T], \quad \bar <\Omega>= \ < x\ : \ 0 \leq x \leq l \>$$ рассматривается уравнение $$ \begin \tag <1>\frac<\partial u> <\partial t>= \frac<\partial > <\partial x>\left( k(x) \frac<\partial u> <\partial x>\right) + f(x,t), \quad x \in \Omega,\ t \in (0, T], \end $$ дополненное граничными (первого рода) $$ \begin \tag <2>u(0, t) = u(l, t) = 0, \quad t \in (0, T], \end $$ и начальными $$ \begin \tag <3>u(x, 0) = I(x), \quad 0,\quad x \in \bar\Omega, \end $$ условиями.

Для простоты изложения мы рассматриваем простейшие однородные граничными условия и зависимость коэффициента \( k \) только от пространственной переменной, причем \( k(x) \geq \kappa > 0 \).

Вместо условий первого рода (1) могут задаваться другие граничные условия. Например, граничные условия третьего рода: $$ \begin \tag <4>-k(0) \pd (0, t) + \sigma_1(t) u(0,t) &= \mu_1(t),\\ \tag <5>k(l) \pd (l, t) + \sigma_2(t) u(l,t) &= \mu_2(t), \quad t \in (0, T] \end $$

К нестационарным уравнениям математической физики также относится гиперболическое уравнение второго порядка. В одномерном по пространству случае ищется решение уравнения $$ \begin \tag <6>\frac<\partial^2 u> <\partial t^2>= \frac<\partial > <\partial x>\left( k(x) \frac<\partial u> <\partial x>\right) + f(x,t), \quad x \in \Omega, \ t \in (0, T]. \end $$ Для однозначного определения решения этого уравнения помимо граничных условий (2) задаются два начальных условия $$ \begin \tag <7>u(x, 0) = I(x), \quad \frac<\partial u><\partial t>(x, 0) = V(x), \quad x \in \bar\Omega. \end $$

Разностные схемы для параболического уравнения

Явная схема

Первый шаг процедуры построения разностной задачи состоит в дискретизации области \( \bar_T \). Будем использовать равномерную пространственно-временную сетку $$ \bar\omega = \bar<\omega>_h \times \bar<\omega>_\tau, \quad \bar<\omega>_h = \< x_i = i h, \ i = 0, 1, \ldots, N_x \>, $$ $$ \bar<\omega>_\tau = \ < t_n = n\tau, \ n = 0, 1, \ldots, N_t \>$$

Для аппроксимации уравнения (1) используем явную разностную схему $$ \begin \tag <8>\frac — y^> <\tau>= — Ay^ + \varphi^n, \quad n = 0, 1, \ldots, N_t. \end $$ Здесь через \( y^n = y_i^n \) обозначено приближенное решение в узле сетки \( (x_i, t_n) \in \omega \), сеточный оператор \( A \) определен для сеточных функций \( y = 0 \) при \( x \in \partial\omega_h \) (граничные узлы сетки \( \omega_h \)) следующим образом: $$ Ay^n = -(ay_<\bar>^n)_ = -\frac<1> \left( a_ \frac^n — y_i^n> — a_i \frac^n> \right), \quad x_i \in \omega_h. $$ Для задания коэффициента \( a \) можно использовать выражения $$ \begin a_i &= k(x_i — 0.5h), \\ a_i &= 0.5(k(x_i) + k(x_)),\\ a_i &= \left( \frac<1> \int_>^ \frac \right)^<-1>. \end $$

Ключевое свойство явной разностной схемы (8) заключается в том, что она легко решается. Считаем, что приближенное решение \( y^n \) на временном слое \( t_n \) известно для всех \( x\in \bar\omega_h \). Следовательно, неизвестным в (8) является \( y^ \) для \( x\in\omega_h \). Решая (8) относительно \( y^ \) получим рекурретное соотношение $$ \begin \tag <9>y_i^ = y_i^n + \frac<\tau> \left( a_ \frac^n — y_i^n> — a_i \frac^n> \right) + \tau \varphi_i^n \end $$

Замечание

В случае постоянного коэффициента \( k(x) = a \) соотношение (9) запишется в виде $$ \begin \tag <10>y_i^ = y_i^n + F \left( y_^n — 2y_i^n + y_ \right) + \tau \varphi_i^n. \end $$ Здесь введен параметр (сеточное число Фурье) $$ \begin \tag <11>F = \frac \end $$

Таким образом, вычислительный алгоритм следующий:

1. Вычисляем \( y_i^0 \) для всех \( i = 0, 1, \ldots, N_x \)
2. Для \( n = 1, 2, \ldots, N_t \):
3. применяем (9) для всех внутренних пространственных узлов (\( i = 1, 2, \ldots, N_x — 1 \))
4. устанавливаем граничные значения \( y_i^ = 0 \) для \( i = 0 \) и \( i = N_x \).

Реализация

Файл parab1d.py содержит функцию solver_Ex_simple для численного решения одномерного параболического уравнения с однородными граничными условиями и постоянным коэффициентом \( k(x) = a \) по явной схеме:

Для ускорения расчета можно воспользоваться векторными вычислениями, т.е. заменить явный цикл

арифметикой над срезами

\( F \) — ключевой параметр при дискретизации параболического уравнения

Параметр \( F \) — безразмерная величина, которая объединяет физический параметр задачи, \( a \), и параметры дискретизации \( h \) и \( \tau \). Свойства явной разностной схемы зависят от величины \( F \).

Функция solver_Ex реализует как скалярную так и векторизованную версии. Кроме того, функция solver_Ex имеет в функцию обратного вызова так что пользователь может обрабатывать решение на каждом временном слое. Функция обратного вызова имеет вид user_action(u, x, t, n) , где u — массив с решением на n -ном временном слое, x — пространственная сетка, t — временная сетка.

Верификация

Возьмем постоянный коэффициент \( k(x) = \alpha \). Подставляя его в уравнение (1), получим функцию источника $$ f(x,t) = 10 \alpha t + 5x(l-x) $$ Легко проверить, что, если подставить \( u(x,t) \) и \( f(x,t) \) в явную схему (8), то получится тождество.

Вычисление правой части легко автоматизировать, используя sympy :

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

Для использования функций при численном решении, требуется, чтобы u и f были функциями Python. Например, точное решение должно быть функцией u_exact(x, t) , а правая часть f(x, t) . Параметры a и l в определении u и f являются символами и должны быть заменены на объекты типа float в функциях Python. Это можно сделать, определяя a и l как числа и подставляя их вместо символов. Соответствующий код может выглядеть следующим образом:

Здесь мы также создали функцию начальных данных I .

Идея теста заключается в том, чтобы наше решение точно (с машинной точностью) воспроизводилось реализованным расчетным кодом. Для этих целей в файле test_solver_ex.py мы создадим тестовые функции для сравнения точного и приближенного решений в конце временного интервала. Одну для тестирования скалярной реализации:

Вторую для тестирования векторизованной реализации:

Для запуска тестов можно воспользоваться утилитой py.test . Для этого в каталоге где собраны файлы с именами начинающимися с test_ достаточно запустить команду

В этом случае будут выполнены все функции с именами вида test_*() из всех файлов, которые также начинаются с test_ . Тестовые функции (их имена имеют вид test_*() , можно создавать в файлах с любыми названиями. В этом случае или если мы хотим запустить тесты из одного файла, команда запуска тестов будет выглядеть так

В случае удачного прохождения тестов вывод будет примерно таким

В случае провала теста будет выведено сообщение которое передается в качестве второго параметра команды assert .

Проверка порядка сходимости.

Если выбранное решение не является точным для разностной схемы, мы можем проверить порядок сходимости. Для параболического уравнения порядки сходимости по времени и пространству разные. Для численной погрешности вида $$ E = C_t \tau^r + C_x h^p $$ имеем \( r=1 \), \( p = 2 \) (\( C_t \) и \( C_x \) — неизвестные постоянные). Для упрощения, введем один параметр дискретизации \( \delta \): $$ \delta = \tau, \quad h = K \delta^, $$ где \( K \) — некоторая постоянная. Это позволяет нам выделить только один параметр дискретизации $$ E = C_t \delta^r + C_x(K\delta^)^p = \tilde \delta^r, \quad \tilde = C_t + C_xK^p. $$ В нашем случае \( r=1 \).

Используя точное решение дифференциальной задачи, мы вычислим погрешность \( E \) на последовательности уточненных сеток и посмотрим, выполняется ли \( r = 1 \). Мы не будем оценивать константы \( C_t \) и \( C_x \), так как нас не интересуют их значения.

В случае постоянного коэффициента \( k(x) = \alpha \) для выбора \( K \) воспользуемся условием устойчивости явной схемы: $$ \frac<\alpha \tau> \leq \frac<1> <2>\Rightarrow h \geq \sqrt <2\alpha>\delta^ <1/2>$$ Поэтому, выбирая \( K = \sqrt <2\alpha>\), мы обеспечим в эксперименте выполнение условия устойчивости.

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

1. Вычисление погрешностей на последовательности сеток. Для этого зададим начальный параметр дискретизации \( \delta_0 \), и выполним эксперимент уменьшая \( \delta \): \( \delta_k = 2^k \delta_0 \), \( k=0, 1, \ldots, m \). Для каждого эксперимента мы должны сохранять \( E \) и \( \delta \). Стандартно для оценки погрешности используются нормы сеточных функций: $$ \begin \tag <12>E &= \| e^n \|_2 = \left( \sum_^ \sum_^ (e_i^n)^2 h\tau, \right)^<1/2>, \quad e_i^n = u_<\rm>(x_i, t_n) — y_i^n,\\ \tag <13>E &= \| e^n \|_\infty = \max_ |e_i^n|. \end $$

Можно также использовать норму на одном временном шаге, например, в конечный момент времени \( T \) (\( n = N_t \)): $$ \begin \tag <14>E &= \| e^n \|_2 = \left( \sum_^ (e_i^n)^2 h \right)^<1/2>,\\ E &= \| e^n \|_\infty = \max_ <0 \leq i \leq N_xi>| e_i^n |. \tag <15>\end $$

2. Пусть \( E_k \) — погрешность, полученная на эксперименте с номером \( k \), а \( \delta_k \) — соответствующий параметр дискретизации. Положив \( E_k = D \delta_k^r \), мы можем оценить \( r \) сравнивая последовательные эксперименты: $$ \begin E_ &= D \delta_^r,\\ E_k &= D \delta_k^r. \end $$ Отсюда получим $$ r = \frac<\ln(E_/E_)><\ln (\delta_/\delta_)>. $$

Так как \( r \) зависит от \( k \) добавим индекс к \( r \): \( r_k \), где \( k = 0, 1, \ldots, m-2 \), если проведено \( m \) экспериментов: \( (\delta_0, E_0) \), \( (\delta_1, E_1) \), \( \ldots \), \( (\delta_, E_) \).

В нашей рассматриваемой аппроксимации параболического уравнения \( r=1 \) и отсюда значения \( r_i \) должны стремиться к \( 1 \) с ростом \( k \).

Численный эксперимент

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

Рассмотрим задачу, где \( x/L \) — новая пространственная переменная, а \( at/L^2 \) — новая временная переменная. Правая часть (источник) \( f \) отсутствует, а \( u \) ограничена величиной \( \max_ |I(x)| \). В этом случае параболическое уравнение принимает вид $$ \frac<\partial u> <\partial t>= \frac<\partial^2 u> <\partial x^2>$$ на отрезке \( [0, 1] \) с граничными условиями \( u(0) = 0 \) и \( u(1) = 0 \). Протестируем два начальных условия: разрывную площадку $$ I(x) = \begin 0, & |x — L/2| > 0.1 \\ 1, & |x — L/2| \leq 0.1 \end $$ и гладкую функцию Гаусса $$ I(x) = e^<-\frac<1><2\sigma^2>(x — L/2)^2> $$

Функции plug и gaussian в файле parab1d.py запускают решение двух соответствующих задач:

Эти функции используют функцию viz для запуска солвера и визуализации решения посредством использования функции обратного вызова для рисования:

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

Задачи

Задача 1: Задача третьего рода

Написать программу для численного решения краевой задачи для одномерного параболического уравнения с краевыми условиями третьего рода: $$ \frac<\partial u> <\partial t>= \frac<\partial^2 u> <\partial x^2>+ f(x,t), \quad 0 —>