∂∂ut + α(x, t)∂∂ux = f (x, t).

Постановка корректных задач для последнего уравнения принципиально отличается от постановки задачи (8.4.20), (8.3.2), (8.3.3). Например, при α(x,t)> 0 , кроме начального ус-

ловия (8.3.2), необходимо задать краевые условия только на левой стороне прямоугольника D , а правая сторона должна быть свободна от каких-либо условий на искомое решение u(x,t). В связи с этим могут возникать существенные трудности при численном решении

уравнения переноса с большими значениями α a02 . В частности, аппроксимация первой про-

изводной решения по x центральным разностным отношением	∂u			=	ui+1, j −ui−1, j	+ O(h2 )
	∂x		i, j		2h

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

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

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

∂2u

−a2 ∂2u

= f (x,t),

(8.5.1)

∂t 2

0 ∂x2

которое описывает поперечные колебания упругой натянутой струны, где x

является про-

дольной координатой, t - временем, а u(x,t) - смещением участков струны в поперечном на-

правлении.

Коэффициент a0

определяется

материалом струны и силой

ее натяжения:

a

0

=

T

, где

T - сила натяжения, Н

м2

, а ρ- плотность материала, кг

м

.

ρ

Требуется найти решение уравнения колебаний струны (8.5.1) в прямоугольнике

D

=

(x, t)

при начальных условиях (8.3.5) и краевых условиях (8.3.3).

0

< x < l, 0 < t < T

Начальные условия означают, что задано смещение и скорость участков струны в начальный момент времени. Краевые условия означают, что задано смещение начала и конца струны при всех t ≥ 0 ; в частности, если u(0,t)= u(l,t)= 0 , то концы струны являются жестко закрепленными.

Теорема 8.2. Задача (8.5.1), (8.3.3), (8.3.5) имеет единственное решение, которое

непрерывно и дважды непрерывно дифференцируемо в замкнутой области D , если за-

данные функции ϕ0 (x), ϕ1

(x), ψ1 (t), ψ2 (t),

f (x,t) являются достаточно гладкими и, кроме

того, выполнены условия согласования начальных и краевых данных:

ψ

1

(0)= ϕ

0

(0), ψ

2

(0)= ϕ

0

(l), ψ/

(0)= ϕ

(0), ψ/

(l)= ϕ

(l),

1

1

2

1

ψ1// (0)= a02 ϕ0// (0)+ f (0,0), ψ2// (0)= a02 ϕ0// (l)+ f (l,0).

Последние два условия следуют из уравнения (8.5.1) при x = 0, t = 0 и x = l, t = 0 . Для

нахождения приближенных значений решения задачи (8.5.1), (8.3.3), (8.3.5) используем метод конечных разностей. Как и в подразд. 8.4, разобьем отрезок 0 ≤ x ≤ l оси абсцисс на n

частей длины: h = l n , отрезок 0 ≤ t ≤T оси t - на m частей длины τ = T m и построим в об-

ласти D прямоугольную равномерную сетку (см. рисунок на с. 200). В узлах этой сетки точное решение задачи, существование и единственность которого доказана теоремой 8.2, принимает значения u(xi ,t j ).

Получим систему алгебраических уравнений относительно u(xi ,t j ). Для аппрокси-

204

мации производных от решения в узлах

(xi ,t j ), i =1,2,..., n −1, j =1,2,..., m −1

используем

центральные разностные отношения:

∂2u

=

u(xi+1 ,t j

)− 2u(xi ,t j

)+ u(xi−1 , t j )

+ O(h2 ),

(8.5.2)

∂x2

h2

i, j

∂2u

=

u(xi , t j+1 )− 2u(xi ,t j

)+ u(xi , t j−1 )

+ O(τ2 ).

(8.5.3)

∂t 2

i. j

τ2

Соответствующий пятиточечный разностный шаблон показан на рисунке.

Подставляя (8.5.2) и (8.5.3) в уравнение (8.5.1), находим

