CC BY
5МАТЕМАТИКА
УДК 517.6:517.958
К ВОПРОСУ О РЕДУКЦИИ НЕОДНОРОДНОЙ КРАЕВОЙ ЗАДАЧИ ДИРИХЛЕ ДЛЯ ВОЛНОВОГО УРАВНЕНИЯ НА ОТРЕЗКЕ
канд. физ.-мат. наук, доц. Д.Ф. ПАСТУХОВ, канд. физ.-мат. наук, доц. Ю.Ф. ПАСТУХОВ (Полоцкий государственный университет);
Н.К. ВОЛОСОВА
(Московский государственный технический университет им. Н.Э. Баумана)
Предложен алгоритм решения общей начально-краевой задачи неоднородного волнового уравнения на отрезке с неоднородными краевыми условиями. Определено понятие краевой функции. Исходная задача с неоднородными краевыми условиями редуцируется к двум простым модифицированным задачам, т.е. к задаче с модифицированной правой частью и к задаче с модифицированными начальными условиями, но с однородными граничными условиями. Получено разложение невязки задачи в самом общем виде для оптимального параметра аппроксимации разностной схемы г = 1. Формула невязки определяется производными четного порядка по времени и координате от правой части уравнения и производными четного порядка по времени от краевой функции. Написана программа на основе алгоритма редукции, решены точно и численно три тестовых примера, показывающие, что неоднородные краевые условия Дирихле сохраняют все свойства задачи с однородными краевыми условиями при использовании модифицированных условий и краевой функции.
Ключевые слова: краевая функция, модифицированные начальные условия и правая часть уравнения, неоднородно-краевая задача волнового уравнения на отрезке, согласование начальных и краевых условий.
Введение. Рассмотрим произвольную неоднородную начально-краевую задачу для волнового уравнения на отрезке с краевыми условиями Дирихле. Обычно задачи математической физики решают для однородных граничных условий, подразумевая, что возможна редукция задачи с неоднородными граничными условиями к задаче с однородными граничными условиями [1-3, 5, 14]. С другой стороны, отметим важность особенностей граничных условий в задачах математической физики, если границы имеют сложную геометрию или граничные функции принадлежат к более слабому классу функций, например, классу кусочно-непрерывных функций на отрезке. Неоднородная краевая задача эквивалентна однородно-краевой задаче с модифицированными неоднородными начальными условиями и модифицированной правой частью волнового уравнения, тогда в случае кусочно-непрерывных краевых функций можно использовать методы управляемых систем ОДУ [9, 10].
Участники семинара «Возобновляемые источники энергии» в МГУ, основанного профессором В.В. Алексеевым, д. ф. м. н., академиком РАЕН, считали, что в задачах математической физики определяющую роль играют граничные условия. В данной статье, продолжающей работу [1] на случай неоднородных краевых условий на отрезке, волновое уравнение на отрезке рассматривается в самой общей постановке задачи.
Отметим также, что уравнения математической физики, используемые в данной работе и в работах [1, 5], решаемые предложенными в них алгоритмами с двойной точностью, открывают путь для применения уравнений математической физики в совершенно новой области - стеганографии. Метод впервые предложен Н.К. Волосовой, сотрудником МГТУ им. Н.Э. Баумана, которой принадлежит смелая идея отображения пространства кусочно-непрерывных функций в пространство бесконечно-дифференцируемых функций на прямоугольнике, применяя координатные функции специального вида, но при этом сохраняются особенности постановки краевой задачи для уравнения Пуассона [4, 6-8].
В уравнениях математической физики, являющихся следствием вариационных методов и методов решения ОДУ, могут использоваться новые результаты, полученные в работах [11-13].
Постановка задачи. Рассмотрим начально-краевую задачу для неоднородного волнового уравнения на отрезке:
^ = а2 д2Щ + f (х г), х е (а, Ъ), г е (0, Т) Эг2 дх2
и (х,0) = ф(х), х е [а, Ъ] (1)
щ (х,0) = у( х), х е [а, Ъ]
и(а, г) = |(0, и(Ъ, г) = |2(г), г е [0,Т],
где а2 - квадрат фазовой скорости волны;
f (х, г) - неоднородная правая часть волнового уравнения на отрезке [а, Ъ];
ф(х), у (х) - неоднородные начальные условия;
(г), ^(О - неоднородные граничные условия; (х, г) - координата и время соответственно.
Проведем линейную редукцию задачи (1). Построим простейшую линейную функцию, удовлетворяющую краевым условиям:
Определение 1. Пусть ЭО- кусочно-гладкая граница области О , О- компакт. Дважды непрерывно дифференцируемую функцию V(О,Т) е С2 (О,Т), V : ОхТ ® Я1, непрерывную на границе области V (ЭО)е С(ЭО) и определенную во всей области задачи О :
V (х, г) = (х-^ МО + МО, (2)
^ Ъ - а ) ^ Ъ - а ;
назовем краевой функцией для исходной задачи (1), если
1) V(а, г) = | (г), V(Ъ, г) = т2 (г), г е [0, Т], V(ЭОг, г) = | (г, г), г = 17,
то есть граничная функция V(О, г), на г-й компоненте ЭОг- границы с переменными г•, равна г-й функции граничных условий V(ЭОг-, г) = | (г, г) , г = 1,1, где I - параметр связности;
2) V(х, г) является одним из 3 слагаемых численного решения задачи (1) в замкнутой области О.. Преимущество линейной краевой функции V(х, г) по координате х заключается в том что,
2 Э^(х, г)
Эх2
= 0 , т.е. краевая функция дает нулевой вклад в волновое уравнение задачи (1).
Замечание 1. Задание краевой функции в виде (2) для задачи (1) не единственно, например, опре-
\п /, \m
х - a 1 , . ( b - х I
делению (1) удовлетворяет функция V(х, t) = I-I ) +1-I МО, n, m е N.
^ b - a ) ^ b - a )
Утверждение 1. На гладких элементах границы размерности dim > 1 существует производная краевой функции по нормали щ и вдоль касательного вектора к границе г, равная соответствующим производным от граничных условий:
Эт.- (r, t) ЭV (r, + en,, t) . Эт.- (r,, t) ЭV (r, + en,-, t) _ . —
^ " ' = lim-^-^, V е С-1 (ЭО,-), ; = lim-^-^, V е Сп1 (Эй,-), - = 1, l.
Эг e®0 Эг г v Эп, e®0 Эп, n\ и> >
Доказательство. В силу того, что краевая функция дважды дифференцируема по определению 1, разложим краевую функцию в ряд Тейлора, сохраняя слагаемые со вторыми производными включительно
о. ч / ч C.ЭV(r + en,,t) §2 Э2V(r + e0n,-,t)
V(ц +щ + 5г,t) = V(ц + щ,t) + 5—-^+ ---^—^<^(0е (0,1))
Эг 2 Эг2
V (r + en,- + t5, t)-V (r + en,-, t) ЭV (ц + en,-, t) 5Э^ (ц + e0n,-, t)
^-=--1-----^
5 Эг 2 Эг2
V(r, + en, +г5,t)-V(r, + en,-,t) ЭV(r, + en,-,t) Эт.-(r-,t) _ . i _ .
lim —^----^-^ = lim-^-^ = ^ ' « V(ЭО,)е Сг1 (ЭО,),
e,5®0 5 e®0 Эг Эг
где ц +en,- +г5 - точка, удаленная от граничной точки по внутренней нормали щ (ц) на расстояние e и вдоль касательного вектора г(ц- ) в граничной точке на расстояние 5 , щ(ц ),- L г(ц ).
Аналогично получим
У ( г + гп1 + 5пг, t) = У ( г +enг•, t ) + 5
ЭУ (- + £пг, t) 52 Э2У (- + е8пг, t) Эп- 2 Эп2
2
«(9е (0,1))
У (г- + £пг- + 5пг-, t)-У (г- + епг-, t) ЭУ (г + £пг-, t) 5Э2У (г- + £0пг, t)
^-=--1-----^
5 Эщ 2 Эпг-
У (г- + еп,- + 5п,-, t)-У (г; + еп-, t) ЭУ (г + т-, ^ Эт.- (г, t) _ , 1 _ ,
11т —^---^-^-^ = 11т-^-^ = ^ " ' «У(ЭОг)е С,,1(ЭОг).
" " 5 е®0 Эпг Эпг
е,5®0 5 е®0
В точках нарушения гладкости границы векторы щ(г--) и х(г) не определены, не определены Эт (г, t) Эт (г, t)
и-,-, поэтому можно требовать только более слабых условий для краевой функции на
Эп Эх
границе У(ЭО-,0 = т- (1,t), У(ЭО)е С(ЭО). Доказательство утверждения 1 завершено.
Решение системы уравнений (1) с краевой функцией (2), согласно определению 1, разложим в сумму
и(х, г) = У(х, t) + и (х, 0 .
Подставим формулу (3) в систему уравнений (1):
иа + Уя = а2 (ихх + Ухх) +1(х,0 = а2ихх + /(х, t), хе (а,Ъ), t е (0,Г) У (х,0) + и (х,0) = ф(х), х е [а, Ъ] У(х,0) + и{(х,0) =у(х), хе [а,Ъ]
У (а, t) = М!«, и (а, t) = 0, У (Ъ, t) = |т2«, и (Ъ, t) = 0, t е [0,Т ]. Выражая из последней системы уравнений неизвестную функцию и (х, t), получим
и„ = а2и хх - У и + / (х, t) = а2ихх - Г 1 m'2(t) - I 1 т"х (t) + / (х, t), х е (а, Ъ), t е (0, Т)
(3)
Ъ-а
Ъ-а
и (х,0) = Ф(х) - У (х,0) = Ф(х) -1 | т2(0) - I тх(0), х е [а, Ъ]
Ъ - а
Ъ - а
(4)
и{ (х, 0) =у( х)-У (х, 0) = у( х )-
х - а
Д2(0) -
Ъ - х
т1(0), х е [а, Ъ]
Ъ - а ) V Ъ - ау
и (а, t) = 0,и (Ъ, t) = 0, t е [0, Т ].
Как видно из системы уравнений (4), она линейна относительно и (х, t), что дает возможность провести редукцию линейной задачи и(х,t) = и^(х,^ + и2(х,t) :
= а2Щхх - | I М^И -1 I МО + /(х, t), х е (а, Ъ), t е (0,Т) V Ъ - а ) V, Ъ - а )
и1(х,0) = 0, х е [а, Ъ]
и1{ (х,0) = 0, х е [а, Ъ]
и1(а,t) = 0, и1(Ъ,t) = 0, t е [0,Т],
и2а = а и2хх, х е (а, Ъ), t е (0, Т)
и2 (х, 0) = Ф(х) - [х-^1 Д2 (0) - I1 т (0), х е [а, Ъ]
(5)
Ъ - а
Ъ-а
иЙ(х,0) = у(х)-| IД2(0)-1 ^ IД1(0), хе [а,Ъ]
(6)
Ъ- а) VЪ- ау и2(а,t) = 0, и2(Ъ,t) = 0, tе [0,Т].
Непосредственно подстановкой убеждаемся, используя формулу и(х, t) = и1(х, t) + и2(х, t), что построчно сумма соответствующих уравнений систем (5) и (6) дает такое же уравнение из системы урав-
нений (4), в той же строке. То есть линейная редукция задачи (4) свелась к решению задач (5) и (6), а сумма решений задач (5) и (6) равна решению системы уравнений (4) ввиду ее линейности.
Таким образом, исходная задача (1) с неоднородными граничными условиями имеет решение и(х,г) = V(х,г) + Ul(х,г) + х,г), где краевая функция V(х,г) определена формулой (2) и на концах отрезка равна граничным условиям в (1), а П^х,г), ^(х,г) - решения линейных систем уравнений (5) и (6).
Отметим, что в работе [1, формулы (15) (16), (18), (19), (22), (24)] предложен алгоритм инициализации задачи на втором временном слое по начальным условиям задачи (1). Предложен и алгоритм для численного решения задач (5), (6) в [1, формулы (20), (27)] с однородными краевыми условиями на основном временном отрезке с бесконечным порядком аппроксимации путем выбора оптимального параметра аппроксимации разностной схемы, а также алгоритм масштабирования задачи, сокращающий число элементарных математических операций в [1, формулы (21), (28)] в I раз.
Задача (5) имеет однородные условия, но волновое уравнение содержит модифицированную правую часть, измененную краевыми условиями задачи (1). По сравнению с задачей (1), задача (6) имеет однородное уравнение и однородные краевые условия, но начальные условия модифицированы краевыми условиями исходной задачи (1). Таким образом, мы свели общую задачу (1) с неоднородными краевыми условиями к двум известным задачам (5), (6) с модифицированной правой частью волнового уравнения и модифицированными начальными условиями соответственно. Но задачи (5), (6) решены алгоритмами, предложенными в работе [1]. Тем не менее нужно подробнее рассмотреть особенности разностных задач (8), (9).
Для создания численного алгоритма неоднородной краевой задачи по аналогии с формулой (2) построим дискретный аналог краевой функции на равномерной сетке для неоднородной краевой задачи:
VП = V (хт, п) = Iт2(гп)+I т1(1п) =
Ъ-а
Ъ-а
тН ] (Ь-а-тН] ,п.
112 (пт) +1 —--11 (пт), (7)
Ъ-а
Ъ-а
тттт ¡г—- , (Ъ - а) Т
т = 0, М, п = 0, N, Н =-, т = —, хт = а + тН, г„ = пт.
М N т
Используя вспомогательную функцию (7), аналогично задачам (5), (6) получим их разностный аналог с модифицированными неоднородными условиями (для шаблона крест)
П п+1 + П п-1 2П п
и1т + П1т - 2П
'1т
1т = а2 П1т+1 + П1т-1 '
п
1т
Н2
+ f (а + тН, пт), т = 1, М -1, п = 1, N -1
f (а + тН, пт) = f (а + тН, пт) -
Пш = 0, Пьт = 0, т = 0М
щп = 0,ПМ = 0, п = 0,N.
тН Ъ - а
12(пт) -
Ъ - а - тН Ъ - а
т1(пт) ° /т, Н=
(Ъ - а) М
-, т = — N
(8)
П п+1 + П п-1 2П п
П2т + П- 2П-
2т
--т = а2 П2т+1 + П2т~1 -2П2т , Н = , т = Т, т = 1,М-1, п =
т2 Н2 М N
П2т =ф(а + тН) = ф(а + тН)-(_т^]|2(0)-[ -Ъ—а—— ||х(0) = фт, т = 0,М
Ъ - а
Ъ-а
Ъ-а
П2т = 0, у(а + тН) = у(а + тН)-1 112(0)-( Ъ ^ тН 1М0) = ут, т = 0,М
Ъ - а
П20 = 0, П2М = 0, п = 0, N.
(9)
щ п+1 + щ п-1 - 2, п ,п + п -2, п - (и— \ Т _ _
ит + ит 2ит _ „2 ит+1 + ит-1 2щт + /п н = (Ъ а) т = — т = 1М - 1п = 1N -1
т, М ' N' , , ,
т
2
=а
Н
2
хт = а + тН, гп = ит, ^ = Да + тН,пт)-[ тН |^(пт)-[Ъ—а—— ||1(пт)
^ Ъ - а) ^ Ъ - а )
и(а + тН, 0) = и0т = ф(а + тН) = ф(а + тН) -(112 (0) -[ —— 111(0) = фт, т = 0, М
Ъ-а
Ъ-а
и(а + тН,т) = ° ф1(а + тН), т = 0, М, ф1(х) = ^(ф(х), у(х)) и(а, пт) = и° = 0, и(Ъ, пт) = иМ = 0, п = 0, N
у(а + тН) = у(а + тН) - [ 112 (0) -( Ъ ? тН 111 (0) = ут, фт(х) = ^(фт(х),ут(х)).
Ъ - а
Ъ-а
2
т
1 п = Уп + ип = Уп + и п + и п
ит = Ут + ит = Ут + и 1т + и 2т,
где Ут, ит, ит, и2т - решения задач (7), (10), (8) и (9) соответственно.
Второй временной слой ит , согласно (10), с использованием модифицированных фт, ут условий требует задания и1т, ит аналогично формулам (15), (18) из [1, с. 175, 176], связывая три (четыре) последовательных временных слоя. Для однородного волнового уравнения (9) получим систему линейных уравнений с трехдиагональной симметрической матрицей, решаемой методом прогонки с оптимальным параметром аппроксимации г = 1 [1, с. 175, формула (15)]:
и2т-1 - 4и2^ + и^ = "2фт " 2хут ° Рт + О (х3 ), Ат = 1, Вт = 1, Ст = 4, т = 1, М -1, и20 = 0, и2М = 0, формулами прогонки вперед [1, 5]
В А V - V _
1т = С -Ат. , Пт = "-"А1 - т , т = 1, М -1, 10 = 0, П0 = и 20 = 0,1м = 0, V М = и2М = 0
Ст Ат1т-1 Ст Ат1т-1
и формулами прогонки назад в [1, с. 175, формула (16)]
и2т =1ти2т+1 +Пт, т = М -1,1. (11)
Для оптимального параметра г = 1, связывая пять последовательных временных слоев, выражая их друг через друга, оставляя в записи только 2 первых слоя, аналогично формуле (18) из [1, с. 176] получим пятидиагональную СЛАУ:
и2т-2 - 2 и 2т-1 +10и 2т - 2 и 2т+1 + и 2т+2 = Фт + Фт-1 + Фт+1 + 3хУт + О (^ )° Рт + О (х4 ),
99 _ _ _ _ -
Ат1 = 1, Ат2 = - ^ Вт1 = - Ст = -10, Вт2 = 1, рт = Фт + Фт-1 + Фт+1 + 3хУт, т = ^ М - 2. Получим аналогично формуле (19) [5; 1, с. 176]) формулы прогонки вперед
Т = В1т + А2т12т-1 + А1т11т-212т-1 у = В2т
11т = „ . л л .л .л , 12т ="
Ст А1т11т-211т-1 А1т12т-2 А2тУ1т-1 Ст А1т11т-211т-1 А1тУ2т-2 А2т11т-1
V = А1т11т-2пт-1 + А1тпт-2 + А2тпт-1 рт ^ = 2 М - 2 Ст - А1т11т-211т-1 - А1т12т-2 - А2тУ1т-1
110 = 0, 120 = 0, 1ц = 0, 121 = 0, П0 = и20 = 0,
(12) 1 (11) 1 1 (12) 1 (11) 1 п1 = и 21 = и 21, VМ = и2М = 0, ПМ -1 = и 2М -1 = и 2М -1 ,
(12) 1 (11) 1
где запись и 21 = и 21 читается так: значению и121 в формуле (12) присваивается старое значение и121 из формулы (11).
Тогда формулы прогонки назад аналогично [1, с. 176, формула (19б)] имеют вид
11М -1 = 12М -1 = 11М = 12М = 0 ,
и2т =11ти2т+1 +12ти2^+2 +пт, т = М - 2,2. (12)
Проверим, что указанные коэффициенты на границах дают краевые условия Дирихле:
и20 =110и21 +120и22 +п0 = U2°, и21 =111и22 + 121и2З +П1 = и2Ъ и2М-1 = 11М-1и2М + 12М-1и2М +1 +ПМ-1 = и2М-1, и2М-2 = 11М-2и2М-1 +12М-2и2М +пМ-2,
и2М = 11Ми2М+1 +12Ми2М+2 +п М = и2М .
9
Коэффициенты матриц СЛАУ Ат = 1, Вт = 1, Ст = 4 в формулах прогонки (11) и Ат1 = 1, Ат2 = —, 9
Вт1 =-~, Ст =-10, Вт2 = 1 в формулах (12), удовлетворяют условию корректности формул прогонки
([1, с. 177, утверждение 3]) и совпадают с аналогичными коэффициентами, корректность которых доказана в работе [1].
Запишем явную разностную схему для решения однородного уравнения, учитывая оптимальный параметр г = 1 из [1, с. 178, формула (20)] с бесконечным порядком аппроксимации:
и2пт+1 = и2пт+1 + ^1-1 - ^Ц-, п = , т = 1, М -1 . (13)
Возможен также алгоритм масштабирования (укрупнения) ячеек сетки с параметром г = 1 и коэф-
2
фициентом масштабирования I, сокращающим число и время вычислений в I раз.
Сначала инициализация проводится с шагом х по времени по формулам (11), а затем по формулам (12), далее по явной формуле (13), с минимальными шагами сетки (к,х) решается волновое уравне-
ние и2"т+ = и2т+1 + и2""— - и2"1, п = 2,1, т = 1, М -1 во временном промежутке t = [2х, I * х]. Среди решения в слое ^(тк, I * х), т = 0, М и среди начального слоя и2(тк,0), т = 0, М выбираются узлы
более редкой сетки хт1 = а + к * I * т1, т1 = 0, М /1, и решение на редкой сетке и2(а + т1* к * 1,0), и2(а + т1* к * I, I * х), т1 = 0, М /1.
Далее используется формула (13) с крупным вектором шага (I * к, I * х)
1 = г = х2а2 /к2 ^х= к/а, и"1 = и2"11+1 + и2"11—1 -и"1, п1 = 2,N/1, т1 = 1,М /1-1. (14)
Рассмотрим вторую задачу инициализации г = 1 для системы уравнений (8) свяжем условия для 3 временных слоев х,0) = ф(х) ° 0, и^ (х,0) = у(х) ° 0, / (х, t) Ф 0, тогда на этапе инициализации по трем (четырем) временным слоям (формула (22) в работе [1, с. 179]) имеем
1 1 п" - 2 0 - 2 - 2 -
и1т-1 - 4и1т + и1"+1 = -"х + и1т = - 1х , Ат = 1 Вт = 1 Ст = 4, ¥т = -"х , 1 = 1 М - 1 (15)
Воспользуемся затем формулами прогонки (11) с учетом (15) для решения и", т = 1, М -1,
и110 = и11М = 0 . Перейдем ко второй задаче инициализации для 4 временных слоев, аналогично формуле (24) из [1, с. 180]:
" 9 " " 9 " " (- 9- - -1 2 I 4\
и1т-2 - 2 и1т-1 + 10и1т - 2и1т+1 + и1т+2 =-1 /т,2 - 2 /т,1 + /т-1,1 + -/1+1,1 Iх + О(х ) (16)
с коэффициентами
9 9 ___ „ „ (-- 9
Ат1 = 1 Ат2 =-^ Вт1 =-2, Ст =-10, Вт2 = 1 рт =-^ /т,2 - 2 /т,1 + -/т-1,1 + -/т+1,1) х2, 1 = 2,М - 2,
" (16) 1 (15) 1 " (16) 1 (15) 1
110 = 0,120 = 0,= 0,121 = 0, п0 = и"0 = 0, V" = и 11 = и 11, VМ = и11М = 0, VМ —1 = и 1М-1 = и 1М-1,
используя формулы прогонки (12), получим второй временной слой задачи (8) и", т = 0, М .
Наконец, нужно написать основную рекуррентную разностную формулу для решения неоднородного уравнения с нулевыми модифицированными начальными условиями. Для параметра г = 1 имеем по формуле (25) из [1, с. 180])
иС1 = и".! + и 1-1 - и"- + /т;х2 + Л (и"), п = , т = 1, М -1, (17)
Э2У„
/ = / -(у ) " 'т'п = 0
./т,п ./т,п \ут,п)^' э 2
Утверждение 2. Остаточный член погрешности в формуле (17) зависит от правой части волнового уравнения и краевой функции и представим в виде
¥ о ( к-1 л2(к-1) Г
Я(Пп) = V 2 т2к^а2' Э /тп
Я(Щт) = 5(2к)![т 5а Эх2'Эг2(к-т-1)
)
2т2к Э2к^
к=1 (2к)! Эг2к
- 5
Доказательство. Воспользуемся формулой невязки в общем виде для одномерного волнового уравнения на отрезке, полученной в работе [1, с. 180, формула (26)], с учетом формулы (17), имеем
, 2=1 тп+1 Пп
ттпп = т тп +т т" т тп-1 + г т2 + 2т
П1т = П1т+1 + П1т-1 П1т + 1т,пт + 4!
4 ( Э2 Э2
2 Э /т,п Э /т,п
а ----+
2~] 2т6
Эх'
2
Эг
2
(4
6!
Э4 / Э4 / Э4 / ]
4 Э /т,п 2 Э /т,п Э /т,п
Эх4
+а
Эх2Эг2
Эг4
8!
(
Э° / э6 / Э6 / Э6 / ]
6 Э /т,п 4 Э •/m.n 2 Э /т,п э /т,п
- + а
+ а
Эх6 Эх4Эг2 Эх2Эг4 Эг6
2т
10 ( . Э8
10!
Э8 / э8 / Э8 / Э8 / Э8 / ]
8 Э /т,п 6 Э •/m.n 4 Э •/m.n 2 Э /т,п э /т,п
+а
+а
+а
Эх8 Эх6Эг2 Эх4Эг4 Эх2Эг6 Эг;
2т
12 ( э™/ эю / Э10/ Э™/ ЭЮ / Э™7 ]
Э /т,п . 8 Э /т,п . 6 Э /т,п . 4 Э /т,п . 2 Э /т,п . Э /т.:
12!
10 ./т,« 8 ./т,« 6 ./т,« 4 ./т,« 2 а -тт;—+ а —-—— + а —-—— + а -;—— + а
Эх
10
Эх8Эг2 Эх6Эг4 Эх4Эг6 Эх2Эг8 Эг
10
+... =
= 77п + п" - ип~1 + / т2 | * = П1т+1 + П1т-1 П1т + /m.nт + 4!
2т4
' Э2 / э2/
2 Э /т,п Э /т,п
а -г—
Эх2
Эг2
2Г_ 6!
6 ( Э4 / Э4 / Э4
4 т,п 2 т,п т,п
Эх4
+а
Эх2Эг2 Эг4
2т
.8 ( Э6 / Э6 / Э6 / Э6 /
6 т,п 4 т,п 2 т,п т
8!
+а
+а
6 ^ ] т,п
Эх6 Эх4Эг2 Эх2Эг4 Эг6
2т
10
(8
10!
8 Э fm.n 6 Э fm.n 4 Э fm.n 2 Э fm.n Э fm.n
а --—+ а —-—- + а —;—- + а —-—- + ■
Эх8 Эх6Эг2 Эх4Эг4 Эх2Эг6 Эг8
2т
12
12!
Э10 / эш / Эш / Эш / Эш / Эш /
10 Э /m.n 8 Э /m.n 6 Э /m.n 4 Э /m.n 2 Э /т,п Э Л а -тт—+ а —-—— + а —-—— + а —;—— + а —-—— +
10
10
10
10
10
Эх
10
Эх8Эг2 Эх6Эг4 Эх4Эг6 Эх2Эг8 Эг
10
-т
2 Э_^тп Эг2 '
2т
4 Э4V 2т6 Э6V 2т8 Э8V 2т10 Э^ 2т12 Э12V
т,п 2 т,п 2 т,п 2 т, п 2 т
4! Эг4 6! Эг6 8! Эг8 10! Эг10 12! Эг
12
= П1т+1 + П1пт-1 - ПЪ- + /т,пт2 + 5 / ог ч, т2к 5
2 Г 2к к-1 21 Э2(к-1) /т
=2 (2к)!
I=0
Эх 21 Эг2(к-
-5
2т2к Э2^т
к=1 (2к)! Эг2к '
(18)
где /тп = /тп -1 —- 112 (гп) + [ Ъ хт 1| (гп) - модифицированное узловое значение правой части.
Ъ - а ) Утверждение 2 доказано.
Ъ-а
Аналогично выражению (14) можно провести укрупнение шага в неоднородной модифицированной правой части волнового уравнения для системы уравнений (8). Сначала инициализация проводится
+
+
+
+
+
+
+
+
+
+
с шагом t времени по формулам (15), а затем по формулам (16), далее по явной формуле (18) с минимальными шагами сетки (h,t) решается волновое уравнение с невязкой общего вида (18)
= U1L1 + Um-1 - U^"1 + JZtt2 + R(щт), n = 2N, m = 1,M -1, (19)
во временном промежутке [0, l * t], n e 2, l.
Среди решения в слое U\(a + mh,l*t), m = 0,M и среди начального слоя U\(a + mh,0), m = 0,M
выбираются узлы более редкой сетки xmi = a + hl * mi, mi = 0, M /l, tn\ = t * l * ni, ni = 0,1, тогда получаем решение на редкой сетке
{^(a + mi* h*l,0), U1(a + mi* h*l,l*t), mi = 0,M /l} = {u1(a + mi* h*l, tn1), mi = 0,M /1, ni = 0,i}.
2 2 2
Далее используется формула (19) с крупным шагом (l * h, l * t),1 = z = t a / h ^ t = h / a :
Uimr=Uim\+i+Ui^i-i - ui-1+Tmiinií212+R (ui), ni=xNTI, m=i, M / I -1. (20)
В последних формулах (19), (20) невязка R (uI ) определяется двумя последними суммами из выражения (28) в [1, с. 181].
Рассмотрим тестовый пример 1:
utt = uxx + sin(t)sin(x), x e (0, p), t > 0 u (x,0) = cos(x), x e [0, p]
(2i)
ut (x,0) = sin (2x), x e [0, p]
u(0, t) = cos(t) = m1(t), u(p, t) = - cos(t) = m2 (t), t > 0. u( x, t) = У (x, t) + U1 (x, t) + U2 (x, t) = cos( x) cos(t) + i (sin(t) -1 cos(t)) sin( x) + sin(2x)sin(2t). (22)
Можно проверить, что решение системы уравнений (21) описывается формулой (22), т.е. удовлетворяет волновому уравнению и 4 условиям (21).
К неоднородной краевой задаче (21) можно применить редукцию к однородной задаче, описанную формулами (1)-(6). Воспользуемся формулами (2), (5), (6) с краевой функцией:
У (x, t) = (l-a) m2(t) + (b-a) mi(t) = (Pp ) (- cos(t) ) + ^j cos(t) = - 2xj cos(t) и с модифицированной правой частью
2 x j
f (x, t) = f (x, t) - Vtt (x, t) = - ^—J cos t + cos t + sin x sin t, - Vtt (x, t) = V(x, t)
и точным решением - формула (22), используя формулу (18), получим
( 1 л J6 „8 10 £ 12 Л ( _ 2 „4 „6 Л
nn+l = jjn + nn _ rrn-1 + 2 f U1m ~U1m+1~r U1m_1 U1m ^ 2fm,n
1 2 0t4 ^t" ct10 6t —t2 _ 2—+ 3--4—+ 5---+
2! 4! 6! 8! 10! 12!
v /
+ 2Vm,n
t t t
2! 4! 6!
V J
1 ^ t2 t4 t6 t d. , t . 12 0t4 _ t6 T8 T10 6t12
1 _cos(t) =---+--... ^--(1 _cos(t)) = — sint =—t2 _2— + 3--4— + 5---+...
2! 4! 6! 2 dtK ' 2 2! 4! 6! 8! 10! 12!
UC1 = Ufm+1 + U1m_1 _ U^1 +1 sin(t)fmn + 2Vm,n (1 _ COS(t)) =
= U1nm+1 + U1nm_1 _ U1n_1 +1 sin(t) sin(mh) sin(n t) + 2 _ ^ ^^ Jj J cos (nt) (1 _ cos(t)), (23)
n = 2, N, m = 1, M _ 1.
В первом примере по формуле (23) для бесконечного ряда получена производящая функция т 8ш(т)/тп + 2Vm n (1 - еов(т)), состоящая из двух слагаемых. Первое совпадает с аналогичным примером
с той же неоднородной правой частью волнового уравнения и однородными краевыми условиями в работе [1, с. 181]. Второе слагаемое своим существованием обязано краевой функции из первого примера. Предположим, что аналитическое решение каждой частной задачи (18), (19) представимо в виде бесконечного медленно сходящегося ряда для каждой из подзадач с модифицированными неоднородным уравнением или с модифицированными начальными условиями. Тогда экономнее численно решать каждую из модифицированных подзадач, сравнивая их сумму и сумму с краевой функцией с известным общим точным решением. В таблице 1 по порядку записаны численные решения частных задач (8), (9), краевая функция и сумма всех трех величин соответственно
П1(х,г), П2(х,г), V(х,г), П1(х,г) + П2(х,г) + V(x,г), г = 25.1327412287183.
Таблица 1
П1( х, г) П2( х, г) V (х, г) П1 (х, г) + П 2 (х, г) + V (х, г)
0.000000000000000Е+00 0.000000000000000Е+00 1.00000000000000 1.00000000000000
-3.88322207745093 0.151056516295153 0.800000000000000 -2.93216556115578
-7.38632732196183 0.209016994374947 0.600000000000000 -6.57731032758688
-10.1664073846305 0.187785252292473 0.400000000000000 -9.57862213233804
-11.9513286589662 0.109016994374947 0.200000000000000 -11.6423116645913
-12.5663706143592 -7.216449660063518Е-16 0.000000000000000Е+00 12.5663706143592
-11.9513286589662 -0.109016994374948 -0.200000000000000 -12.2603456533412
-10.1664073846305 -0.187785252292473 -0.400000000000000 -10.7541926369230
-7.38632732196184 -0.209016994374948 -0.600000000000000 -8.19534431633678
-3.88322207745093 -0.151056516295154 -0.800000000000000 -4.83427859374609
0.000000000000000Е+000 0.000000000000000Е+000 -1.00000000000000 -1.00000000000000
В первом столбце таблицы 2 с равномерным шагом указаны координаты узлов. Во втором столбце записаны значения решения согласно формуле (22) в узлах координатной сетки
и(х, г) = V(х, г) + П1 (х, г) + П2 (х, г) = ео8(х) ео8(г) +1 (81п(г) - г ео8(г)) 8ш(х) + 51п(2х)51п(2г). Таблица 2
х ехасг (П^х, г) + П2(х, г) + V(х, г)) питепса1
0.000000000000000Е+000 1.00000000000000 1.00000000000000
0.314159265358979 -2.93216556115578 -2.93216556115578
0.628318530717959 -6.57731032758688 -6.57731032758688
0.942477796076938 -9.57862213233805 -9.57862213233804
1.25663706143592 -11.6423116645913 -11.6423116645913
1.57079632679490 -12.5663706143592 -12.5663706143592
1.88495559215388 -12.2603456533412 -12.2603456533412
2.19911485751286 -10.7541926369230 -10.7541926369230
2.51327412287183 -8.19534431633677 -8.19534431633678
2.82743338823081 -4.83427859374609 -4.83427859374609
3.14159265358979 -1.00000000000000 -1.00000000000000
В третьем столбце программа вычисляет сумму краевой функции и численных решений задач (8), (9). Программа с параметрами п = 100, I = 10, т = 8, г = п-т-т = 25.1327412287183 возвращает норму Че-бышева для невязки задачи (7). Норма относительной погрешности имеет порядок 1Е-15, что соответствует двойной точности решения. Последние столбцы таблиц 1 и 2 совпадают.
Относительная норма Чебышева равна 1.119445273655962Е-015. Среднее арифметическое от модуля численного решения в конечный момент времени по всем узлам составляет 7.93409415003782. В примере 1 краевые условия Мг) = ео8(г), ^(г) = - со8(г) линейно зависимы.
Тестовый пример 2.
utt = uxx +1 sin(x), x e (0, p), t > 0 u (x,0) = sin(x /2), x e [0, p]
ut (x,0) = "~c°s(x/2), xe [0,p] u(0, t) = sin(t /2) = mj(t), u(p, t) = c°s(t /2) = m2(t), t > 0. Точное решение примера (24) есть
u( x, t) = (t - sin(t)) sin( x) + sin( x /2) c°s(t /2) + c°s( x / 2) sin(t / 2)
с краевой функцией
v (x, t) = ( J^ j M2(t) + [ j mi(t) = ( p j c°s(t/2) + (^ j sin(t/2)
и модифицированной правой частью
f (x) = f (x, t) - Vtt (x, 0) = t sin( x) +—c°s(t / 2) +
4p
p-x 4p
sin(t /2),
а также с модифицированными начальными условиями:
p-x
(24)
j(x) = j(x) - V(x, 0) = sin(x /2) —, y(x) = y(x) - Vt (x,0) = -c°s (x / 2) -
p 2 v 2p
Тогда модифицированные правую часть волнового уравнения и начальные условия для второго тестового примера подставим в формулу (18):
UC = Uim+1 + U1m-1 - U1m + 2 fm
\
(„. 2 „.4 „.6 „.8 10 _12
t t t t t t
~V. Ti "6! 8 10! 12! ...
v у
+ 2Vm
( х^ -,2
468 t + t t +
468
222. 244. 266. 288.
10
12
21010. 21212.
+...
=uim+i+urn- - urn+2 mn (1 - c°s(t))+2m - c°s [=
= Uim+1 + Uim-1 - ит + 2(nt) sin(mh) (1 - c°s(t)) + 2 ((—j c°s(nt / 2) + ^ p mhj sin(nt / 2)J - cos (^J
2||, (25)
n = 2, N, m = 1, M -1.
Производящая функция 2 /тп (1 - ео8(т)) + 2Утп(! - -Г || в формуле (25) состоит из двух слагаемых: первое слагаемое обусловлено правой частью волнового уравнения, второе определяется краевой функцией.
В таблице 3 собраны результаты численного решения примера 2 с параметрами программы п = 100, I = 10, т = 8, г = п• т-т = 25.1327412287183.
Таблица 3
x exact (Ui(x, t) + U2(x,t) + V (x, t)) numerical
0.000000000000000E+000 -4.898425415289509E-016 -4.898425415289509E-016
0.314159265358979 7.92287861994210 7.92287861994209
0.628318530717959 15.0816716382986 15.0816716382986
0.942477796076938 20.7868052690006 20.7868052690006
1.25663706143592 24.4904425702249 24.4904425702249
1.57079632679490 25.8398480099049 25.8398480099049
1.88495559215388 24.7116743123074 24.7116743123074
2.19911485751286 21.2238212934494 21.2238212934494
2.51327412287183 15.7237111602188 15.7237111602188
2.82743338823081 8.75413249549701 8.75413249549700
3.14159265358979 1.00000000000000 1.00000000000000
Относительная норма Чебышева равна 3.004681541631892Е-015. Среднее арифметическое от модуля численного решения в конечный момент времени по всем узлам составляет 16.5534985368844.
Норма погрешности имеет порядок 10-15, т.е. программа работает с использованием алгоритма согласно формулам (11)-(20) с двойной точностью. В примере 2 краевые условия ^(г) = со8(г /2), ^(О = /2) линейно независимы.
Тестовый пример 3.
utt = uxx + cos(3t)sin(3x), x e (0, p), t > 0 u (x,0) = sin x, x e [0, p]
ut (x,0) = "3"sin(J, xe [0,p]
u(0, t) = 0 = m1 (t), u(p, t) = - sin |J = m2(t), t > 0.
(26)
Его точное решение есть
u( x, t) = sin (^TJ sin ( 3tJ +1 sin (3x) sin (3t) + sin x cos t
с краевой функцией V(x,t) = —-JMt) + —-JMi(t) = Jsin^и модифицированной правой
--9x (3t Л
частью f (x) = f (x, t) - Vtt (x,0) = cos(3t)sin(3x)--sin I — I, а также с модифицированными начальными
4p
3 . (3xЛ 3x
условиями j(x) = j(x) - V(x,0) = sin x, y(x) = y(x) - Vt (x,0) = — sin I — I +--.
2 ^ 2 ) 2p
Тогда найденные модифицированные правую часть волнового уравнения и начальные условия для третьего тестового примера подставим в формулу (28):
Т7"+1_гrn ,Tjn _rjn-l ,r)f U1m ~U1m+1^~ U1m-1 U1m ^ 2 Jm,n
(/,Л2 2
- 2 2! 4!
32t4 34t6 , 36t8 38t10
+3
6!
-4
8!
- + 5-
10!
-6
310t12 л -+...
12!
+2Vm
3 ]212 - (3J4t4++(316t6 - (3]8 tL+(3 j10 -13 j12^ + Л 2j 2! 12j 4! + +12J 6! Uj 8! + UJ 10! UJ 12! +.
: Ufm+1 + U1nm-1 - un- + fm,n (3t sin(3t) - 8t2 ) + 2Vm,n (1 - COs (yJJ
(27)
Поскольку
(3t)2 (3t)4 (3t)6 t — , . 1 . .2 (3t)4 (3t)6 (3t)8
1 - cos(3t) = ^—-—-+-—-—...« -—(1 - cos(3t)) = — (3t)2 - 2^-+3^— 4^- + 2! 4! 6! 2 dtK ' 2V ' 4! 6! 8!
(3t)10 (3t)12 2 t2 0 3214 , 34t6 „ 3618 . 38110 Л 310112 3t
+5-—---6-—— +... = 4t2 +--2-+ 3--4-+5--6-+... = — sin(3t) ^
10! 12! 2! 4! 6! 8! 10! 12! 2
t2 32t4 34t6 36t8 38t10 ^--2-+ 3--4-+ 5--6
2!
4!
6!
8!
10!
10 12
3 t 3t 2
-+... = — sin(3t) - 4t2.
12! 2
3 Г i! - f 3 f +f 3 Г ± - f 3 )8 i* + f 3 f - f 3 f il+..., - cos f 3t
2J 2! ^2J 4! f 2) 6! f 2J 8! f 2J 10! ^2J 12! f 2
Производящая функция fmn (3t sin(3t) - 8t2) + 2Vmn - cos f3j)j в формуле (27) состоит из двух слагаемых, первое слагаемое обусловлено правой частью волнового уравнения, второе соответствует
+
краевой функции. В таблице 4 показаны результаты численного решения примера 3 с параметрами программы п = 200, I = 20, т = 4, г = птт = 12.5663706143592.
Таблица 4
x exact (U^x,t) + U2(x,t) + V(x,t)) numerical
0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000
0.314159265358979 0.309016994374945 0.309016994374945
0.628318530717959 0.587785252292470 0.587785252292472
0.942477796076938 0.809016994374946 0.809016994374947
1.25663706143592 0.951056516295155 0.951056516295152
1.57079632679490 1.00000000000000 0.999999999999998
1.88495559215388 0.951056516295155 0.951056516295152
2.19911485751286 0.809016994374947 0.809016994374948
2.51327412287183 0.587785252292471 0.587785252292473
2.82743338823081 0.309016994374946 0.309016994374948
3.14159265358979 8.572244476756629E-016 7.347638122934264E-016
Относительная норма Чебышева равна 7.209524148017449E-015. Среднее арифметическое от модуля численного решения в конечный момент времени по всем узлам составляет 0.631375151467504.
Норма погрешности имеет порядок 10-15, т.е. программа работает с использованием алгоритма согласно формулам (11)-(20) с двойной точностью. В примере 3 задано только одно ненулевое краевое условие
m\(t) = 0, m2(t) = -sin^~~j . Таким образом, три рассмотренных примера исчерпывают все различные
случаи задания неоднородных условий, показывающие, что программа, написанная по алгоритму формул (11 )—(20), численно решает начально-краевую задачу волнового уравнения на отрезке с двойной точностью.
Для классического решения задач математической физики необходимо также согласовать начальные и краевые условия [3; 4, с. 43], которые заключаются для постановки задачи (1) в выполнении условий
У (a, t) = mi(t), У (b, t) = m2(t), j(a) = u(a,0) = mi(0), j(b) = u(b,0) = m2(0). (28)
Отметим, что в первых двух примерах условие согласования классического решения выполнено на двух концах отрезка и, как следствие, достигнуто меньшее значение нормы погрешности, чем в третьем примере, в котором условие согласования выполнено на левом конце отрезка, но не выполнено на правом его конце: j(a)-|J.1(0) = 0, j(b)-|J.2(0) = 1 ^ 0. В результате норма погрешности, полученная
программой, больше примерно в 10 раз, чем в первых двух примерах, тем не менее устойчивость всех указанных алгоритмов разностных схем имеет место и в случае рассогласования начальных и граничных условий.
Программа написана на языке FORTRAN и использует входные функции тестового примера 2. Опишем входные функции программы, численно решающей произвольную неоднородную начально-краевую задачу для волнового уравнения на отрезке:
fy(x, t) = f (x, t) - vtt (x, t) = f (x, t) -1 -—a I|m2 (t) + | -—- IM-1 (t) - модифицированная правая часть
V b - a ) V b - a )
волнового уравнения;
u0( x) = j( x) -1 I ^2 (0)+1 I im (0), u1(x) = y( x) -1 I m 2 (0)+1 I im (0) - модифици-V b - a ) V b - a ) v b - a ) v b - a )
рованные условия для начального смещения точек отрезка струны и модифицированные условия для
начальной скорости точек струны.
Сумма функций fan(x, t) + ff0(x, t) + f 2(x, t) в программе равна известному точному решению. Если указанные функции найдены, то их записывают в программу (для определения нормы невязки задачи по Чебышеву). Первые два слагаемых представляют частные решения задачи (1), соответствующие неоднородным начальным условиям и неоднородным граничным условиям, третье слагаемое соответствует неоднородной правой части волнового уравнения.
Программа написана на языке FORTRAN, ее особенностью является m = 4k, k е N, быстродействие (время работы) составляет 0,01 с.
program wave;integer(8), parameter::n=100,n1=10,ll=n/n1;integer(8)::i,j,k;integer(8),parameter::m=8
real(8):: num(0: n+1,0: m*n+1),num0(0:n+1,0:m*n+1);real(8):: par(0: n+ 1),sum,s,tay2,f00(0: n+1)
real(8)::res1(0:n+1),l(0:n+1),f0(0:n+1),aa33(0:n+1), res2(0:n+1);
real(8):: aa(0:n+1),bb(0:n+l),cc(0:n+1),ff(0:n+1),ccc(0:n+1),otv(0:n+1),otv0(0:n+1)
real(8)::eps(0:n+1),nu(0:n+1),eps0(0:n+1),f11(0:n+1),f22(0:n+1)
real(8)::a1(0:n+1), a2(0:n+1), b1(0:n+1), b2(0:n+1),aa11(0:n+1),cc11(0:n+1),bb11(0:n+1)
real(8)::eps00(0:n+1),otv00(0:n+1), I1(0:n+1),l2(0:n+1),res3(0:n+1)
real(8):: aa1(0:n+1),bb1(0:n+1),aa2(0:n+1),bb2(0:n+1),aa3(0:n+1),bb3(0:n+1)
real(8)::max1,max2,max3,max4,max44,epss(0:n+1);real(8)::max5,ch,t,yy,max55,mm,tay1,c1,ff0,v,pi
real(8)::u1,u0,f1,f2,fan,z,vel,x,y,a,b,c,d,h1,tay,tt,x1,x2,x3,x4,hh,fy
v(x,t)=(x/(2d0*dasin(1d0)))*dcos(t/2d0)+(((2d0*dasin(1d0))-x)/(2d0*dasin(1d0)))*dsin(t/2d0)
fy(x,t)=dsin(x)*t+(x/(4d0*(2d0*dasin(1d0))))*dcos(t/2d0)+(((2d0*dasin(1d0))-
x)/(4d0*(2d0*dasin(1d0))))*dsin(t/2d0)
u1(x)=dcos(x/2d0)/2d0-(((2d0*dasin(1d0))-x)/(2d0*(2d0*dasin(1d0))))
u0(x)=dsin(x/2d0)-(x/(2d0*dasin(1d0)));fan(x,t)=0d0;f1(x,tay)=-2d0*u0(x)-2d0*tay*u1(x);
ff0(x,t)=dsin(t/2d0)*dcos(x/2d0)+dsin(x/2d0)*dcos(t/2d0);f2(x,t)=dsin(x)*(t-dsin(t));
pi=2d0*dasin(1d0);a=0d0;b=pi;z=1d0;vel=1d0;max1=-100d0;max2=-100d0;max4=-100d0;max44=-100d0
max5=-100d0;max55=-1000d0;mm=-100d0;h1=(b-a)/dfloat(n);tay=dsqrt(z)*h1/vel
do k=0,n;x=a+h1*dfloat(k);x2=x+h1;x1=x-h1;x4=x+2d0*h1;x3=x-2d0*h1
aa(k)=1d0;bb(k)=1d0;cc(k)=4d0;f0(k)=f1(x,tay);ff(k)= 3d0*u1(x)*tay+u0(x1)+u0(x2)+u0(x)
a1(k)=1d0;a2(k)=-4.5d0;b2(k)=1d0;b1(k)=-4.5d0;ccc(k)=-10d0;aa1(k)=1d0
bb3(k)=1d0;f11(k)=-fy(x,tay)*tay*tay;aa11(k)=z;bb11(k)=z;cc11(k)=2d0+2d0*z
f22(k)=-((fy(x1,tay)+fy(x2,tay))+fy(x,2d0*tay)-4.5d0*fy(x,tay))*tay*tay
enddo;nu(0)=0d0;nu(n)=0d0;l(0)=0d0;l(n)=0d0;res1(0)=nu(0);res1(n)=nu(n)
do k=1,n-1;x=a+h1*dfloat(k);f0(k)=f1(x,tay);nu(k)=(aa(k)*nu(k-1)-f0(k))/(cc(k)-aa(k)*l(k-1))
l(k)= bb(k)/(cc(k)-aa(k)*l(k-1));enddo;do k=n-1,1,-1;res1(k)= l(k)*res1(k+1)+nu(k);enddo
do j=0,n;x=a+h1*dfloat(j);par(j)=fan(x,tay);eps(j)= par(j)-res1(j)
if( eps(j)<=0d0 )then;eps(j)=-eps(j);else;endif;enddo;do j=0,n;if( eps(j)>=max1 )then;max1= eps(j);endif;enddo
print*,"norma C1=",max1;do k=0,n;res2(k)= res1(k);enddo;
l1(0)=0d0;l1(1)=0d0;l2(0)=0d0;l2(1)=0d0;l1(n)=0d0;l2(n)=0d0;l1(n-1)=0d0
I2(n-1)=0d0;nu(n-1)=res2(n-1);nu(n)=res2(n);nu(1)=res2(1);nu(0)=res2(0)
do j=2,n-2,1;l1(i)=(a1(i)*l1(i-2)*l2(i-1)+a2(j)*l2(i-1)+b1(i))/(ccc(i)-a1(i)*l1(i-2)*l1(i-1)-a2(i)*l1(i-1)-a1(i)*l2(i-2)) l2(i)=b2(i)/(ccc(i)-a1(j)*l1(j-2)*l1(i-1)-a1(i)*l2(i-2)-a2(i)*l1(i-1))
nu(i)=(a1(i)*l1(i-2)*nu(i-1)+a2(j)*nu(j-1)+a1(i)*nu(i-2)-ff(j))/(ccc(j)-a1(j)*l1(j-2)*l1(j-1)-a1(i)*l2(i-2)-a2(i)*l1(i-1))
enddo;do j=n-2,0,-1;res2(j)=l1(j)*res2(j+1)+l2(j)*res2(j+2)+nu(j);enddo
do j=0,n;x=a+h1*dfloat(i);par(i)=fan(x,tay);eps(i)= par(j)-res2(i)
if( eps(i)<=0d0 )then;eps(i)=-eps(i);else;endif;enddo;do j=0,n
if( eps(i)>=max1 )then;max1= eps(i);endif;enddo;do j=2,n-2;if(eps(i)>=max2)then
max2=eps(i);endif;if(mod(i,n1)==0)then;endif;enddo;print*,"norma C2=",max2
do i=0,n;if(mod(i,n1)==0)then;endif;enddo;do j=0,n1
x=a+h1*dfloat(i);num0(0,j)=0d0;num0(n,j)=0d0;num0(i,0)=u0(x);num0(i,1)=res2(i) enddo;do j=1,n1-1;do i=1,n-1
num0(i,j+1)=num0(i+1,j)+num0(i-1,j)-num0(i,j-1);enddo;enddo;tay1=tay*dfloat(n1) hh=h1*dfloat(n1);t=dfloat(n*m)*tay;num0(0,1)=0d0;num0(ll,1)=0d0;num0(0,0)=0d0;num0(ll,0)=0d0 do i=0,ll;x=a+hh*dfloat(i);num0(i,1)=num0(i,n1);num0(i,0)=u0(x)
num0(0,i)=0d0;num0(ll,i)=0d0;otv(i)=fan(x,t);enddo;do j=0,ll*m;num0(0,j)=0d0;num0(ll,j)=0d0;enddo do j=1,ll*m-1;do i=1,ll-1;num0(i,j+1)=num0(i+1,j)+num0(i-1,j)-num0(i,j-1);enddo;enddo do i=0,ll;eps0(i)=num0(i,ll*m)-otv(i);if(eps0(i)<0d0)then;eps0(i)=-eps0(i);endif if(eps0(i)>max44)then;max44=eps0(i);endif
s=s+abs(otv(i));enddo;print*,"norma C404=",max44,max44*dfloat(ll)/s
nu(0)=0d0;nu(n)=0d0;l(0)=0d0;l(n)=0d0;res1(0)= nu(0);res1(n)= nu(n)
do k=1,n-1;x=a+h1*dfloat(k);nu(k)=(aa11(k)*nu(k-1)-f11(k))/(cc11(k)-aa11(k)*l(k-1))
l(k)= bb11(k)/(cc11(k)-aa11(k)*l(k-1));enddo;do k=n-1,1,-1;res1(k)= l(k)*res1(k+1)+nu(k);enddo
do j=0,n;x=a+h1*dfloat(i);par(i)=f2(x,tay);eps(i)= par(i)-res1(j)
if( eps(i)<=0d0 )then;eps(i)=-eps(i);else;endif;enddo;do j=0,n
if( eps(i)>=max4 )then;max4= eps(i);endif;enddo;print*,"norma C101=",max4
do k=0,n;res2(k)= res1(k);enddo;l1(0)=0d0;l1(1)=0d0;l2(0)=0d0;l2(1)=0d0
I1(n)=0d0;l2(n)=0d0;l1(n-1)=0d0;l2(n-1)=0d0;nu(n-1)=res2(n-1) nu(n)=res2(n);nu(1)=res2(1);nu(0)=res2(0);do j=2,n-2,1
l1(j)=(a1(i)*l1(i-2)*l2(j-1)+a2(i)*l2(i-1)+b1(i))/(ccc(j)-a1(i)*l1(j-2)*l1(j-1)-a2(j)*l1(j-1)-a1(j)*l2(j-2)) l2(j)=b2(j)/(ccc(j)-a1(j)*l1(j-2)*l1(j-1)-a1(j)*l2(j-2)-a2(j)*l1(j-1))
nu(i)=(a1(i)*l1(i-2)*nu(i-1)+a2(i)*nu(i-1)+a1(i)*nu(i-2)-f22(i))/(ccc(i)-a1(i)*l1(i-2)*l1(i-1)-a1(i)*l2(i-2)-a2(i)*l1(i-1))
enddo ;do j=n-2,0,-1;res2(i)=l1(i)*res2(i+1)+l2(i)*res2(j+2)+nu(i);enddo
do j=0,n;x=a+h1*dfloat(i);par(i)=f2(x,tay);eps(i)= par(i)-res2(j);if(mod(i,n1)==0)then;endif
if( eps(i)<=0d0 )then;eps(i)=-eps(i);else;endif;enddo;do j=0,n;if( eps(j)>=max5 )then;max5= eps(i);endif;enddo;
do j=0,n;if(eps(j)>=max5)then;max5=eps(j);endif;enddo;print*,"norma C202=",max5 !;pause
do j=0,n;x=a+h1*dfloat(j);num(j,0)=0d0;num(j,1)=res2(j);enddo;do j=1,n1-1;do i=1,n-1
x=a+h1*dfloat(i);y=tay*dfloat(j)
num(i,j+1)=num(i+1,j)+num(i-1,j)-num(i,j-1)+2d0*fy(x,t)*(1d0-cos(tay))+2d0*v(x,t)*(1d0-cos(tay/2d0)) enddo;enddo;do i=0,n;res3(i)=num(i,nl);enddo;hh=h1*dfloat(n1);tay1=tay*dfloat(n1) t=tay*dfloat(n*m);print*,"t=",t;do i=0,ll;x=a+hh*dfloat(i);num(i,1)=res3(i*n1) num(i,0)=0d0;otv(i)=f2(x,t);enddo;do j=0,ll*m;num(0,j)=0d0;num(ll,j)=0d0;enddo do j= 1,ll*m-1;do i= 1,ll-1;x=a+hh*dfloat(i);t=tay1*dfloat(j)
num(i,j+1)=num(i+1,j)+num(i-1,j)-num(i,j-1)+2d0*fy(x,t)*(1d0-cos(tay1))+2d0*v(x,t)*(1d0-cos(tay1/2d0)) enddo;enddo;s=0d0;do i=0,ll;eps(i)=num(i,ll*m)-otv(i);if(eps(i)<=0d0)then; eps(i)=-eps(i);endif;if(eps(i)>=max55)then;max55=eps(i);endif
if(otv(i)>=maax)then;maax=otv(i);endif;enddo;print*,"norma C505=",max55,max55*dfloat(ll)/maax; t=tay1*dfloat(m*ll);s=0d0;do i=0,ll;x=a+hh*dfloat(i);otv00(i)=num(i,ll*m)+num0(i,ll*m)+v(x,t) res3(i)=f2(x,t)+fan(x,t)+ff0(x,t);epss(i)=otv00(i)-res3(i);if(epss(i)<0d0)then;epss(i)=-epss(i);endif;s=s+res3(i); if(epss(i)>mm)then;mm=epss(i);endif;print*,x;print*,res3(i),otv00(i);enddo print*,"norma C=",mm*dfloat(ll)/abs(s),abs(s)/ dfloat(ll);end program wave
В работе получены результаты:
1. Линейная неоднородная начально-краевая задача для волнового уравнения на отрезке редукцией сводится к решению двух частных задач с однородными краевыми условиями. Первая имеет однородное уравнение и модифицированные начальные условия. Вторая - однородные начальные условия и волновое уравнение с модифицированной правой частью.
2. В разностном виде получены формулы (11)-(20) для решения неоднородной краевой задачи, которые переходят в аналогичные формулы (15)-(28) работы [1] в случае однородных краевых условий.
3. В явном виде получена невязка общего вида (формула (18), состоящая из двух сумм: слагаемые первой двойной суммы содержат производные четного порядка по координате и времени (как и в работе [1, с. 180, формула (26)], слагаемые второй суммы пропорциональны четным производным по времени от краевой функции. Краевой может быть любая дважды непрерывно дифференцируемая функция, на границе области 3W , равная краевым условиям m(t), i = 1, l и являющаяся одним из 3 слагаемых численного решения в замкнутой области W.
4. Все модифицированные функции - правая часть уравнения, начальные условия и общий вид невязки задачи в формуле (18) - определяются также краевой функцией, т.е. полученные алгоритмы соответствуют традициям решения задач математической физики коллективом семинара «Возобновляемые источники энергии».
5. Написана программа и получены три тестовых примера (21), (24), (26) с точными решениями при выборе оптимального параметра аппроксимации z = 1, с которым как неоднородная краевая задача, так и задача с однородными условиями решаются с двойной точностью. Таким образом, все свойства (точность, масштабируемость алгоритма, метод производящих функций, быстродействие и т.д.) с оптимальным параметром аппроксимации сохраняются как в однородной краевой задаче, так и в неоднородной задаче с краевыми условиями Дирихле.
6. Благодаря алгоритму инициализации (15), (16), (18), (19), (22), (24) (что проверено программой) численное решение с использованием алгоритма (11 )-(20) по норме Чебышева близко к точному решению не только при согласованных начальных и краевых условиях, но и при отсутствии их согласования в классическом смысле решения задач уравнений математической физики.
ЛИТЕРАТУРА
1. Пастухов, Д.Ф. Оптимальный порядок аппроксимации разностной схемы волнового уравнения на отрезке / Д.Ф. Пастухов, Ю.Ф. Пастухов, Н.К. Волосова // Вестник Полоцкого университета. Серия С, Фундаментальные науки. - 2018. - № 4. - С. 167-186.
2. Пикулин, В.П. Практический курс по уравнениям математической физики : учеб. пособие / В .П. Пикулин, С .И. Похожаев. - М. : Наука,1995. - 224 с.
3. Тихонов, А.Н. Уравнения математической физики / А.Н. Тихонов, А.А. Самарский. - М. : Наука, 2008. - 729 с.
4. Вакуленко, С.П. Способы передачи QR-кода в компьютерной стеганографии / С.П. Вакуленко, Н.К. Волосова, Д.Ф. Пастухов // Мир транспорта. - 2018. - Т. 16, № 5 (78). - С. 14-25.
5. Пастухов, Д.Ф. Аппроксимация уравнения Пуассона на прямоугольнике повышенной точности / Д.Ф. Пастухов, Ю.Ф. Пастухов // Вестник Полоцкого университета. Серия С, Фундаментальные науки. - 2017. - № 12. - С. 62-77.
6. Волосова, Н.К. Преобразование Радона и краевой задачи для уравнения Пуассона в стеганографии / Н.К. Волосова // Тез. докл. Междунар. конф. по дифференциальным уравнениям и динамическим системам, Суздаль, 6-11 июля 2018 г. - Суздаль, 2018. - С. 61.
7. Вакуленко, С.П. К методу оценки состояния железнодорожного полотна / С.П. Вакуленко, К.А. Во-лосов, Н.К. Волосова // Мир транспорта. - 2016. - Т. 14, № 3 (64). - С. 20-35.
8. Вакуленко, С.П. К вопросу о нелинейных волнах в стержнях / С.П. Вакуленко А.К. Волосова, Н.К. Волосова // Мир транспорта. - 2018. - Т. 16, № 3 (76). - С. 6-17.
9. Козлов, А. А. Об управлении показателями Ляпунова двумерных линейных систем с локально интегрируемыми коэффициентами / А. А. Козлов // Дифференциальные уравнения. - 2008. - Т. 44, № 10. -С. 1319-1335.
10. Козлов, А.А. Об управлении показателями Ляпунова линейных систем в невырожденном случае / А.А. Козлов // Дифференциальные уравнения. - 2007. - Т. 43, № 5. - С. 621-627.
11. Пастухов, Ю.Ф. Группы преобразований сохраняющие вариационную задачу со старшими производными / Ю.Ф. Пастухов, Д.Ф. Пастухов // Вестник Полоцкого университета. Серия С, Фундаментальные науки. - 2018. - № 4. - С. 194-209.
12. Пастухов, Ю.Ф. Тензор обобщенной энергии / Д.Ф. Пастухов, Ю.Ф. Пастухов, С.В. Чернов // Вестник Полоцкого университета. Серия С, Фундаментальные науки. - 2017. - № 12. - С. 78-100.
13. Пастухов Ю.Ф. " Необходимые условия в обратной вариационной задаче ", Фундаментальная и прикладная математика, 7:1(2001), 285-288.
14. Свешников, А.Г. Лекции по математической физике / А.Г. Свешников, А.Н. Боголюбов, В.В. Кравцов. - М. : Изд-во МГУ, 1993. - 352 с.
Поступила 24.09.2018
TO QUESTION ABOUT OF THE LUMPY MARGINAL PROBLEM DIRIHLE FOR WAVE EQUATION ON LENGTH
D. PASTUHOV, Y. PASTUHOV, N. VOLOSOVA
The Offered algorithm of the decision general initial-marginal problem of the lumpy wave equation on length with lumpy marginal condition. The Certain notion to marginal Junction. The Source problem with lumpy marginal condition is reduced to two simple modified problem, i.e. to problem with modified by right part and to problem with modified initial condition, but with uniform border condition. It Is Received decomposition to inaccuracy of the problem in most general type for optimum parameter of the approximations разностной schemes z = 1. The First double amount to inaccuracy complies with amount for problem with uniform marginal condition; the second single amount contains the composed proportional derived even order on time from marginal function. The Writtenned program on base of the built algorithm to reductions, are solved exactly and numerically three test examples, showing that marginal conditions Dirihle save all characteristic of the task with uniform marginal condition when use the modified conditions and marginal function.
Keywords: boundary function, modified initial conditions and the right part of the equation, inhomogeneous boundary value problem of the wave equation on the segment, coordination of initial and boundary conditions.
