ОПТИМАЛЬНЫЙ ПАРАМЕТР АППРОКСИМАЦИИ РАЗНОСТНОЙ СХЕМЫ ВОЛНОВОГО УРАВНЕНИЯ НА ОТРЕЗКЕ Текст научной статьи по специальности «Математика»

МАТЕМАТИКА

УДК 517.6: 517.958

ОПТИМАЛЬНЫЙ ПАРАМЕТР АППРОКСИМАЦИИ РАЗНОСТНОЙ СХЕМЫ ВОЛНОВОГО УРАВНЕНИЯ НА ОТРЕЗКЕ

канд. физ-мат. наук, доц. Д.Ф. ПАСТУХОВ, канд. физ-мат. наук, доц. Ю.Ф. ПАСТУХОВ (Полоцкий государственный университет);

Н.К. ВОЛОСОВА

(Московский государственный технический университет им. Н.Э. Баумана)

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

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

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

Напомним, что порядком аппроксимации дифференциального оператора разностным оператором [4, с. 102] называется максимальное положительное число p, если существуют положительные числа p,С > 0, не зависящие от шага сетки h такие, что норма невязки (разности дифференциального и разностного оператора) не превышает ||( Lu — LhUh < Chp.

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

Нами построен алгоритм инициализации задачи, т.е. аппроксимация второго временного слоя решения по начальным данным задачи. Алгоритм сводится к прогонке либо трехдиагональной системы линейных алгебраических уравнений (СЛАУ), либо к прогонке сначала трехдиагональной, а затем пятидиагональной СЛАУ. Найдены также достаточные условия корректности формул прогонки, более слабые, чем диагональное преобладание элементов матрицы. Алгоритм инициализации дает приближение

второго слоя решения по начальным условиям с относительной погрешностью не хуже чем 10—10 при числе узлов 500. Инженерный американский продукт ANSYS Fluent завершает решение задач с относи-

—3

тельной точностью 10 , и с более грубой точностью на процессе инициализации задачи, т.е. с точно-—1 —2

стью 10 —10 . Благодаря применению спектрально устойчивых разностных схем относительная погрешность начальных данных уменьшается от значения 10—8 —10—10 до величины 10—15. В работе нами

построен алгоритм укрупнения ячеек сетки (масштабирования) с коэффициентом масштабирования

2

l > 1, позволяющий сократить число вычислений в сотни раз (l раз). Неоднородная начально-краевая задача при числе узлов 300 на редкой координатно-временной сетке решается с помощью указанных

здесь методов с двойной точностью 10—15 (т.е. первые 15 десятичных знаков аналитического и численного решения на последнем временном слое во всех координатных узлах совпадают).

Некоторые двухмерные и трехмерные колебания в приложениях часто сводятся к одномерным колебаниям [5, с. 181], например, в прямоугольном кристалле колебания атомных параллельных плоскостей порождают одномерные звуковые волны, поэтому решаемая нами численно задача полезна и в многомерных случаях.

Постановка задачи. Рассмотрим начально-краевую задачу для неоднородного волнового уравнения на отрезке [а, Ь] .

'2 2

= a2 Э^ + a2 >0, (л,г)е (а,Ь)х(0,Т),

Эг2 Эл2

и(л,0) = ф(л), л е [а, Ь], (1)

щ (л, 0) = у(л), л е [а, Ь],

и(а,г) = и(Ь, г) = 0, ле [а, Ь], где и(л, г) - точное аналитическое решение задачи (1); f (л, г) - функции внешнего источника, ф(л) -начальное смещение и у(л) - начальная скорость точек струны с закрепленными концами на отрезке "(л, г) е [а, Ь] х [0, Т].

Аналитической задаче (1) сопоставим разностную задачу на равномерной сетке, используя пятиточечный шаблон - крест:

и п+1 + и п—1 _ 2ип ип + щп _ 2щп

ит + ит 2ит _ „2 ит+1 + ит—1 2ит

+ f ( лт, гп ),

(2)

= а

х2 к2

лт = а + тк, гп = пх, [т,п] е [1,М — 1]х[1, N — 1], и(тк,0) = ф(тк), те [0,М],

и(тк,х) = ф^тк), те [0, М], л) = ¥ (ф^л),у(л)), и(0, пх) = и(N, пх) = 0, п е [0, N].

, Ь — а Т „ 2 л, "

где к =-, х = — - временной и координатный шаги сетки; а - квадрат фазовой скорости волны.

N М

Задание начальных условий ф(л), у(л) в линейной аналитической задаче (1) эквивалентно заданию значений двух начальных временных слоев решения и(тк,0) = ф(тк) = и°, и(тк,х) = ф^(тк) = и1т, т е [0, М] в разностной задаче (2), поскольку их запрашивает узловое значение ит в первой рекуррентной формуле системы (2). Зависимость второго временного слоя ф^л) = ¥(ф^л),у(л)) в системе уравнений (2) как функции начальных условий определяется с помощью предварительного численного алгоритма инициализации задачи методом прогонки.

а2х2

Обозначим параметр 7 = —— и перепишем разностное уравнение (2) в эквивалентном виде

и п+1 + и п—1 — 21п = ит + ит 2ит =

к2

2 (ит+1 + ит—1 — 2ит ) + f(лm, гп )х = 2 (ит+1 + ит—1 — 2ит ) + /т,пх . (3)

Потребуем, чтобы разностное уравнение (3) аппроксимировало первое уравнение системы (1) с максимальным алгебраическим порядком, далее разложим узловые значения Цт^1, ит—1, ит—1, ит+1 в ряд Тейлора в центральном узле (т, п) «(лт, гп) с вектором шага соответственно (к, х):

п ¥ / V п (я=2к) ¥ о -\2к п

1 Э ит х 1) Э ит '^ 2 Э ит

тх = ^^^ ~ =

п+1 + п—1 — 2 п = —2 п + 2 п +V"1 1 Э ит хя + V-1 ( 1) Э ит „.я 4 _ ' ^ 2 Э ит _2к ит + ит 2ит = 2ит + 2ит + ^ / ^ я х +2-1

я=1 (я)! Эгя (я)! Эгя

к=1 (2к)! Эг2к

= X'

;>2, п

2 о ит

+ I

о -\2к п 2 0 ит х2к

а»2 к=2 (2к)! а»2к

Меняя в последней формуле буквенные обозначения (т « п), (х «»)(А « х), получим

индексов

переменных

/ п + п _ 2тп\ + г х2 = 2 ^ит+1 + ит_1 2и^ + /т,пх = 2

^ д2, п

2 0 ит

+ I

о -\2к п 2 0 ит ,2к

дх2 к=2(2к)! дх2к

+ ^ Х2

1

Приравняем правые части последних двух выражений, используя тождество

2 п

д2и:

д» 2

2

а х

= а

2 д2и

2п

т + /" 1 ¿т

А

2

дх2

д2 2 д и

имеем х

2п

2 д ит

+ I

а»2 к=2(2к)! а»2к

2к п 2 п

2 д ит х2к __2 2 д ит

х =х а о

дх2

оо о ^2к. п

+ х2 +v 2 д и + /т,пХ + ¿2(2к)! д»2к

т х2к =

■+ I

2к п

2 д ит 2к

дх2 (2к)! дх2к

+ /т,пх ^

« I

2к п о

2 д ит х2к =2 I

к=2 (2к)! д»2к

2к п 2 д ит_А2к

к=2(2к)! дх2к

(4)

Невязка первого уравнения системы (2) равна невязке уравнения (4), т.е. разности левой и правой частей уравнения (4):

о2

^) = ¿2 М!

(а 2к

2к п

д2кип

д ит х2к _ 2 " -т а

2к п

д ит 2к

д»

2к

дх

2к

(5)

J

2 2 2 2

г, , ~ а X . д и 2 д и

Замечание 1. Если параметр 2 =—— = 1 и волновое уравнение —— = а —— однородное, то

А2 д»2 дх2

2п

из

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

2 ди

т = а2

д2ип и

2к п 2к п

д ит 2к д ит

д»

2к

=а

дх

2к

в

2к п 2к п

д ит 2к д ит 2к

-х _ 2——— А = А

тогда

а2к 2к д и

д»

2к

дх

2к

дх

2к

д»2 дх2

формуле (5) преобразуем выражение

2к 2к \ д2к п

X а ,2к д и,

следует тождество в скобках

А

2к

2к д ит к

--2 = А —— (2 _2) = 0. Другими словами, для пара) дх 1 '

метра 2 = 1 однородное одномерное волновое уравнение имеет бесконечный порядок аппроксимации! В этом случае разностная схема (3) с /тп °0,2 = 1 точна для многочленов произвольной степени. Равномерная норма погрешности аппроксимации первого уравнения задачи (1) первым уравнением разностной задачи (2) не зависит от шага сетки А (аппроксимация происходит с двойной точностью даже для крупного шага сетки).

Для упрощения формулы (5) докажем первое утверждение.

2 п 2 п

о ит 2 о ит

Утверждение 1. Для неоднородного одномерного волнового уравнения —= а -т + }т

д»

справедлива формула

д2 р

п д 2 Р--п