u(xi ,t j+1 )−2u(xi ,t j )+u(xi ,t j−1 )

−a02

u(xi+1 ,t j

)−2u(xi ,t j )+u(xi−1 ,t j )

+O(τ2 + h2 )

= f (xi ,t j ).

(8.5.4)

τ2

h2

j+1

Отбрасывая

неизвестную

погрешность

O(τ2

+ h2 ), получаем приближенное уравнение

ui, j+1

−2ui, j

+ui, j−1

−a2

ui+1, j −2ui, j

+ui−1, j

= f (x ,t

),

(8.5.5)

τ2

0

h2

i

j

j

где

i =1,2,..., n −1, j = 0,1,..., m −1 .

К

этому

уравнению

следует

присоединить условия,

которые

следуют из

краевых и начальных данных (8.3.3) и (8.3.5):

u0, j

= ψ1 (t j

), un, j

= ψ2 (t j ),

j = 0,1,..., m .

(8.5.6)

j-1

ui,0

= ϕ0 (xi ),

ui,1 − ui,0 = τϕ1 (xi ),

i = 0,1,..., n .

(8.5.7)

Уравнение (8.5.5) представляет собой явную раз-

i-1

i

i+1

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

колебаний струны, так

как с учетом (8.5.6) и (8.5.7) оно позволяет найти решение последовательно на всех слоях по времени t = t j , j =1,2,..., m . Действительно, на первом слое j =1 значения ui,1 легко вычис-

ляются с помощью (8.5.7). Далее перепишем (8.5.5), оставляя в левой части только значение

искомого решения на слое

j +1:

ui, j+1 = 2ui, j

−ui, j−1

−

τ2 a02

(ui+1, j − 2ui, j

+ ui−1, j

)+ τ2

f (xi ,t j ).

(8.5.8)

При j =1 получаем

h2

u

= 2u

− u

− τ 2 a02

(u

− 2u

+ u

)+τ 2 f (x , t

), i =1,2,..., n −1.

i,2

i,1

i,0

h2

i+1,1

i,1

i−1,1

i 1

В этом равенстве ui,1 , ui−1,1 ,

ui+1,1 ,

ui,0 известны в силу (8.5.7), поэтому легко найти ui,2 при

i =1,2,..., n −1. Значения ui,2

в крайних точках i = 0 и i = n

определяются из краевых условий

(8.5.6). Поэтому значения приближенного решения оказываются известными на всем слое t = t2 . Теперь, полагая j = 2 , из (8.5.8) можно найти ui,3 на следующем слое и так далее ана-

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

	Чтобы увеличить точность нахождения приближенного решения на первом слое
	t = t1 , вместо формулы ui,1						− ui,0	= τϕ1 (xi ), которая получена из аппроксимации начального
	условия ϕ	(x	)= ∂u			ui,1	−ui,0	+ O(τ), можно использовать другое соотношение, которое
					=
	1	i	∂t				τ
					t=0

следует из разложения точного решения по формуле Тейлора:

u(xi ,t1 )= u(xi ,0)+ τ	∂u		+	τ2	∂2u		+ O(τ3 ).	(8.5.9)
	∂t			2	∂t 2
		x=x	,t =0				x=x ,t =0
		i					i
			205

В этой формуле первая производная равна ϕ1 (xi ). Для нахождения второй производной используем уравнение (8.5.1):

∂2u

= a2 ∂2 u

+ f (x

,0)= a2 ϕ// (x

)+ f (x

,0).

∂t 2

x=x ,t=0

0 ∂x

2

x=x ,t

=0

i

0 0

i

i

i

i

Подставляя это выражение в (8.5.9) и опуская погрешность аппроксимации, получаем

окончательно

u

= u

+ τϕ

(x

)+

τ2

[a2

ϕ// (x

)+ f (x

,0)].

i,1

i,0

i

1

i

2

0

0

i

Свойства полученного приближенного решения ui, j существенно зависят от того, вы-

полнено или нет условие устойчивости

τa0

≤1.

