Ращиков В.И

	ВАРИАНТЫ ЗАДАНИЙ
0	< αi <βi < lx , 0 <μi < νi < ly (i = 0 ÷4);
0	< pi , qi , ri , si ≤ 3 (i = 0 ÷4);

γi = 0.1÷10 (i =1÷5), ci = 0.1÷2 (i =1÷5).

ВАРИАНТЫ НАЧАЛЬНЫХ УСЛОВИЙ

1)u0 (x, y) =1;

2)u0 (x, y) = c0η(x −α0 )η(β0 − x)η( y −μ0 )η(ν0 − y);

3)u0 (x, y)

4)u0 (x, y)

5)u0 (x, y)

					p0
= c0			x			1
	lx
					p0
= c0			x				1
	lx
= c	x				p0
						1
0
	lx

x

q0

r0

y

s0

−

y

1

−

;

lx

ly

ly

x

q

y

s0

0

−

sinr0

π

;

l

x

l

y

y

q0

x

s0

−

sinr0

π

.

ly

lx

ВАРИАНТЫ КРАЕВЫХ УСЛОВИЙ

(i=1,2)

1)ψi ( y,t) = η( y −μi )η(νi − y)e−γit ;

2)ψi ( y,t) = ciη( y −μi )η(νi − y)(1 − e−γit );

3)ψi ( y,t)

4)ψi ( y,t)

5)ψi ( y,t)

y pi

y qi

= ci

1

−

(1 −e−γit );

ly

ly

y pi

y qi

= ci

1

−

e−γit ;

ly

ly

y

si

= c sinri

π

e−γit ;

i

l

y

111

			y		si		(1 −e−γit ).
	6) ψi ( y,t) = ci sinri	π
		l		y

ВАРИАНТЫ КРАЕВЫХ УСЛОВИЙ

1)ψi (x,t)

2)ψi (x,t)

3)ψi (x,t)

4)ψi (x,t)

5)ψi (x,t)

6)ψi (x,t)

=ciη(x −αi )η(βi −

=ciη(x −αi )η(βi −

		x	pi				x	qi
	= ci				1	−
							lx
	lx

(i=3,4)

x)e−γit ;

x)(1−e−γit );

e−γit ;

pi

x

qi

= ci

x

1

−

(1−e−γit );

lx

lx

x

si

= ci sinri

π

e−γit ;

lx

x

si

= ci sinri

π

(1

−e−γit ).

lx

ВАРИАНТЫ ИСТОЧНИКОВ ТЕПЛА

1)f (x, y,t) =c5u0 (x, y)e−γ5 ,

2)f (x, y,t) =c5u0 (x, y)(1−e−γ5 ),

где u0(x,y) – одно из возможных начальных условий, указанных выше.

	ПРОГРАММИРОВАНИЕ
Функции f, u0,	ψ целесообразно оформить в виде подпро-
грамм – функций.	Для хранении u ji и u ji можно отвести два дву-

мерных массива размерности (N+1)(M+1), хотя при более искусном программировании достаточно одного двумерного массива

112

размерности (N+1)(M+1) и одного - двух одномерных рабочих массивов размерности N+1. Число пространственных шагов N, M должно составлять несколько десятков, временной шаг τ следует выбирать равным (0.2÷0.8)τmax, где τmax дается формулой (16.6). Число временных шагов Nt следует выбирать таким, чтобы процесс успел установиться, полагая, например,

γ5 Nt τ 10 ÷100.

Для наблюдения за процессом целесообразно через несколько десятков шагов τ выводить на экран распределение температуры u(x,t) в исследуемой области. Как альтернативу можно вывести распределение u(x,y,t) вдоль средней линии y=ly/2. Для этого надо описать вспомогательный одномерный массив v , имеющий N +1 элементов, перед выдачей пересылать в него значения u(j,M/2) (j = 0,1,...N).

Для проверки программы рекомендуется предварительно решить тестовую задачу, полагая

a2 =1, lx = ly =1, ψ1 = ψ2 = ψ3 = ψ4 = 0;

u0 (x, y) =1, f ( x, y,t) = 32[x(1 − x) + y(1 − y)].

При этих данных установившееся распределение температуры (при больших t ) имеет вид

u(x, y,t) =16x(1− x) y(1− y) ,

а вдоль средней линии y = 0.5 решение имеет вид параболы

u(x,0.5,t) 4x(1− x)

с вершиной u=1 в точке x = 0.5.

Блок-схема программы представлена на рис. 16.1.