д» 2 Р

ит = а2Р

Р_1

дх2 Р

т + I а I=0

21

д2(Р_1) Г

т, п

дх21 д»2( Р_1 _1)

, Р > 2.

дх

2

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

а4, п

д ит д»4

э_ д»2

2 ( Г|2, п

а2 ^+т

дх

2

=а

д! дх2

2п

д 2и

дх2

+ 1т

д2 Г д !т,п

4п 4 д ит 2

+- = а4-т + а2

д 1т,п + д 1т,п

д» 2

дх4

дх

2

д»

2

что верно. Следовательно, получаем формулу (6), если р = 2 .

и

а

Пусть справедлива также формула

э2( Р 1)и1 = а2( р—1)

э2(р—V 2, э2(р—2)^

"т + V"1 р—1)

, тогда

э2рит Э2

Эг2 р =Э?

Эг2(Р—1 Эл2(Р—1) £0 Эл21Эг2(Р—1—2)

Г Э2(Р—1)ит . Р^2 2, Э2(Р—2)fтп 1 2(р—1) Э2(Р—1) Э2ит . Р—2 2, Э2(р—1Ь

,2( р—1):

т + ^ „2,

Эл2(Р—1) £0 Эл21Эг2(Р—

■2)

= а 2( Р—1)

__ит + 'V а2,_

Эл2(Р—1) Эг2 ¿0 Эл21Э12(Р—1—1

= а2(р—1) Э2(р—1) Г 2 Э2

Эл2(Р—1)

2 Э ит + Ь

2 ./т,п

Эл2

р—2 Э2(Р—11) Ь Э2Р„п Э2(Р—11) Ь

21 Э /т,п = а2 р Э ит + а2( р—1) Э_Ьт

+ V

а

1=0

Эл2,Эг2(Р—,—1) Эл2Р

Эл2( Р—1)

р—2 Э2(р—1) г

V- 21 Э Л

+ V

-,2 р п р—1 Э2( Р—1) Ь т,п = а2р Э ит , ^ „21 ->т

+

,=0 Эл2,Эг2(Р—,—1) Эл2Р ,=0 Эл2,Эг2(Р—,—1)

т.е. формула (6) (утверждение 1) доказана для произвольного целого р > 2 .

Используя формулу (6), преобразуем невязку разностного уравнения (3), подставив ее в формулу (5):

д(«т)=V

к=2

(2к)!

Га2к, и л2к, и 1 ¥ о Г Г д2к, п к—1 Э2(к—1) Ь 1 з2к, п

Э Э . 2 2к 2к Э Ит * 1 21 Э /, -1 "

х а тт + / а

2к

х2к итк2к

Эг

2к

-Т — 2-

Эл

2к

2

У

£ (2к)!

¥

= V — к=2 (2к)!

¥

= V — к=2 (2к)!

ГГ 2к

х2ка2к

2к

2

V V

Э2к, п к—1 Э2(к—1) Ь Л

Э ит + х2к 21 Э Л 2к +х

+ V а21_

Эл2к Й Эл21 Эг2(к ~1—1)

Э ит 2к

V V

У

Эл2

,=0

Эл2,Эг2(к—,—1)

г

д2к, п к —1 Э2(к —1) Ь

Э ит , „2кк—1 21 Э Н

к2к (2к — 2)^ит + х2к V

,=0

(2=1) ™ 2 = У—=-к=2 (2к)!

Эл21 Эг2(к—1—1 Г к —1 Э2(к—1) г 1

^ 21 " ¿т,п

Й а Эл21 Эг2(к—1—1,.

(2=1)

Эл

2к

(7)

Выпишем несколько первых слагаемых для невязки Я(ит) из формулы (7):

2х

к(ит)

4 Г э2 Н Э2 Н 1

2 Э ¿т,п Э ¿т,п

Эл2

Эг2

2х_ 6!

6 Г э4 f Э4 f Э4 f 1

4 Э /т,п 2 Э ¿т,п Э ¿т,п

+а

Эл4 Эл2Эг2 Эг4

2х'

" 8!

8 Г Э6 Ь Э6 Ь Э6 Ь Э6 Ь 1

6 Jm,n 4 Jm,n 2 " ^т,п " ./т,п а --—+ а —;—— + а —-—Г +

Эл6 Эл4Эг2 Эл2Эг4 Эг6

2х

10

Г8

10!

81

8 Э ¿т,п 6 Э Ьт,п 4 Э Ьт,п 2 Э /т,п Э /т,п

а -^—+ а —-—— + а —;—- + а —-—- + -

Эл8 Эл6Эг2 Эл4Эг4 Эл2Эг6 Эг8

2т

/

.12 Г Э10 Ь Э10 Н Э10 Ь Э10 Ь Э10 Ь Э10 Ь 1

10 Э Ьт,п 8 Э Ьт,п 6 Э Ьт,п 4 Э Ьт,п 2 Э !т,п Э !т,п

а -ГТГ^ + а8—^—+ а6——+ а —;—^ + а —+

12!

Эл

10

Эл8Эг2 Эл6Эг4 Эл4Эг6 Эл2Эг8 Эг

10

(8)

Определение 1. Одномерное уравнение в частных производных

2 Э2 Ь Э2 Ь _ 2 _ а2—+ —= 0, а2 > 0 Эл2 Эг2

(9)

Э2Ь 2 Э2Ь п 2 п

по аналогии с волновым уравнением —— — а —— = 0, а > 0 назовем волновым уравнением комплекс-

Эг

Эл

2

ного аргумента. Рассмотрим свойства решений (9).

+

2

+

+

+

+

+

Утверждение 2 (свойства волнового уравнения комплексного аргумента).

1) Уравнения xi + at = Q = const, xi — at = C2 = const являются уравнениями характеристик для (9). Действительно,

2 Э2f Э2f - ( э .э Y э .э . . _

a —— +—— = 0 ^ I a — + i— II a--i— I f = 0 ^

Эх2 Эt2 l Эх ЭtJl Эх Эt 1 f

3f .

a — + i — = 0 Эх Эt

3f -

a—— i — = 0 Эх Эt

dx = dt ai

dx = dt . a —i

xi — at = Q ^ х + iat = — iC\ = Q xi + at = C2 ^ х — iat = — C = C2

2) Решением волнового уравнения комплексного аргумента является линейная комбинация произвольных дважды дифференцируемых функций на комплексных характеристиках

f (х, t) = Qu(x + iat) + C2v(х — iat), u(w), v(w) e C2 (Z).

(10)

Действительно,

2 ^ + = a2 Эх2 Эt2

(Qu (х + iat) + C2V (х — iat)) + (ia)2Qu (x + iat) + (—ia)2C2V (х — iat) = 0.

3) Пусть функция / (х± ш») = /1 (х,») + г/2 (х,») комплексного аргумента (х ± га») — решение уравнения (9), тогда функции действительного аргумента /[(х, г), /2(х,»), т.е. действительная и мнимая части / (х± Ш), также являются решениями уравнения (9). Действительно, имеем

2 Э2f (х± iat) Э2f (х±iat) 2 Э2 (f1(х,t) + if2(x,t)) Э2 (fx(x,t) + гf2(x,t)) a ----1----= 0 = a ----1----= 0 ^

Эх2

Эt2

Эх2

Эt2

2 Э2 ( f (x, t) + f2(x, t) ) Э2 ( f (x, t) + if2 (x, t) ) a ----1----= 0 +1 • 0 ^

Эх2

Эt2

a2 Э2f1(x, t) + Э2fj(x, t) = 0

Эх2

Эt2

a2 Э2 f2 (x, t) + Э2 f2( x, t) = 0

Эх2

Эt2

Таким образом, показано, что функции /[(х,»), /2(х,») — решения уравнения (9). Последняя система уравнений следует из определения равенства двух комплексных чисел. Утверждение 2 доказано.

Приведем примеры из класса решений волнового уравнения комплексного аргумента, используя утверждение 2.

2 2 2 2

Пример 1. /[(х,») = Яе(х + га») = х _ а » . Действительно,

-42/2 2.2\ -ч2/ 2 2.2 \

Э (х — a t ) Э (х — at )

Эх 2

Эt2

= 2a2 — 2a2 = 0.

4 3 3

Пример 2. f2(х, t) = Im (х + iat) = 4х at — 4xa t. Действительно,

Э2 (4x3at — 4xa3t3) Э2 (4x3at — 4xa3t3) —^---'- + —^---'- = 24a3 xt — 24a3 xt = 0.

Эх 2

Эt2

a

a

a

Пример 3. /э(x, t) = Im (exp(x + iat)) = exp (x) sin(at). Действительно,

2 d2 (exp(x)sin(at)) Э2 (exp(x)sin(at)) ¡ 2 2 a —----- +------- = ( a - a

dx2

dt2

(a2 - a2) exp (x) sin (at) = 0.

Замечание 2. Если правая часть Ь(л, г) неоднородного волнового уравнения действительного аргумента задачи (1) удовлетворяет однородному волновому уравнению мнимого аргумента, то порядок невязки разностной схемы задачи (2) увеличивается со второго порядка до четвертого. Более того, невязку в этом случае удается выразить через частные производные по одной независимой переменной г. Нужно учесть, что невязка уравнения (3) на 2 порядка больше невязки уравнения (2). Действительно, из формулы (8) получим

R(u") = 2t6 ( Э4/

R(um ) = "6Т

п2

m,n 2 d

+ a —у

dt4 dx2

( n2 .

2 d fm,n d

a ----+

2, AA

dx2

at2

2t' " 8!

8 ( ^4 ( а2 , A

dt

V v

2 d fm,n d fm

a -г--

dx2

dt2

+a

л4 ( a2 f d2 / A A

2 d /m.n d Jr.

dx4

dx2 dt

6 ( 2

+

2t

10 (d8 /

10!

dt8

m,n 2 d

+ a

6 ( d2 / d2 / 2 /m, n /m

dx2dt4

2A

m,n

V dx2 dt

\

2t6 d4/m,n 2t10 d8 f

+a

d

2

dx

a

2 d2 /m,n d2 /m

dx

2

dt2

+... =

6! dt4 10! dt8

+... = 21

t4k+2 d4k/m,

(4k + 2)! dt4k

(11)

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

разностное уравнение (3) на спектральную устойчивость. Учтем, что численное решение ит (л, г) удовлетворяет также уравнению (3), как и проекция аналитического решения на узлы сетки (и(л, г) )т = ит, ет = Ит — ит Ц<т+1 + Ц^т 1 — 2ит = 2 (Ц<т+1 + Ц^т 1 — 2ит) + Ьт пх2. Вычитая из последнего

выражения уравнение (3) получим однородное уравнение относительно невязок в каждом внутреннем узле сетки (т, п):

en+1 +pn-1

mm

en+1 +pn-1

mm

Подставляя в последнее уравнение ошибку округления вида ет =1п(ф)егтф, получим спектраль-

ное уравнение

1 +1/1-2 = z(eij + e-ij - 2) = -4zsin2 (ф| l2 - 2- 2zsin2 (1 +1 = 0 .

(12)

Как следует из определения спектральной устойчивости [4, с. 125],

® 0 Ы1(ф)| < 1 "фе [0,2я].

Используя коэффициенты квадратного уравнения (12), получим =

=1/1=1.

Если корни действительные, то при < 1 ^ > 1. Тогда возможен единственный результат

|1l| = = 1 [0,2р] , что заведомо невозможно, т.к. = -2zsin2 ^ jjj+ ^-2z зависит от j. Следовательно, имеем пару комплексно сопряженных корней:

1 - 2z sin21 ф|| -1

+

+

+

2

(

m

n

u

m

• 2 | j

1U =|1 -2zsin21 ^ I I±üИ-I 1 -2zsin21 ^

• 2 | j

X12r =|1 - 2z sin21j II +1 -|1 - 2 z sin2 I^II = 1, Vje[0,2p].

• 2 ( j

Выясним условия, при которых дискриминант отрицательный. 1 -2zsin2(j||2 -1 <0«1 -4zsin2(jj + 4z2sin4(jj< 1«4sin2(jjz(z-1)<0«

^ z e (0,1], z Ф 0, Vje [0,2p], то есть с параметром z = 1 разностное уравнение (3) является спектрально устойчивым.

Тестовый пример 1

utt = uxx + sin(x) sin(t), x e (0, p), t > 0 • u(x, 0) = sin(x), ut (x, 0) = sin(2x), x e [0, p] u(0, t) = u(p, t) = 0, t > 0 Сведем исходную задачу к трем простым:

u(x, t) = u1 (x, t) + u2 (x, t) + u3( x, t), где %(x, t), u2( x, t), u3( x, t) - решения первой (a), второй (b) и третьей (с) задач соответственно: ultt = u1xx + sin(x)sin(t), x e (0, p), t > 0

a) • u1 (x, 0) = 0, u1t (x, 0) = 0, x e [0, p]

u1(0, t) = u1(p, t) = 0, t > 0;

u2tt = u2xx, x e (0, pXt > 0

b) • u2(x,0) = sin(x), u2t (x,0) = 0, x e [0, p]

u2(0,t) = u2(p, t) = 0, t > 0;

u3tt = u3xx, x e (0, p), t > 0

c) • u3 (x, 0) = 0, u3t (x, 0) = sin(2x), x e [0, p]

u3(0,t) = u3(p,t) = 0, t > 0.

Выбираем решение задачи (a) в виде %(x, t) = _/1 (t) sin( x), удовлетворяющее граничным условиям %(0, t) = %(p, t) = 0, подставляя в первое уравнение системы, получим

Л + Í1 = sin t

/1(0) = Л'(0) = 0

общее решение однородного уравнения (с индексом оо) есть f[oo(t) = A sin t + B cos t. Частное решение неоднородного уравнения (с индексом ч) ищем в виде

/1ч(t) = (a + bt)cos t, /1ч = bcos t- (a + bt)sin t, / 1ч = -2bsin t- (a + bt)cos t, подставляем значения /1ч (t), / 1ч в первое уравнение системы:

2

2

2

tcost '

-2b sin t = sin t ^ b = -1/ 2, a = 0, /1(t) = /1oo (t) + /1ч (t) = A sin t + B cos (t)--—, / (0) = /1 (0) = 0

(sin t -1 cos t) (sin t -1 cos t) sin(x) B = 0, A -1/2 = 0, f1(t) = ±---, u1(x, t) = f1(t)sin(x) = ±-^——

Решение задачи (b) ищем в виде щ,(x,t) = /2(t)sin(x), которое подставим в уравнения системы (b)

/2'' + /2 = 0

[ /2(0) = 1, /2(0) = 0 f2(t) = A1 sin t + B1 cost,/2(0) = 1 ^ B1 = 1,/ 2(0) = 0 ^ A = 0 Получим u2 (x, t) = /2 (t)sin x = cos t sin x .

Решение задачи (c) ищем в виде U3 (x, t) = /3 (t)sin (2x), которое подставим в уравнения системы (c)

/3'' + 4 /3 = 0

/з(0) = 0, /3'(0) = 1

/3 (t) = A2 sin (2t) + B2 cos (2t), /3(0) = 0 ^ B2 = 0, / 3 (0) = 1 ^ 2A2 = 1

. . sin (2t) sin (2 x) U3 (x, t) = /3(t) sin (2x) = 2.

Решение исходного примера есть

(sin (t)-1 cos (t)) sin(x) sin (2t) sin (2x) u( x, t) = u1( x, t) + u2( x, t) + u3( x, t) =----+ cos (t) sin (x) +----.

Инициализация задачи. Инициализацией задачи назовем определение второго временного слоя решения по начальным данным задачи, т.е. определения функции x) в задаче (2)

9u( x,0)N

"эГ

u(mh,t) = ф1(отй), mе [0, M], j1(x) = FI j(x) ° u(x,0),y(x) °

Аппроксимируем первую производную решения в начальный момент времени у(х). Рассмотрим два случая инициализации разностной задачи (2):

1) в задаче (2) первое уравнение однородное ф(х) Ф 0, у(х) Ф 0 ;

2) в задаче (2) первое уравнение неоднородное ф(х) = 0, у( х) = 0.

В этих двух случаях инициализация задачи существенно отличается.

Построим квадратурную формулу для функции у(х), связывающую узловые значения функции в одном узле в трех временных слоях, с момента времени t = 0 точную для многочленов максимальной степени

щ (0) =1 (С0и(0) + С1м(х) + С2И(2т) ), и(0 ° 1: ых (0) = 0 =1 (С0 + С1 + С2)

X X

ы(г) ° t: их(0) = 1 =1 (0С0 + хС1 + 2хС2 ), ы(г) ° t2 : их(0) = 0 =1 (02 С0 + х2С1 + (2х)2 С2).

х ^ '

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

fC0 =-3/2

C0 + C1 + C2 = 0 C1 + C22 = 1 ^ C1 + C24 = 0

C1 = 2

C2 =-1/2

1 ( 3

ut (0) I — u0 + 2u1 ^ — u2 I, откуда

1

11 2

2

u(2t) ° u2 = -3u0 + 4u1 - 2tut.

(13)

Формула (13) имеет третий алгебраический порядок погрешности, т.е. погрешность имеет вид О(х3). Получим аналогичную формулу, связывающую четыре временных слоя:

ых(0) =1 (с0ы(0) + С1ы(х) + С2ы(2х) + С3ы(3х)), и^) ° 1: ых (0) = 0 =1 (С0 + С1 + С2 + С3)

X X

и(г) = г: их(0) = 1 =1 (С00+хС1 + 2хС2 + 3хС3), и(г) = г2 : их(0) = 0 = -(02С0 + х2С1 + (2х)2 С2 + (3х)2 С3) х х \ '

и(г) = г3: их(0) = 0 = -(03С0 +х3С1 + (2х)3 С2 + (3х)3 С3) .

Решаем систему линейных уравнений:

С0 + С1 + С2 + С3 = 0 С1 + 2С2 + 3С3 = 1 С1 + 4С2 + 9С3 = 0 С1 + 8С2 + 27С3 = 0

С0 =—11/6 С1 = 3 С2 =—3/2

С = 1/3

их (0) =1 [ —11 и0 + 3и1 — 3 и2 +1 и3 ], х V 6 2 3 у

откуда

3 11 0 1 9 2

и(3х) ° и =— и — 9и +—и + 3хих 2 2 х

(14)

Формула (14) имеет четвертый порядок погрешности, т.е. погрешность имеет вид О(х4). В первой задаче инициализации с использованием формул (3), (13) Ьтп ° 0 для трех временных

слоев имеем:

п+1 + п—1 — 2 п = Лп + п — 2 п \ + г ит + ит 2ит = 21 «т+1 + ит—1 2ит I + Л

+2 (1—2) ит+2 (ит+1+ит—1 )=

= —3ит + 4ит — 2хих ^

^ 2ит—1 — 2 (1 + 2) и/ + 2ит+1 = —2ит — 2хих, т = 1, М — 1.

(15)

Система линейных уравнений (15) представляет трехдиагональную матрицу с коэффициентами

Ат = 2, Ст = 2 (1 + 2), Бт = 2, ¥т = —2ит — 2хих

относительно неизвестных «т— 1, ит, ит+1.

Решить систему уравнений (15) можно методом прогонки, например, с помощью формул [3, с. 44], [7, с. 68], прогонки вперед

Бт

"т_ .. _ Атпт—1 ¥т ____1 л,? 1 1 _ п ,, _ ,1

*т ~ С — Т^Г-, пт = С — , л-, т = 1М — 1 Л0 = 0, п0 = «0

-т АтЛ'т—1

и формул прогонки назад

«т = 1т«т+1 + пт, т = М —1,1

(16)

Условие корректности формул прогонки (условие абсолютного диагонального преобладания) выполнено, так как |Ст| > |Ат| + |Бт| > 0 ^ 2 (1 + 2 )> 22 > 0 [3, с. 35]. Поскольку формула(13) точна для

всех многочленов второй степени, т.е. имеет погрешность О (х3) с третьим алгебраическим порядком, то

и погрешность решения системы уравнений (16) также О (х3), т.е. формулу (15) можно считать начальным этапом инициализации задачи. Используя формулы (3), (14) Ьтп °0, построим второй этап инициализации задачи (2)

п+1 + п—1 — 2 п = Л п + п — 2 п \ + г х2 ^ щ 3 = 1 + 2/'1 — \ 2 + Л 2 + « 2 \ = ит + ит 2ит = 2\Цт+1 + «т—1 2ит I + /т,пх ^ ит = «т + 211 2!ит + 2!«т+1 + ит—11 =

= —ит — 9ит + "2 «т + 3хитх ^ 8ит + 22«/и + 2 («т+1 + «т—1) = _ит + 3хитх =

=8«т+V-2—22 1 (—«т+2 (1—2) «т+2 («т+1+«т—1))+

2

2

0

Z 2 и1 z um-

11 - 5

+z (-um+i+2 (1 - z) um+i+z (um+2+um)- um-i+2 (1 - z) um-i+z (um+um-2)) °

+ (2z(1 - z) - z [5 + 2z))um-1 + ulm |8 - [5 + 2zJ 2 (1 - z) + 2z2 j + (2z(1 - z) - z |5 + 2z^j)ulm+1 + z2ulm+2

- 2z i um+3tumt+zum_i+zum+i о z2um_2 - [ ^ z+4z21 um_!+um (3+z+6 z2) -1—z+4z21 um+i+z2u1

m-1 ^ um | 2, 1

2 J m /и* m—1 tn~ 1 m— 2 1 ^

2 I 1 - . и -2 1..1 2 J m1 ' | 2

Для оптимального параметра z = 1 формула (17) перейдет в формулу

)-(1 z + 4 z 2 j

2

о ап-2-|1 z+4 z 2 j um-1+u\n (3+z+6z2)-1 2z+4z2 | u/+1+z 2u/+2=(3 - 2z)uí+zu/-1+zum+1+t • (17)

m+1 z um+2

0 , O..0

um-2 2 um-1 + 10um 2 um+1 + um+2 = um + um-1 + um+1 + 3tum 1

(18)

Отметим, что трехдиагональная матрица (15) и пятидиагональная матрица (18) систем линейных уравнений в данном случае используются для аппроксимации второго временного слоя, а не для решения основной задачи, как в работах [7, 8].

Линейная система уравнений (19) имеет пятидиагональную матрицу прогонки относительно неизвестных ы^-2, ыи ы]п, ы^+1, ыт+2 , обозначим ее коэффициенты

Ат1 =1, Ат2 = -9, Ст = -10, Вт1 = -Вт2 =1, Рт = ыт-1 + ыт + ыт+1 + 3хытх.

Для решения системы уравнений (18) используем формулы прогонки [7, с. 70], формулы прогонки

вперед:

11m ="

B1m + A2m 12m-1 + A1m11m-212m-1

Cm A1m 11m-2 ^m-1 A1m 12m-2 A2m11m-1

12m =

B

2m

Cm A1m 11m-2 11m-1 A1m12m-2 A2m11m-1

(19a)

A1m 11m-2Vm-1 + A1mVm-2 + A2mVm-1 Fm m = 2 M - 2

Cm - Am 11m-211m-1 - A1m12m-2 - A2m 11m-1

и формулы прогонки назад

im = 11mum+1 + 12mum+2 + vm , m = M - 2,2 (196)

u0 = 110ul +120u2 +v0, ul = 111u2 +121u3 +v1, uM-2 = 11M-2uM-1 +^2M-2uM +vM-2.

Из последней формулы видно, что u0, u! принимают фиксированные значения (краевое условие Дирихле - краевые значения u^, uj, uM -1, uM заданы) при любых соседних узловых значениях, если положить V0 = u0,110 = I20 = 0, V1 = uj, 1ц = I21 = 0 •

Коэффициенты последней системы уравнений не имеют даже нестрогого (абсолютного) диаго-

..........[ 9 1

нального преобладания, так как \Cm\ = 10<\Am1\ + |+ |Bm1| + |Bm^ = 211 + — J = 11. Однако недиагональные элементы матрицы Am1,Am2,Bm1,Bm2 являются знакопеременными и удовлетворяют требованию условного диагонального преобладания:

\Cm\ = 10 - Am1 + Am2 + Bm1 + Bml\ = 2

1 - 9 2

= 7.

Определение 2. Говорят, что квадратная матрица коэффициентов j имеет условное диагональное преобладание, если для любой ее строки

\аи\-

I

j=1, j*i

"í. j

, i = 1, n •

Оказывается, что для корректности формул прогонки вперед (19а) при решении системы уравнений (18) достаточно выполнить условие более слабое, чем абсолютное диагональное преобладание эле-

ментов матрицы |аг-,г-| > V ^ |Ст| > 2(|А2т| + |Л1т|) ^ |—10| V2(|—9/2 +1) = 11, но более сильное, чем условное диагонального преобладания элементов матрицы (для коэффициентов левой части уравне-

п

ния (18)) \а: А >

I

а

j =1, j *i

i, j

Cm\> 2 (| Am| "I Am|)H"1°I > 2 (l"9/2-1) = 7,

1-10 = \Cm\> 2| A2J + 1 Ami = 21-9/2 +1 = 10,

а именно

Утверждение 3. Пусть квадратная пятидиагональная матрица коэффициентов линейной системы уравнений

AmuJ-2 + A2muJ-1 " CmuJ + B1muJ+1 + B2muJ+2 = Fm, m = 2, ^ " 2

является:

1) симметрической: Aim = B2m, A2m = Blm, со знакочередующимися недиагональными коэффициентами для определенности

Am > 0, A2m < 0, Cm < 0;

2) удовлетворяет условию (условию подчинения коэффициентов):

|Cm| > 2|A2m| + |A1m| ^ Cm <-2|A2m|-|A1m| , |A2^ > 2IAm\ ■

Тогда в краевой задаче Дирихле Vo = u0,1ю = I20 = 0,V1 = uj, 1ц = I21 = 0 :

1) 0<11m < 1, -1/2<12m <0, m = 0,n-2;

2) формулы прогонки вперед (19а) для коэффициентов , ^2m, Vm корректны. Доказательство утверждения 3 проведем по индукции:

1) для базы индукции m = 2 имеем:

^12 =■

^22 =

В\2 + A22121 + A12110121 C2 - A12110111 - A12120 - A22111 C2

= A22 > 0,

A

22

C2

A

22

A

22

B

22

= A^ < 0,

A12

2| A22 + | A1 2 2IA22 IA121 ' '

< 1,

A2 1 1 < - < -

C2 - A12110111 - A12120 - A22111 C2 2\A221 + | A2I 5|A2I 5 2

Sign(l12) = Sign

A

a22 V C2 у

= 1, Sign(l22) = Sign

A

a12 V C2 у

= -1 .

База индукции проверена. Пусть верно

\hk | < 1, I < 2, Sign(l1k) = Sign

A

A2k V Ck у

Ak

Тогда =1

\B1m + A2m12m-1 + Am11m-212m-11

= 1, Sign(l2k) = Sign I I = -1, "k = 2, m -1,

I A2m + A2m12m-1 + A1m11m-212m-11

|Cm Am 11m-211m -1 A1m12m-2 A2m11m-1| |Cm A1m11m -211m-1 A1m12m -2 A2m11m-1\

-1/2 < l2m-1 < 0,1/2 < l2m-1 +1 < 1, 0 < 1^-1 < 1, - | A2m| < A2m (1 + W1) < -1 A2m | /2,

-1 Am\ / 2 < Am11m-212m-1 < 0 , -1 A2m | -1 A1m | / 2 < A2m (1 + 12m-1 ) + Am11m-212m-1 < -1 A2m | / 2 < 0 , -| Am| < -A1m11m-211m-1 < 0 , 0 < -A1m12m-2 < \Am\ /2, 0 < -A2m11m-1 < | A2^ , Cm - A1m < Cm - Am11m-211m-1 - A1m12m-2 - A2m11m-1 < Cm + Am / 2 + A2m < 0 ,

Am /2

|Cm| +| Am\/2

<|11m| =

A2m(1+) + A1m11m-212m-1

Cm A1m11m-211m-1 A1m12m-2 A2m11m-1

\A2J +| Am|/2 < \A2J +| Am I /2 = ^

|Cm| |Am|/2 |A2m| |A2j +| Am| I2

h I \B2m\ . Am Hlm 2 1

i2m = 1-Г — |-i—i-i-£-i-Г = _ < _ •

\Cm - A1m11m-211m-1 - Am12n-2 - A2m11m-l\ \A2m\ + |Am| /2 5/2|Am| 5 2

Из предыдущих оценок следует, что знаменатель и числители формул заключены в пределах Cm - Am11m-211m-1 - Am^2m-2 - A2m11m-1 — Cm + | Am | / 2 +1 A2m | — -1 A2m |-1 Am | / 2 < 0 ,

A2m (1 + ^2m-1) + A1m11m-212m-1 — -1 A2m 1/2 < 0, B2m = Am > 0 •

Тогда с учетом предыдущих неравенств получим, что 0 — ^1m — 1, -1/2 — 12m — 0 ; 2) ввиду доказанной первой части утверждения 3 имеем

0<11т < 1, -1/2<12т <0, т = 0,п-2

Ст -|Ат| /2 < Ст - Ат 11т-211т-1 - Ат^2п-2 - А2т^1т-1 < Ст + |Ат| /2 + | А2т| < - | А2т\-|Ат| /2 < 0

Поэтому знаменатель во всех трех формулах прогонки (20) , %2т, пт строго меньше нуля, т.е. в ноль не обращается, а формулы прогонки корректны. Утверждение 3 доказано.

Замечание 3. В работе [4] показано, что достаточными условиями корректности формул прогонки симметричной пятидиагональной матрицы является полуторное абсолютное диагональное преобладание ее элементов. В данном случае утверждение 3 указывает достаточные условия корректности прогонки той же матрицы, но более слабые, чем даже абсолютное диагональное преобладание ее элементов.

Рассмотрим равномерную норму погрешности процесса инициализации, т.е. разности численного

решения ит, полученного по формулам (13), (15), (16), и точного решения первой части тестового при-

1 , , sin (2t) sin (2hm) --1 1

мера 1 (для однородного уравнения) um = costsin(hm)+--, m = 0,M, A — -■

Норма Чебышева определяется формулой Д= max |Dm| = max

m=1,M m=1,M

-1 . 1

, программа возвращает

норму погрешности norma C = 5.174140550234796E-006 при M = 200 и norma C = 6.465012549750071E-007

I# лпп Д(200) 5.174140550234796E-006 0 nn„ о3

при M = 400, тогда—;-- =-= 8.003 » 2 алгебраическии порядок погрешно-

Д(400) 6.465012549750071E-007

сти алгоритма(13), (15), (16) равен p = 3, т.е. подтверждает приближение формулы (13) O(t3).

Последовательный алгоритм инициализации сначала на первом этапе по формулам (13), (15), (16), а затем по формулам (14), (18), (19) с помощью программы дает значение равномерной нормы norma C = 8.299974534053955E-008 приM = 200 и norma C = 5.173976233563415E-009 M = 400, тогда Д(200) 8.299974534053955E-008 лгпл „4 ^

—;-- =-= 16.04 » 2 . То есть алгебраический порядок погрешности алгоритма

Д(400) 5.173976233563415E-009

(13), (15), (16) и (14), (18), (19) равен четырем ( p = 4 ), что подтверждает формулу (14), а именно O(t4).

Построим явную разностную схему для решения однородного уравнения, учитывая (3) для параметра z = 1

n+1 о n . Í n . n П—1 n . n П—1 О ЛГ 1 Л/f 1 /ОПЧ

= 2Um + Z Mm+1 + um—1 — 2Um — um = Um+1 + um—1 — um , n = 2, N, m = 1, M — 1. (20)

Отметим, что формула (20) имеет бесконечный порядок аппроксимации, т.е. не вносит дополнительной погрешности аппроксимации дифференциального уравнения разностным уравнением.

Программа с учетом алгоритма инициализации по формулам (13), (15), (16), а затем по формулам (14), (18), (19) явной формуле (20) для первой части задачи (однородного уравнения) тестового примера 1

Щи = и1хс, * е (0, я), г > 0 и1 (х, 0) = sin(х), иц (х,0) = sin(2х), * е [0, я] и1(0, г) = и1(я, г) = 0, г > 0

sin (2t) sin (2x)

с точным решением x,0) = cos tsin x +----, xe [0, p], t = 3p» 9.42477796 возвращает равномерную норму относительной погрешности 7.214666558239686E-016 при M = 50 и норму 5.275261571720079E-015 при M = 100, t = 9.42477796076938.

Это означает что, во-первых, норма погрешности по формуле (20) не зависит от шага сетки при z = 1, во-вторых, с увеличением вычислений при меньшем шаге накапливается только большая абсолютная погрешность. В программе предусмотрен масштабный параметр l = 10 (в приведенном примере).

2

Данный параметр позволяет экономить время счета в l = 100 раз, а, главное, число арифметических

2

операций в l2 = 100 , т.е. уменьшить абсолютную погрешность, пропорциональную числу арифметических операций.

Сначала инициализация проводится с шагом t по времени по формулам (13), (15), (16), а затем по формулам (14), (18), (19), далее по явной формуле (20) с минимальными шагами сетки (h, t) решается

волновое уравнение ufi^1 = м^+х + U\m-\ - un-1, n = 2, l, m = 1, M -1 во временном промежутке [0, l -t],

п е 0,1. Среди решения в слое щ (шк, I • х), т = 0, М и среди начального слоя щ (шк, 0), т = 0, М выбираются узлы более редкой сетки хш1 = к • I • ш1, ш1 = 0, М /I, п = х • I • п1, п1 = 0,1, и решение на той же

сетке щ (ш1 • к • I,0), щ (ш1 • к • I, I • х), ш1 = 0, М /1.

Далее используется формула (20) с крупным вектором шага (I• к,I•х),1 = 7 = х2а2 /к2 ^ х = к/а

= Mm+1 + um- - мт1, n1 = 2, N / -l, m1 = 1, M /1 -1. (21)

Рассмотрим вторую задачу инициализации c использованием условий на трех временных слоях j(x) = 0, x e [a, b], ut (x,0) = y(x,0) = 0, f (x, t) Ф 0,

U2tt = u2xx + sin(x)sin(t), x e (0, p), t > 0

(sin t -1 cos t) sin( x)

u2(x,0) = 0, u2t(x,0) = 0, xe [0,p] u2(x,t) = ---.

u2(0, t) = u2(p, t) = 0, t > 0, Формула (13) перейдет в u(2t) ° u2 = -3u0 + 4u1 - 2tut = 4u1, с учетом (3) получим

n+1 _ n-1 . o fi \n in n \ , x „.2 2 _ 0 ~ (i 1 i ( 1 1 ^ i ^ т2 a 1 - v

um =-um + 2 (1 - z) um + z\um+1 + um-1 ) + fm,nx , um =-um + 2 (1 - z) um + z\um+1 + um-1 J + fm,1x = 4um ^

^ zum-1 - 2 (1 + z) um + zum+1 =- fm,1t + um =-fm,1t2. (22)

В трехдиагональной матрице системы уравнений (22) коэффициенты

Am = z, Cm = 2 (1 + z), Bm = z, Fm = um - fm,1t

обеспечивают условие корректной прогонки вперед |Cm| > |Am| + |Bm| > 0 ^ 2(1 + z) > 2z > 0 [3, с. 45].

Программа с использованием метода прогонки для алгоритма инициализации (23) относи-

(sin t-t- cos (t)) sin( x)

тельно точного решения x, t) =--- дает равномерную норму погрешности

norma C = 1.291625408111898E-006 при M = 200 и norma C = 1.614815612291833E-007 при M = 400,

Д(200) 1.291625408111898E-006 _QQ_ . ,3 -

откуда —;-- =-= 7.9986 » 2 алгебраический порядок погрешности p = 3.

Д(400) 1.614815612291833E-007

Рассмотрим вторую задачу инициализации на четырех временных слоях, с учетом формул (14), (3) получим

3 11 0 1 9 2 1 9 2

u(3t) ° u = — u - 9u1 + - u + 3tut = -9u1 + - u , 2 2 t 2

откуда

i n+1 = -, n-1 + 9/1 - \.n + in + n \ + f t2 ^^ ,3 = - 1 + 2/1 - \. 2 + Í 2 + u 2 \ + r t2 = um = um + 2li zJum + zlum+1 + um-1 I + Jm,nt ^ um = um + ^um + ^um+1 + um-1 I + fm,2t =

192 1 I 5 | 2 /2 2 \ 2

= -9um + ^Um ^ 8um -1 ^ + 2z I Um + z (um+1 + um-1) = ~fm,2t ^

8um- | 2+2z j (-um+2 (1 - z) um+z (um+1+um-1)+fm,1t2)+ +z (-um+1+2 (1 - z) um+1+z (um+2+um)+fm+1,1t2 - um -+2 (1 - z) um-1+z (um+um-2)+fm-ut2)=- fm,2t

21 1 1 2) 1 1 / 2 \ 1 1 2) 1 21

^ z um-2 -1 ^z + 4z Ium-1 + um (3 + z + 6z )-1 ^z + 4z Ium+1 + zum+2 =

= - | fm,2 - ^| + 2z j fm,1 + z ( fm-1,1 + fm+1,1)|t2 . (23

Для параметра z = 1 уравнение (23) переходит в уравнение (24)

9

1 91 1П 191 1 = \ f 9 f f f j 2

um-2 - 2 um-1 + 10um - 2 um+1 + um+2 = -1 fm,2 - ^ fm,1 + Jm-1,1 + Jm+1,1 It .

(24)

2 m—1 m 2 m-r 1 m-r2 I ^ m,2 2 -

Отметим, что коэффициенты пятидиагональной матрицы системы уравнений (24) Am1, Am2,

C B B

Cm, Bm1,

т, Вт1, Вт2 •

= =9 = =9 = = Г 9 ^ 2

Ат1 = 1, Ат2 = - Ст = -10, Вт1 = -2, Вт2 = 1, *т = -1 /т,2-^ ?т,1 + fm-1,1 + /т+1,1 11

совпадают с коэффициентами матрицы системы (18) за исключением ¥т. Однако корректность формул прогонки для пятидиагональной матрицы (18), которая определяется только коэффициентами Ат1, Ат2, Ст, Вт1, Вт2 , доказана с помощью утверждения 3, следовательно, корректность формул прогонки для системы уравнений (24) также доказана.

Программа с использованием формул (16), (22), а затем (19), (24) в двухэтапном алгоритме инициали-

зации относительно точного решения u2 (х, t) =

(sin t-t- cos (t)) sin( x)

2

возвращает равномерную норму

погрешности norma C = 2.028797013322345E-008 при M = 200, а norma C = 1.268260177248376E-009 при

/rnn Д( 200) 2.028797013322345E-008 , 4 -

M = 400, откуда —;-- =-= 15.997 » 2 = 16 , т.е. алгебраический порядок

Д(400) 1.268260177248376E-009

погрешности для указанного алгоритма инициализации методом прогонки равен четырем ( p = 4 ).

Наконец, нужно написать основную рекуррентную разностную формулу для решения неоднородного уравнения с нулевыми начальными условиями. Для параметра z =1 преобразуем формулу (3)

1 n+1 = -, n-1 + 2/1 - \,п + (,п + ,п \ + f t2 = -. n-1 + п + п + j t2 + d(, n ) (25)

um = um +2Л1 z)um + zlum+1 + um-1 I + Jm,nt = um + um+1 + um-1 + Jm,nt + R(um). (25)

В формуле (25) невязку аппроксимации неоднородного волнового уравнения (1) разностным уравнением (2) запишем согласно формуле (8) для z = 1 n = 2, N, m = 1, M -1

u n+1 = un + u n - u n-1 + j t2 + 2t um = um+1 + um-1 um + Jm,nt + 41

2 Э fm,n Э fm

a ----

Эх2

Эt2

6!

4 " Jm,n 2 Э fm,n Э fm

a ----+ a --—— + -

Э4 fm

Эх4

22

Эх2Эг

Эt4

2t_

" 8!

\ Э6 f э6 f э6 f э6 f

6 Э Jm,n 4 Э Jm,n 2 Э Jm,n Э Jm,n

a --—+ a —;—— + a —-—т+-

Эх

Эx4Эt2

Эx2Эt4

2t

10 \ э8 f Э8 f Э8 f

" Vm,n 6 Jm,n . 4 Jm

10!

8 Э Jm,n 6 Э Jm,n 4 Э Jm,n 2 Э fm,n Э fm,n

Эх8

+a

+a

+a

Эх6Э^ Эx4Эt4 Эx2Эt6 Эt8

2t

121 эЮ f Э10 f Э10 f Э10 f Э10 f Э10 f

10 Э fm,n 8 Э fm,n 6 Э fm, n 4 Э fm, n 2 Э fm, n Э fm

12!

Эх

10

+a

Эx8Эt2

+a

Эx6Эt4

+a

Эx4Эt6

+a

Эх2Э^

Эг

10

+...

2

+

+

+

Для примера 1 /(х, г) = х)8ш(г), а = 1 согласно формуле (26) получим явную разностную схему

^2 о„4 „„6 /|_8 с„10 .„12 Л

n+1 = п + п - u n-1 + 2 f 4m = um+1 + um-1 um + 2 fm

- 3t_ - 4t_ 5t_ -

2! 4! 6! 8! 10! 12!

/

, ,, t2 t4 t6 t d,, t2 2t4 3t6 4t8 5t10 6t12

Поскольку 1 - cos(t) =----1---... ^--(1 - cos(t)) =----1-----1-----+..., то

2! 4! 6! 2 dtK ; 2! 4! 6! 8! 10! 12!

um+1=um+1+um_1 - um-1+1 sin (t) fm n=um+1+um_1 - um1+1 sin (t) sin(mh) sin(nt), (27)

n = 1, N, m = 1, M-1.

В формуле (27) для бесконечного функционального ряда получена производящая функция т Sin (т) fm,n .

Аналогично алгоритму (21) можно провести укрупнение шага в неоднородном уравнении. Сначала инициализация проводится с шагом т по времени по формулам (16), (22), а затем по формулам (19), (24), далее по явной формуле (26) с минимальными шагами сетки (h, т) решается волновое уравнение

иПт1 = u"m+i + M"m-i - u2m + т sin (т) sin(mh) sin(nt), n = 2,11, m = 1, M -1 во временном промежутке [0, l • т],

n e 0, l. Среди решения в слое u2 (mh, l -1), m = 0, M и среди начального слоя u2 (mh, 0) ° 0, m = 0, M вы-

бираются узлы более редкой сетки xm1 = h -1 - m1, m1 = 0, M /1, tn1 = t -1 - n1, n1 = 0,1 и решение на ней

Далее используется формула (27) с крупным шагом (l-h,l-1),1 = z = t2a2 /h2 ^ t = hIt

u2(m1- h -1,0) ° 0, u2(m1- h -1, l -t), m1 = 0, M /1

2 2 2

' a

urn = u2nm1+1 + u2nm1_1 - unmt +1 -1 - sin (t -1) sin(m1 - h -1) sin(n1 -1 -1), n1 = 2, N /1, m1 = 1, M /1 -1. (28) Программа с учетом алгоритма укрупнения шага для тестового примера 1 u2tt = u2xx + sin(x) sin(t), x e (0, p), t > 0,

(sin t -1 cos t) sin( x)

u2(x,0) = 0, u2t (x,0) = 0, x e [0, p], u2(x, t) ^^--

u2(0,t) = u2(p,t) = 0, t >0,

(sin t -1 cos t) sin( x)

с точным решением u2(x,t) =---, xe [0,p], t = 3p» 9.42477796 возвращает равномерную норму относительной погрешности norma C = 3.769546289141946E-015 при M = 100,l = 10, t = 9.42477796076938.

(sin t -1 cos t) sin( x) sin (2t) sin (2x) Для суммы решений u( x, t) =----+ cos t sin x +---- программа возвращает равномерную норму относительной погрешности norma C = 7.578604648156152E-016. Следовательно, мы можем сказать, что приведенные нами алгоритмы численного решения волнового неоднородного уравнения на отрезке имеют двойную точность относительной погрешности 10-16.

Таблица 1.

x exact numerical

0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000

0.314159265358979 1.14719128466915 1.14719128466915

0.628318530717959 2.18208749344321 2.18208749344321

0.942477796076938 3.00338577486150 3.00338577486150

1.25663706143592 3.53069173081718 3.53069173081718

1.57079632679490 3.71238898038469 3.71238898038469

1.88495559215388 3.53069173081718 3.53069173081718

2.19911485751286 3.00338577486150 3.00338577486150

2.51327412287183 2.18208749344321 2.18208749344321

2.82743338823081 1.14719128466915 1.14719128466915

3.14159265358979 4.546215133239269E-016 0.000000000000000E+000

В таблице используются следующие обозначения: x - координата узла; exact, numerical - точное и численное решения для примера 1.

Тестовый пример 2.

utt = uхх + sin(x)t, х е (0, я), t > 0

• u(х,0) = sin(х), ut (х,0) = sin(2x), х е [0, я]

u(0, t) = u(p, t) = 0, t > 0. Можно проверить, что решением последней задачи является функция

sin (2t) sin (2х)

u( х, t) = cos t sin x + -

2

- + sin x (t - sin t).

Применяя формулу (27) для неоднородной части f (х, t) = sin(x)t, a = 1, получим

n+1

= un + un — , n-1 + 2 f = um+1 + um—1 um + 2 Jr.

\

( „. 2 „.4 „.6 _8 _10 12 t t t t t t

2! 4! 6! 8! 10! 12!

v

n-1

= um+1 + um—1 — um 1 + 2(1 — cos(t)) fmn = um+1 + um—1 — um 1 + 2(1 — cos(t))sin(mhxnt), m = 1,M — 1, n = 1,N — 1.

В данном случае мы получили другую производящую функцию 2(1 — cos(t)) fm n. Повторяя вычислительные алгоритмы для примера 1, получим, что программа при M = 100, l = 10, t = 9.42477796076938

sin (2t) sin (2 x)

для примера 2 возвращает для суммы решений u(х, t) = cos tsin x + -

2

- + sin x (t — sin t) равно-

мерную норму относительной погрешности norma C = 3.005569483512412E-015.

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

u n+1 = un + un — u n—1 + 2 f um = um+1 + um—1 — um + 2 fm

10

t

12

\

( 2 4 6 8 ttttt

2! 4! 6! 8! 10! 12!

v /

дает значение относительной погреш-

ности norma C = 2.337665153842987E-015, т.е. не хуже, чем с помощью производящей функции. Очевидно, применение производящих функции экономнее для решения волнового уравнения в двухмерном и трехмерном случаях.

Таблица 2. (для примера 2)

m

x exact numerical

0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000

0.314159265358979 2.60339956371325 2.60339956371325

0.628318530717959 4.95196023917890 4.95196023917890

0.942477796076938 6.81578854409794 6.81578854409794

1.25663706143592 8.01243997792951 8.01243997792951

1.57079632679490 8.42477796076938 8.42477796076938

1.88495559215388 8.01243997792951 8.01243997792951

2.19911485751286 6.81578854409794 6.81578854409794

2.51327412287183 4.95196023917890 4.95196023917889

2.82743338823081 2.60339956371325 2.60339956371325

3.14159265358979 1.031703 662030091E-015 0.000000000000000E+000

Тестовый пример 3.

utt = uxx + sin(x)ch(t), xe (0,p), t > 0 • u(x,0) = sin(x), ut (x,0) = sin(2x), x e [0, p] u(0, t) = u(p, t) = 0, t > 0.

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

srn (2t) srn (2 x) + sin (x)(ch (t)-cos (t)) u(x, t) = cos (t) sin (x) + —-———-—- +--—---,

отметим, что неоднородная часть уравнения f (x, t) = sin(x)ch (t), a = 1 удовлетворяет однородному вол-

2 d2 f d2 f 2

новому уравнению (9) комплексного аргумента a —— +--— = 0, a > 0 . В этом случае согласно форму-

dx dt

ле (11) запишем невязку

¥ _4k+2 d4k f

u n+1 = n + n - n-1 + f t2 + _t_ Jm,n = n + n - un-1 + f t2 +

um = um+1 + um-1 um + fm,nt + 4k = um+1 + um-1 um + fm,nt +

k =1 ( 4k + 2)! dt4k

2t

6 d4 f ot10 d8 f ot14 d12 f

m, n 2 t m,n 2 t m

6! dt4 10! dt8 14! dt

12

+ = u n + n - u n-1 + 2 f - +... = um+1 + um-1 um + 2/m

^ 2 _6 _10 14 ^ t t t t

-+-+-+-+...

2! 6! 10! 14!

= um+i + um - um 1 + fmn (ch(t) - cos(t)), так как ch(t) - cos(t) = — (et + e t)- cos(t) =

~(ex + e T)-cos(t):

2 ( t2 t4 t2k - 1 + — + — +... + -—- +. 2 [ 2! 4! (2k)!

(„ 6 „10 14

Л ( t2 t4 t.

i t t , .л t

1--+-+ ... + (-1)k 7-г

2! 4! (2k)!

2k \ л

= t2 + 2

t t t t

—+—+—+...+---

6! 10! 14! (2 + 4k)!

2+4k ^ +...

sin (2t) sin (2x) sin (x) (ch (t) - cos (t)) Для точного решения u( x, t) = cos t sin x +-----+ 5--- программа

в примере 3 для M = 100, l = 10, t = 9.42477796076938 возвращает норму относительной погрешности norma C = 2.557857845110683E-015. Производящая функция равна fm n (ch(t) - cos(t)).

Таблица 3. (для примера 3)

x exact numerical

0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000

0.314159265358979 957.152937976024 957.152937976022

0.628318530717959 1820.61307750630 1820.61307750629

0.942477796076938 2505.85892405305 2505.85892405305

1.25663706143592 2945.81383976773 2945.81383976773

1.57079632679490 3097.41197215405 3097.41197215405

1.88495559215388 2945.81383976773 2945.81383976773

2.19911485751286 2505.85892405305 2505.85892405305

2.51327412287183 1820.61307750630 1820.61307750630

2.82743338823081 957.152937976024 957.152937976023

3.14159265358979 3.793110381505353E-013 0.000000000000000E+000

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

program wave

integer(8), parameter::n=100,n1=10,ll=n/n1,m=3;integer(8)::i,j,k real(8):: num(0:n+1,0:m*n+1),num0(0:n+1,0:m*n+1)

геа1(8)::раг(0:п+1)^ит^Дау2,ГО0(0:п+1)

геа1(8)::^1(0:п+1),1(0:п+1),ГО(0:п+1),аа33(0:п+1), ^2(0т+1)

геа|(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)

геа1(8)::а1(0:п+1), а2(0:п+1), Ь1(0:п+1), Ь2(0:п+1),аа11(0:п+1),сс11(0:п+1),ЬЬ11(0:п+1)

real(8)::eps00(0:n+1),otv00(0:n+1), I1(0:n+1),l2(0:n+1),res3(0:n+1)

геа1(8):: аа1(0:п+1),ЬЬ1(0:п+1),аа2(0:п+1),ЬЬ2(0:п+1),аа3(0:п+1),ЬЬ3(0:п+1)

real(8)::max1,max2,max3,max4,max44,epss(0:n+1);real(8)::max5,ch,t,yy,max55,mm,tay1,c1

real(8)::u1,u0,f1,f2,fan,z,vel,x,y,a,Ь,c,d,pi,h1,tay,tt,x1,x2,x3,x4,hh,fy

!ch(t)=(dexp(t)+dexp(-t))/2d0;

fУ(x,tay)=dsin(x)*dsin(tay);ul(x)=dsin(2d0*x);u0(x)=dsin(x)

fan(x,t)=dsin(x)*dcos(t)+dsin(2d0*x)*dsin(2d0*t)/2d0;f1(x,tay)=-2d0*u0(x)-2d0*tay*u1(x)

f2(x,t)=dsin(x)*(dsin(t)-t*dcos(t))/2d0;pi =2d0*dasin(1d0);a=0d0;Ь=pi;

Z=1d0;vel=1d0;max1=-100d0;max2=-100d0;max4=-100d0;max44=-100d0

max5=-100d0;max55=-1000d0;mm=-100d0;h1=(Ь-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;ЬЬ(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;Ь2(k)=1d0;Ь1(k)=-4.5d0;ccc(k)=-10d0;aa1(k)=1d0

ЬЬ3(k)=1d0;f11(k)=-fy(x,tay)*tay*tay;aa11(k)=z;ЬЬ11(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))

1(к)= ЬЬ(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(i);par(i)=fan(x,tay);eps(i)= par(i)-res1(i)

if( eps(j)<=0d0 )then;eps(j)=-eps(j);else;endif;enddo;do j=0,n

if( eps(i)>=max1 )^п;тах1= eps(j);endif;enddo;print*,"norma С1=",тах1

do k=0,n;res2(k)= res1(k);enddo

I1(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 I1(j)=(a1(j)*I1(j-2)*I2(j-1)+a2(j)*I2(j-1)+Ь1(j))/(ccc(j)-a1(j)*I1(j-2)*I1(j-1)-a2(j)*I1(j-1)-a1(j)*I2(j-2)) I2(j)=Ь2(j)/(ccc(j)-a1(j)*I1(j-2)*I1(j-1)-a1(j)*I2(j-2)-a2(j)*I1(j-1))

nu(j)=(a1(j)*I1(j-2)*nu(j-1)+a2(j)*nu(j-1)+a1(j)*nu(j-2)-ff(j))/(ccc(j)-a1(j)*I1(j-2)*I1(j-1)-a1(j)*I2(j-2)-а2(^)*И(^-1)); enddo

do j=n-2,0,-1;res2(j)=I1(j)*res2(j+1)+I2(j)*res2(j+2)+nu(j);enddo;do j=0,n x=a+h1*dfIoat(j);par(j)=fan(x,tay);eps(j)= par(j)-res2(j);if( eps(j)<=0d0 )then eps(j)=-eps(j);eIse;endif;enddo;do j=0,n;if( eps(j)>=max1 )^п;тах1= eps(j);endif;enddo do j=2,n-2;if(eps(j)>=max2)then;max2=eps(j);endif;if(mod(j,n1)==0)then;endif;enddo print*,"norma С2=",тах2^о i=0,n;if(mod(i,nl)==0)then;endif;enddo

do j=0,n1;x=a+h1*dfloat(j);num0(0,j)=0d0;num0(n,j)=0d0;num0(j,0)=u0(x);num0(j,1)=res2(j) enddo;do j=1,n1-1;do 1= 1,П-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*dfIoat(n1);t=dfIoat(n*m)*tay;num0(0,l)=0d0;num0(II,1)=0d0;num0(0,0)=0d0;num0(II,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,II*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,II;eps0(i)=num0(i,II*m)-otv(i);if(eps0(i)<0d0)then;eps0(i)=-eps0(i);endif if(eps0(i) > max44)then; max44=eps0(i);endif

s=s+aЬs(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*dfIoat(k);nu(k)=(aa11(k)*nu(k-1)-f11(k))/(cc11(k)-aa11(k)*I(k-1))

l(k)= ЬЬ11(k)/(cc11(k)-aa11(k)*I(k-1));enddo;do k=n-1,1,-1;res1(k)= I(k)*res1(k+1)+nu(k);enddo

do j=0,n;x=a+h1*dfloat(j);par(j)=f2(x,tay);eps(j)= par(j)-res1(j)

if( eps(j)<=0d0 )then;eps(j)=-eps(j);eIse;endif;enddo;do j=0,n

if( eps(j)>=max4 )then;max4= eps(j);endif;enddo;print*,"norma C101=",max4;do k=0,n

^2(Ю= res1(k);enddo;I1(0)=0d0;I1(1)=0d0;I2(0)=0d0;I2(1)=0d0

I1(n)=0d0;I2(n)=0d0;I1(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

I1(j)=(a1(j)*I1(j-2)*I2(j-1)+a2(j)*I2(j-1)+Ь1(j))/(ccc(j)-a1(j)*I1(j-2)*I1(j-1)-a2(j)*I1(j-1)-a1(j)*I2(j-2)) I2(j)=Ь2(j)/(ccc(j)-a1(j)*I1(j-2)*I1(j-1)-a1(j)*I2(j-2)-a2(j)*I1(j-1))

nu(i)=(a1(j)*l1Cj-2)*nu(i-1)+a2(j)*nu(j-1)+a1(j)*nu(j-2)-f22(j))/(ccc(j)-a1(j)*l1(j-2)*l1(j-1)-a1(j)*l2(j-2)-a2(j)*l1(j-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(j);par(j)=f2(x,tay);eps(j)= par(j)-res2(j);if(mod(j,n1)==0)then;endif if( eps(j)<=0d0 )then;eps(j)=-eps(j);else;endif;enddo;do j=0,n

if( eps(j)>=max5 )then;max5= eps(j);endif;enddo;do j=0,n;if(eps(j)>=max5)then;max5=eps(j);endif;enddo print*,"norma C202=",max5;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)+fy(x,t)*tay*dsin(tay);enddo;enddo do i=0,n;res3(i)=num(i,n1);enddo;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)+fy(x,t)*tay1*dsin(tay1);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;otv00(i)=num(i,ll*m)+num0(i,ll*m);x=a+hh*dfloat(i);res3(i)=f2(x,t)+fan(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,res3(i),otv00(i) enddo;print*,"eps=(norma C)/abs(m)_=",mm*dfloat(ll)/s,"abs(m)_=",s/dfloat(ll);end program wave

В работе получены результаты:

1. Показано, что явная разностная схема задачи (2) аппроксимирует одномерное однородное волновое уравнение на отрезке задачи (1) с оптимальным параметром сетки z = 1 с бесконечным алгебраическим порядком аппроксимации, при этом явная разностная схема спектрально устойчива.

2. Для неоднородного волнового уравнения с параметром z = 1 невязка разностной схемы может быть учтена и явно выражена только через двойную сумму от частных производных правой части неоднородного уравнения (7).

3. Введено понятие волнового уравнения комплексного аргумента и рассмотрены свойства его решения. Если правая часть неоднородного волнового уравнения является решением однородного волнового уравнения комплексного аргумента, то невязка упрощается от суммы по двум индексам до суммы по одному индексу (11), что качественно уменьшает число вычислений.

4. Построен алгоритм инициализации задачи, то есть аппроксимация второго временного слоя решения по известным начальным данным (начальным смещению и скорости точек отрезка).

5. Инициализация сводится к решению СЛАУ с симметричной трехдиагональной матрицей

с погрешностью аппроксимации O (t3), затем к симметричной пятидиагональной матрице с погрешностью аппроксимации O(t4). Найдены ослабленные достаточные условия корректности формул прогонки для пятидиагональной матрицы (условия более слабые, чем диагональное преобладание элементов матрицы).

6. Предложен алгоритм масштабирования узловой решетки по числу узлов в l раз, что сокращает число вычислений в сотни раз при сохранении вычислений с двойной точностью (15-16 первых десятичных знаков результата).

7. Программа тестировалась тремя построенными аналитическими примерами. Показано, что в примерах бесконечный функциональный ряд в разностной схеме относительно временного шага t можно заменить производящей функцией от t, что качественно снижает время вычислений и увеличивает их точность. С помощью программы показано, что построенный алгоритм дает решение 3 примеров с двойной точностью. В самом общем случае с помощью формулы (26) в виде двойного ряда можно достичь 14-го порядка погрешности (приближение O(t14)), причем формула (26) имеет ту же относительную погрешность решения, как и с применением производящей функции, т.е. с двойной точностью.

8. Полученные результаты можно применить, например, при анализе физических систем, рассмотренных в работе [9].

ЛИТЕРАТУРА

1. Пикулин, В.П. Практический курс по уравнениям математической физики : учеб. пособие / В .П. Пикулин, С .И. Похожаев. - М. : Наука,1995. - 224 с.

2. Методы ускорения газодинамических расчетов на сконструированных сетках / К.Н. Волков [и др.]. -М. : Физматлит, 2013. - 536 с.

3. Самарский, А.А. Численные методы решения обратных задач математической физики : учеб. пособие / А.А. Самарский, П.Н. Вабищевич. - М. : Изд-во ЛКИ, 2014. - 480 с.

4. Бахвалов, Н.С. Численные методы в задачах и упражнениях : учеб. пособие / Н.С. Бахвалов, А.В. Лапин, Е.В. Чижонков. - 2-е изд., перераб. и доп. - М. : БИНОМ. Лаборатория знаний, 2010. -240 с. : ил.

5. Киттель, Ч. Введение в физику твердого тела : учеб. руководство : пер. с англ. / Ч. Киттель. - 4 изд. -М. : Наука, 1978. - 792 с.

6. Вакуленко, С.П. К методу оценки состояния железнодорожного полотна / С.П. Вакуленко, К.А. Во-лосов, Н.К. Волосова // Мир транспорта. - 2016. - Т. 14, № 3 (64). - С. 20-35.

7. Пастухов, Д.Ф. Аппроксимация уравнения Пуассона на прямоугольнике повышенной точности / Д.Ф. Пастухов, Ю.Ф. Пастухов // Вестник Полоцкого университета. Серия С, Фундаментальные науки. - 2017. - № 12. - С. 62-77.

8. Пастухов, Д.Ф. Задача построения поля линий тока по температурному разрезу / Ю.Ф. Пастухов, Д.Ф. Пастухов // Вестник Полоцкого университета. Серия С, Фундаментальные науки. - 2015. -№ 4. - С. 27-36.

9. Пастухов, Д.Ф. Классификация профилей температуры в плюс-минус односантиметровом слое от поверхности раздела геотермального озера / Д.Ф. Пастухов // Вестник Московского университета. Серия 3, Физика. Астрономия. - 1995. -Т. 36, № 6. - С. 84-89.

Поступила 28.02.2018

OPTIMUM PARAMETER TO APROXIMATIONS RAZNOSTNOY SCHEMES OF THE WAVE EQUATION ON LENGTH

D. PASTUKHOV, Y. PASTUKHOV, N. VOLOSOVA

The Offered algorithm of the decision initial- marginal problem for lumpy wave equation on length with double accuracy. The Algorithm is founded on choice of the optimum parameter, providing endless algebraic order to approximations to uniform equation. For decision of the system of the linear equations with five diagonal matrixes with marginal condition Dirihle is proved sufficient conditions to correctness molded racing onward more weak, than condition of the diagonal prevalence her(its) element. Using the algorithm of the integration cell nets and use the method producing function gives double accuracy to relative inaccuracy of the decision even rough net with number of the nodes equal 300.

Keywords: method producing function, initializing the problem, weak sufficient conditions to correctness molded racing symmetrical five diagonal matrixes.