(8.5.10)

h

Если выполнено, то при достаточно малых шагах сетки ui, j

является хорошим приближени-

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

	1	n				1
				2		2
max		∑	[u(xi ,t j	)−ui, j ]			= O(h2	+ τ2 ).
0≤ j≤m n +1 i=0

Заметим, что, кроме колебаний струны, уравнение (8.5.1) описывает распространение плоских волн сжатия в упругой сплошной среде, например в атмосфере или слое жидкости. Параметр a0 имеет смысл скорости распространения волн, которая может достигать очень

больших значений (несколько километров в секунду в задачах гидроакустики). При этом условие устойчивости (8.5.10) накладывает сильное ограничение на выбор величины τ, что существенно увеличивает необходимое компьютерное время и затраты на проведение расчетов. Чтобы избежать ограничения на выбор τ, связанного с необходимостью выполнения условия (8.5.10), можно использовать неявную разностную схему, которая строится следующим образом.

Запишем разностные аппроксимации производных, входящих в уравнение (8.5.1), в

узле (x

,t

), используя центральное разностное отношение для ∂u2

и разностное отно-

i

j+1

∂x2

шение «со смещением» - для ∂u2

∂t

2 :

u(xi+1 ,t j+1 )− 2u(xi ,t j+1 )+u(xi−1 ,t j+1 )

∂2u

j +1

≈

,

(8.5.11)

∂x2

h2

i, j+1

∂2u

≈

u(xi ,t j+1 )− 2u(xi , t j )+ u(xi ,t j−1 )

,

(8.5.12)

j

∂t 2

i, j+1

τ2

где i =1,2,..., n −1,

j = 0,1,..., m −1 . При этом исполь-

зуется

пятиточечный шаблон,

представленный на

j-1

рисунке. Подставляя (8.5.11) и (8.5.12) в уравнение

i-1

i

i+1

(8.5.1), получаем приближенное уравнение

ui, j+1

− 2ui, j + ui, j−1

− a02

ui+1, j+1 − 2ui, j+1

+ ui−1, j+1

= f (xi ,t j+1 ),

(8.5.13)

τ2

j = 0,1,..., m −1 .

h2

справедливое при i =1,2,..., n −1,

К уравнению (8.5.13) следует присоединить соотношения (8.5.6) и (8.5.7), которые следуют из начальных и краевых условий. Полученную систему алгебраических уравнений

206

можно решить, последовательно определяя				значения приближенного решения на слоях
t = t j , j = 2,3,..., m методом прогонки. Действительно, при					j =1 соотношение (8.5.13) прини-
мает вид
	ui,2 −2ui,1 +ui,0	−a02	ui+1,2	−2ui,2 +ui−1,2	= f (xi ,t2 ),
				h2
	τ 2

где ui,0 , ui,1 известны в силу условий (8.5.7), поэтому получаем трехдиагональную систему уравнений относительно неизвестных ui,2 , которая может быть решена методом, описанным в подразд. 5.4. Значения ui,2 в крайних точках i = 0, i = n , как обычно, находятся из краевых условий (8.5.6). После нахождения приближенного решения на слое t = t2 аналогично можно вычислить значения решения на слое t = t3 и так далее. Применение неявной разностной

схемы (8.5.13) не требует выполнения какого-либо соотношения между шагами h и τ, которое обеспечивало бы устойчивость решения алгебраической системы уравнений.

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

Вкачестве примера рассмотрим уравнение

						∂2u +	∂2u = f (x, y),				(8.6.1)
						∂x2	∂y2		f (x, y)≡ 0 , то уравнение (8.6.1)
которое является частным случаем уравнения (8.3.6). Если
называется уравнением Лапласа.
	Задача	нахождения				решения		уравнения	(8.6.1)	в	области
(x, y)			< y < l
D =	0 < x < l , 0				при краевых условиях (8.3.7) и (8.3.8), заданных на всей гра-
		1		2

нице рассматриваемой области, называется задачей Дирихле.