В блоке 2 вводятся параметры αi ,βi , μi ,νi , pi , qi , ri , si , γi , ci . В цикле 3-4-5 задаётся начальное распределение потенциала в области. Цикл 6-13 является основным в программе, в нём с шагом τ осуществляется изменение времени t. В этот цикл вложены циклы 7-8-9 и 10-11-12, в которых соответственно задаются граничные условия и вычисляется распределение потенциала внутри области в текущий момент времени. В блоке 14 выводятся полученные результаты.

113

1 Начало

Началь-

2ные

данные

Цикл

3по j=1,N-1 по i=1,M-1

u0ji

1

4

5				Цикл
				по j,i
6		Цикл
		по k=1,kmax
				Цикл
7		по j=1,N-1
		по i=1,M-1
8
		u0i ,uNi ,u j0 ,u jM
9				Цикл
				по j,i

Цикл

10по j=1,N-1 по i=1,M-1

11

uji

12Цикл по j,i

13Цикл по k

Резуль- 14 таты

15 Конец

Рис. 16.1. Блок-схема программы решения уравнения теплопроводности в прямоугольнике

114

СОДЕРЖАНИЕ ОТЧЁТА

Отчет должен содержать:

•формулы и параметры для конкретного варианта;

•текст программы;

•результаты решения и графики.

КОНТРОЛЬНЫЕ ВОПРОСЫ

1.Каков порядок погрешности разностной схемы (16.4)?

2.Получите условие устойчивости (16.6 ) используя спектральный критерий.

3.Как строится решение разностной схемы (16.4)?

4.Изобразите расположение узлов ("шаблон"), на котором построена разностная схема (16.4).

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

(16.5)?

115

Задание № 17

ЧИСЛЕННОЕ РЕШЕНИЕ ОДНОМЕРНОГО ВОЛНОВОГО УРАВНЕНИЯ

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

ПОСТАНОВКА ЗАДАЧИ

Распространение волн в одномерном приближении описывается уравнением

	d 2u	= a2	d 2u	+ f
	dt2		dx2
	(0 ≤ x ≤ l,0 ≤ t ≤T )

с начальными условиями

u(x,0) = μ1 (x) (0 < x < l),
u ′(x,0) = μ		2	(x) (0 < x < l),
t
и краевыми условиями
	u(0,t) = μ3 (t),
	u(l,t) = μ4 (t)
	(0 ≤t ≤T ).
Здесь х - координата,	t - время, a - фазовая скорость волны;

функция f(x,t) описывает влияние среды, т.е. поглощение или возбуждение волн в среде.

Для решения воспользуемся неявной трехслойной разностной схемой

u

j

− 2u

+u

j

u

j

+u

j

j

= Λ

+ f j

( j =1, 2,..., N −1).

(17.1)

τ

2

2

116

Здесь

uj = u(xj ,t), uj = u(xj ,t +τ), uj = u(xj ,t − τ); f j = f (xj ,t), xj = jh ( j = 0,1,..., N ), h = l / N;

u0 =μ3,uN =μ4 ,

а оператор

Λu j = a2 (u j+1 −2u j +u j−1 ). h2

Разностную схему (17.1) можно преобразовать к виду

u j+1 −σu j

+u j−1 = g j ( j =1,2,..., N −1),

g

= −2bu

−(u

−σu

+u

) −

2h2

f

,

j

j

j+1

j

j−1

a2

j

σ = 2 +b,

b =

2h2

.

a2

τ2

Полученное уравнение решается методом прогонки, как и в выполненных ранее работах. Прогонка строится в два этапа.

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

ξj+1 =1 / (σ−ξj ), ηj+1 = (ηj − g j )ξj +1 ( j =1, 2,..., N −1),

ξ1 = 0, η1 =μ3.

Второй этап - вычисление значений uj , (обратный ход прогон-

ки):

u j =ξj+1u j+1 +ηj+1 ( j = N −1, N −2,...,1), uN = μ4 .

Вместо схемы (17.1) можно также воспользоваться более простой явной трехслойной разностной схемой

u j − 2u j +u j		= a2	u j +1 − 2u j +u j−1	+ f		, ( j =1, 2,..., N −1),
τ2			h2		j
						(17.2)
u0 = μ3 ,	uN = μ4 ,

которая преобразуется в расчетную формулу

117

u j = 2(1 −c2 )u j −u j + c2 (u j+1 +u j −1 ) + τ2 f j , c = (τa / h),

где с - число Куранта.

Данная схема устойчива при выполнении условия Куранта: c <1.

Для начала вычислений необходимо, помимо u(x,0), располагать значением u( x, τ ). Поскольку решение u(x, t ) на каждом временном шаге находится с погрешностью O(τ3) значение u( x, t) также необходимо иметь с той же погрешностью O( τ3). Для нахождения воспользуемся рядом Тейлора:

du

d 2u

τ2