Теорема 8.3. Задача (8.6.1), (8.3.7), (8.3.8) имеет единственное решение, которое непрерывно и дважды непрерывно дифференцируемо в замкнутой области D , если заданные функции ψ1 (y), ψ2 (y), ψ3 (x), ψ4 (x), f (x, y) являются достаточно гладкими и,

кроме того, выполняются условия согласования:

ψ1 (0)= ψ3

(0), ψ1

(l2 )

= ψ4 (0), ψ2 (0)= ψ3 (l1 ), ψ2 (l2 )= ψ4 (l1 ),

ψ1//

(0)+ ψ3//

(0)= f (0,0),

ψ1// (l2

)+ ψ4// (0)= f (0,l2

),

ψ

//

(0)+ ψ// (l )= f (l

,0),

ψ// (l

2

)+ ψ//

(l

)= f (l ,l

2

).

2

3

1

1

2

4

1

1

Докажем единственность решения. Предположим, что существуют два решения

~

~

u (x, y)

(x, y)

является реше-

и u (x, y) поставленной задачи. Тогда их разность u = u (x, y)

−u

нием уравнения (8.6.1) с нулевой правой частью:

2

2

2

~

2

~

∂

u

+ ∂

u

= ∂

(u

= f (x, y)− f (x, y)= 0

(8.6.2)

(u −u )+ ∂

−u )

∂x2

∂y2

∂x2

∂y2

и обращается в нуль на всем периметре прямоугольника D в силу заданных краевых усло-

вий

u(0, y)= 0,

u(l1 , y)= 0,

0 ≤ y ≤ l2 ,

(8.6.3)

u(x,0)= 0,

u(x,l2 )= 0,

0 ≤ x ≤ l1.

(8.6.4)

Покажем, что решение u(x, y) такой задачи тождественно равно нулю в области D . Умножим обе части уравнения (8.6.2) на u(x, y) и проинтегрируем результат по области D :

	∂	2	u			∂	2	u
					+u
∫∫ u	∂x			2		∂y		2	dxdy = 0 .
D
					207

Левую часть полученного равенства преобразуем так:

∫∫

∂

∂u

∂u 2

∂

∂u

∂u 2

u

−

+

u

−

dxdy = 0 .

(8.6.5)

D

∂x

∂x

∂x

∂y

∂y

∂y

Теперь разобьем двойной интеграл на четыре интеграла и преобразуем первый и третий с помощью формулы Грина

∫	Pdx + Qdy =		∂Q	−	∂P
		∫∫
						dxdy ,
Γ		D	∂x		∂y

где Γ - граница области

D и интегрирование вдоль контура Γ происходит в направлении

против часовой стрелки. Тогда из (8.6.5) получаем

∂u

∂u

∫∫

∂u 2

∂u 2

∫

∂y

dx −u

∂x

−

∂x

∂y

dxdy = 0 .

(8.6.6)

u

dy

+

Γ

D

В силу условий (8.6.3) и (8.6.4) функция u равна нулю на контуре Γ , поэтому в (8.6.6) остается только интеграл по области D :

∫∫	∂u 2			∂u 2
		∂x		∂y	dxdy = 0 .
			+
D

Из последнего равенства следует ∂u ∂x ≡ ∂u ∂y ≡ 0 . Поэтому u(x, y)= const , и так как на боко-

вых сторонах эта функция обращается в нуль, то u(x, y)≡ 0 в D . Вспоминая определение u , находим u(x, y)≡ u~(x, y), что и требовалось доказать.

Для нахождения приближенного решения задачи (8.6.1), (8.3.7) – (8.3.8) используем метод конечных разностей. Введем такую же расчетную сетку в области D , как и в предыдущих параграфах. Для аппроксимации производных от решения во внутренних узлах (xi , y j ), i =1,2,..., n −1, j =1,2,..., m −1 области D используем формулы:

∂2u

=

u(xi+1 , y j )− 2u(xi , y j )+ u(xi−1 , y j )

+ O(h2 ),

∂x2

i, j

h2

1

1

∂2u

=

u(xi , y j+1 )− 2u(xi , y j )+ u(xi , y j−1 )

+ O(h2 ).

∂y2

h2

2

i, j

2

Подставляя (8.6.7) и (8.6.8) в уравнение (8.6.1), находим

)+ u(xi , y j−1 )

u(xi+1 , y j )− 2u(xi

, y j )+ u(xi−1 , y j )

+

u(xi , y j+1 )−

2u(xi , y j

+

.

2

= f (x

)

2

h1

+ O(h2 + h2 )

, y

j

h2

1

2

i

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

(8.6.7)

(8.6.8)

(8.6.9)

	j+1			Отбрасывая			неизвестную		погрешность
		O(h12 + h22 ), получаем приближенное уравнение
			ui+1, j	−2ui, j +ui−1, j	+	ui, j+1 −2ui, j +ui, j−1		= f (x , y			), (8.6.10)
				h2			h2		i	j
	j			1			2
		где		i =1,2,..., n −1,			j =1,2,..., m −1. К этому уравнению
		следует присоединить условия, которые вытекают из
	j-1	(8.6.7) и (8.6.8):

i-1	i	i+1	u0, j = ψ1 (y j ), un, j = ψ2 (y j ), j = 0,1,..., m, (8.6.11)
			208

ui,0 = ψ3 (xi ), ui,m = ψ4 (xi ), i = 0,1,..., n .

(8.6.12)

Система алгебраических уравнений (8.6.10) содержит (n −1)(m −1) неизвестных и

столько же уравнений. В отличие от подразд. 8.4 и 8.5 здесь не удается получить решение путем последовательного его нахождения на слоях x = const или y = const . Вместе с тем

можно применить другие методы решения линейных систем уравнений, например метод простых итераций. Для простоты рассмотрим его реализацию в частном случае, когда h1 = h2 = h . Тогда соотношения (8.6.10) принимают вид

u

i, j+1

−2u

i, j

+u

i, j−1

+u

i+1, j

−2u

i, j

+u

i−1, j

= h2 f (x , y

j

),

i

то есть

1

[ui, j+1

+ ui−1, j − h2 f (xi , y j )].

ui, j

=

+ ui, j−1

+ ui+1, j

(8.6.13)

4

Выбирая некоторое нулевое приближение ui(,0j) , удовлетворяющее краевым условиям (8.6.11) и (8.6.12), находим следующие приближения согласно методу простых итераций:

	ui(,kj+1) =	1	[ui(,kj)+1 +ui(,kj)−1 +ui(+k1,)	j +ui(−k1,)	j −h2	f (xi , y j )].	(8.6.14)
		4

Система уравнений (8.6.13) записана в том же самом виде, как в подразд. 5.5, а именно, с использованием двух индексов в обозначении неизвестных ui, j . Можно показать, что

справедлива теорема, аналогичная теореме 5.4, поэтому последовательность приближений решения (8.6.14) системы алгебраических уравнений сходится к точному решению этой системы:

	ui(,kj) −ui, j		∞	= max		ui(,kj) −ui, j		→ 0 при k →∞.
				i, j

Аналогичным образом может быть найдено решение системы (8.6.10) в случае раз-

ных шагов сетки h1 и h2 . Для полученного решения ui, j системы уравнений (8.6.10) доказана следующая оценка:

	ui, j −u(xi , y j )				∞	= O(h12	+ h22 ),

из которой следует сходимость ui, j к решению u(xi , y j ) задачи Дирихле для уравнения

(8.6.1) при неограниченном измельчении сетки: h1 →0 , h2 →0 .

Если область D не является прямоугольником, то задача Дирихле состоит в нахождении функции u(x, y), являющейся решением уравнения (8.6.1) в области D и принимающей

заданные значения u Γ = ψ(x, y) на ее границе Γ . В этом случае алгоритм решения задачи

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

производных оказываются несправедливыми и возникает

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

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

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

209