3

u(x, τ) =u(x,0) +

dt

τ +

dt2

+O(τ

).

2

t =0

t=o

Из исходного уравнения следует, что

d 2u

= a

2

d 2u

+ f (x,0),

dt2

t =0

dx2

t

=0

и из начальных условий

du

= μ

2

(x),u(x,0) = μ (x),

dt

t =0

1

d 2μ

d 2u

=

t =0

1

.

dx2

dx2

Окончательно имеем:

d 2μ (x)

u(x

, τ) = μ (x

) + τμ

(x

) +

τ2

a2

+ f (x

,0)

j

2

j

1

j

1

2

dx2

x=x j

j

( j =1, 2,...N −1),

u(0, τ) = μ3 (τ),

u(l, τ) = μ4 (τ).

Затем из уравнения (17.1) находим значения u(x,t) в моменты времени t=2τ, 3τ и т.д.

В дальнейшем ограничимся случаем, когда f = 0, т.е. среда (внутренние источники) отсутствует. Таким образом, будут рассматриваться лишь свободные колебания и вынужденные колебания под действием краевой возбуждающей силы. Характерными временами задачи являются:

T =	l	, τ		=	h	,
			0
0	a				a
		118

где T0 - время прохождения волны вдоль системы, τ0 - время прохождения волной шага h. Можно показать, что диапазоны длин волн λ и частот ν, которые могут наблюдаться в рассматриваемой системе, составляют:

2h ≤λ ≤ 2l, 0.5 / T0 ≤ ν ≤ 0.5 / τ0 .

	ВАРИАНТЫ ЗАДАНИЙ
	ci =1 ÷10(i =1÷4).
	ВАРИАНТЫ ФУНКЦИЙ μ1(x)
1)	μ1 = c1 (x / l)(1− x / l);
2)	μ1 = c1 sin(πnx / l) (n =1÷3);

3)μ1 = 0;

4)μ1 = c1 (x / l)2 (1 − x / l);

5) μ1 = c1 cos(πnx / l) (n =1 ÷3).

ВАРИАНТЫ ФУНКЦИЙ μ2(x)

1) μ2

= 0;

c (x [α,β]),(0 < α <β <l),

2) μ2

2

=

0(x [α,β];

3) μ2

= c2 (x / l)(1− x / l);

4) μ2

= c2 sin(πnx / l) (m =1 ÷3);

5) μ2

= c2 x / l

(0 < x ≤ l / 2),

c2 (1 − x / l) (l / 2 < x < l).

ВАРИАНТЫ ФУНКЦИЙ μi(x) (i = 3,4)

(рi = 0.2÷2, νi

=1÷10).

1) μ

i

= c e− pit /T0

;

2) μ

i

= c (1 −e− pit /T0

);

i

i

3) μ

= c e− pit /T0

cos

2πt

ν

;

4) μ

= c (1 −e− pit /T0

)cos

2πt

ν

;

i

i

i

i

i

i

T0

T0

119

5) μi = ci (t / νiT0 )(1 −t / νiT0 )			(0 ≤ t ≤ νiT0 ),
0 (t > νiT0 ).
ПРОГРАММИРОВАНИЕ
Функции μi (i = 1 до	4) удобно оформить в виде подпрограмм -
функций. Для хранения	uj ,u j ,u j		целесообразно отвести три одно-
мерных массива по N+1 элементов в каждом, и после завершения
очередного временного шага		с помощью оператора цикла перепи-
сывать uj на место uj	, а uj	на местоuj .

Значения l и a можно положить равными 1, а число пространственных шагов N, временной шаг τ и длительность счета Т рекомендуется выбирать из условий:

N=25100, τ = (0.2- 0.8) τ0 , Т= (5-10) Т0 .

Для наблюдения за ходом вычислений через каждые M =20 ÷ 100 шагов следует выдавать график uj .

Для проверки правильности программы рекомендуется предварительно решить тестовую задачу, полагая

μ = μ	3	= μ	4	= 0,	μ	2	(x) = πa sin πx .
1							l	l

Решение этой задачи должно быть близко к стоячей волне вида

u(x,t) = sin πlx sin πlat .

Блок-схема программы для явной схемы представлена на рис.

17.1. В блоке 2 вводятся параметры μi , νi , pi , ci ,kmax , M . В цикле 3- 4-5 задаются начальные распределения в моменты времени 0 и τ. Цикл 6-15 является основным в программе, в нём с шагом τ осуществляется изменение времени t. В этот цикл вложены циклы 7-8-9 и

12-13-14, в которых соответственно вычисляется значения uj в те-

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

ствляется в блоке 11, если это необходимо. В блоке 16 выводятся полученные результаты.

120