Ращиков В.И

В качестве T(x,t) принять

T (x1 , x2 , x3 ,t) = ρ(x1 + g1t, x2 + g2t, x3 + g3t),

где ρ(х) — одна из функции предыдущего задания, g1=0.l ÷0.9

(i =1,2,3).

3. Вычислить объем |V| тела, ограниченного шестимерным эллипсоидом

	6	xi
	Г = ∑(		)2 −1 = 0	(ci = 0.1÷2).
		ci
	i=1

Вычисление выполнить методом Монте-Карло по формуле

V = ∫ dx = ∫ dx1dx2dx3dx4dx5dx6 .

(V ) (V )

ПРОГРАММИРОВАНИЕ

Вычисление интеграла проведем для нескольких значений N, заданных отдельным массивом N (L) в головной программе. Соответственно, в программе следует организовать выдачи результатов по достижении числа случайных чисел очередного значения N из массива. Это даст возможность наблюдать за изменением результатов и сходимостью интегрирования с ростом N.

Блок-схема программы представлена на рис. 7.2. Основу программы составляет цикл (блоки 3-10) по l от l до L, где L — заданное число вариантов с различным количеством случайных чисел Nl, для которых осуществляется выдача результатов. В блоке 4 происходит обращение к датчику случайных чисел для вычисления ξ. На рисунке j — текущее число случайных точек, М —

число случайных точек в области V, I è V — оценки интеграла I и

объема области V.

В качестве теста необходимо вычислить объём эллипсоида в трёхмерной области в соответствии с пунктом 3 задания. Известно аналитическое решение

V = 43 πc1c2c3.

51

1 Начало

Начальные 2 данные

M=0,j=0,S=0

3

4

6

7

Цикл по l=1,L

ξj

ξj, j=j+1

5

Нет

Г(ξ)≤0

Да

M=M+1,S=S+f(ξ)

Да

j<Nl

Нет

8V = MW / Nl , I ≈ WS / Nl

9I , Nl , M , V

10Цикл по l

Графики 11 Результаты

12 Конец

Рис. 7.2. Блок-схема программы вычисления интеграла методом Монте-Карло

52

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

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

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

•значения интеграла I для нескольких значений N, выбран-

ных в диапазоне N~ 102- 104.

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

1.Какие численные методы называются методами МонтеКар-

ло?

2.Как вычисляются кратные интегралы методом Монте-Карло?

3.Как получена формула (7.3)?

4.Каков порядок погрешности интегрирования методом МонтеКарло?

53

Задание № 8

РЕШЕНИЕ СИСТЕМ ЛИНЕЙНЫХ АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ

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

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

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

a x + a x +... + a x = b ,
	11 1	12		2	1n n		1
a21 x1 + a22 x2 +... + a2n xn = b2 ,
..............................................
a	x + a	n2	x	2	+... + a	x	= b ,
	n1 1					nn n	n
которую можно короче записать в матричном виде:
		Ax = b,					(8.1)

где А — квадратная матрица размера n × n, x и b — n-мерные векторы:

A = (aij; i, j = 1, 2, …, n), x = (xi; i = 1, 2, …, n), b = (bi; i = 1, 2, …, n).

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

Прямой метод

Будем полагать, что матрица А невырожденная, т. е. детерминант A ≠ 0 и, следовательно, решение системы (8.1) существует и единственно.

Преобразуем первое уравнение системы (8.1) к виду

x1 + a12(1) +... + a1(1)n xn = b1(1)

54

a(1)	=	a1 j	( j = 2, 3, ...n), b(1)	=	b
					1	.
1 j		a11	1		a11

Умножая затем это уравнение на – ai1 и складывая с i-м уравнением (i = 2, 3, …, n), исключим х1 из всех уравнений, получая систему (n – 1) уравнений. Далее аналогично исключим x2, x3 и т. д. Этот процесс исключения называется прямым ходом метода Гаусса. Наконец, остается одно уравнение, из которого находим xn, а затем последовательно остальные неизвестные xn–1, xn–2, …, x1; этот этап решения называют обратным ходом.

Одним из наиболее эффективных прямых методов является используемый в данной работе метод факторизации, называемый также схемой Халецкого (Холесского) или компактной схемой Гаусса. В этой схеме матрица А представляется в факторизованном виде

A = LU,	(8.2)

где L (lower) — нижняя треугольная матрица, U (upper) — верхняя треугольная матрица с единицами вдоль главной диагонали:

	l	0	...	0			1	u
	11	l22	...	0				12
L =	l21					, U =	0	1
		...	...	...
	...						... ...
	ln1	ln2	...	lnn			0	0

... u1n

... u2n .

... ...

... 1

Элементы матриц L и U находятся из матричного уравнения (8.1). Матрица U может быть получена из (8.1) как результат прямого хода метода Гаусса, а матрица L хранит коэффициенты, на которые умножались уравнения в процессе исключения неизвестных в прямом ходе.

После этого исходная система (8.2) сводится к последовательному решению двух систем

	Ly = b,	(8.3)
		(8.4)
	Ux = y

с треугольными матрицами.

Введем расширенные матрицы, дополняя А и U еще одним столбцом:

A + b, ai, n+1 = bi, U + y, ui, n+1 = yi.

(8.5)

55

Решение выполняется по формулам:

		j −1
lmj	= amj	− ∑ lmk ukj ( j = 1, 2, ..., n; m = j, j +1, ..., n);	(8.6)
		k =1
	i−1
uim = aim − ∑lik ukm / lii (i =1, 2, ..., n; m = i +1, i + 2, ..., n			+1);
	k =1		(8.7)
		n
		xi = ui ,n +1 − ∑ uik xk (i = n, n −1, ..., 1).	(8.8)
		k =i+1

Формулы (8.6), (8.7) описывают факторизацию матрицы А и решение уравнения (8.3) (прямой ход метода Гаусса), формула (8.8) — решение уравнения (8.4) (обратный ход метода Гаусса).

Решение этим методом встречается с трудностями, если в формуле (8.7) lii 0. Число арифметических операций в методах Гаус-

са имеет порядок O(n3).

Контроль правильности решения осуществляется по величине невязки

ξ = Ax −b,

норма которой									(8.9)
	ξ				= max		( Ax)i −bi
					i

должна быть близка к 0. (В формуле выше (Ax)i – i-я компонента вектора Ах, где х—полученное решение.)

Итерационный метод

Преобразуем исходную систему уравнений к виду

		x1 = (b1 −a12 x2 +... +a1n xn ) / a11
	x2 = (b2 − a21x1 + +a2n xn ) / a22
..............................................
x	n	= (b −a	x +... + a	nn−1	x	n−1	) / a	nn	.
		n	n1 1

Если теперь в правую часть системы в качестве вектора x взять произвольный начальный вектор

x(0) = (x1(0) , x2(0) ,......, xn(0) ) ,

56

то в	левой		части получим новое значение вектора
x(1) = (x(1) , x(1)		,......, x(1) ) .
1	2		n

Таким образом мы получаем итерационный метод Якоби, имеющий в компонентах вид

xi(k ) = −∑aij x(jk −1) / aii +bi / aii (i =1, 2, ..., n).

i≠ j

Метод Якоби сходится, если матрица A обладает диагональным преобладанием, т. е.

aii > ∑ aij (i =1, 2, ..., n)

j≠i

(все диагональные элементы aii > 0, так как матрица А положительно определенная).

Представим матрицу А в виде

A = L + D + U,

(8.10)

где L — нижняя, a U — верхняя треугольные матрицы с нулями вдоль главных диагоналей, D — диагональная матрица:

0	0	...	0		a		0	...	0
						11
L = a21	0	...	0		, D =	0	a22	...	0		,
...	...	...	...		...		...	...	...
	an2	...	0			0	0	...
an1									ann

	0	a	...	a
		12		1n
U =	0	0	...	a2n .
... ...			...	...
	0	0	...	0

Итерационный метод Зейделя имеет вид

x(k) = –D–1(Lx(k) + Ux(k–1) – b)
или в компонентах
xi(k ) = −aii−1	i−1	n
	∑ aij x(jk ) + ∑ aij x(jk −1) −bi .
	j=1	j=i+1

(8.11)

(8.12)

57

Как мы видим, в этом методе в отличие от метода Якоби в правой части используются значения x(jk ) , уже вычисленные в ходе

данной итерации.

Здесь индекс k – номер итерации. Начальное приближение x(0) в итерационных методах может быть выбрано произвольно, например, можно принять x(0) = 0 или x(0) = b. Итерации проводятся до выполнения одного из условий:

max	xi( k ) − xi( k −1)	< δ;		ξ( k )	< δ; k ≥ kmax ,
i

где δ — заданная малая величина, kmax — предельное число итераций. Число арифметических операций в итерационных методах имеет порядок O(n2k), где k — число выполненных итераций.

ВАРИАНТЫ ЗАДАНИЙ

n = 5 20;

p, q =1 4;

r,t = 0.1 1.2.

1) a

= 5ip/ 2 ,

a

= ±0.01(i p / 2 ± jq / 2 )

(

i

≠ j

) , b = 4.5ip/ 2

;

ii

ij

i

2) a

=15ip ,

a

= ±i− p / 2 j−q / 2 (i ≠ j),

b = 3.5i p ;

ii

ij

i

3) a =10i p/ 2

, a = ±10−3 (i / jr )1/ q (i ≠ j), b = 9ip/ 2 ;

ii

ij

i

4) a = 8ip/3,

a = ±10−2 (i p

± jq )−1 (i ≠ j), b = 7ip/3 ;

ii

ij

i

5) a = 9ip/3,

a = ±10−2 (i p / 2 ± jq / 2 )−1 (i ≠ j), b = 8ip/ 3

;

ii

ij

i

6) a

= 8.5ip/3 , a

= ±10−3 exp(−

r

)(i ≠

j), b = 7.5it / 2

;

ii

ij

r

i ± j

t

i

7) a

=12ir ,

= ±0.1(ij)−

b =10.5it / 2 ;

a

ij

2

(i ≠ j),

ii

i

8) a

=14ir / 3

,

a

=10−2 sh

r

(i

≠ j),

b = 9.5it / 2 ;

ii

ij

i ± j

t

i

9) a

=11ir / 2 ,

a

= ±10−2 sin

πr / 2

(i ≠

j),

b = 8.5tt / 3.

ii

ij

i ± j

t

i

58

ПРОГРАММИРОВАНИЕ

Компактная схема Гаусса. Для экономии памяти компьютера треугольные матрицы L, U и вектор у можно объединить в единую матрицу с элементами

l11		u12	u13		u1n		u1,n+1
l		l	u	23	u	2n	u		2,n+1
	12	22
l		l	l		u		u
	31	32	33		3n			3,n+1
	l	l	l		l		u	n,n+1
	1n	2n		3n	nn

(поэтому этот метод и называют компактной схемой Гаусса).

При таком объединении главная диагональ матрицы U в памяти компьютера исчезает, однако равенство единице всех ее элементов уже учтено в формуле обратного хода (8.9) тем, что коэффициент

перед ui,n+1=yi равен 1. Удобнее, однако, сохранить обозначения в формулах (8.2) — (8.8), но совместить матрицы L и U в памяти

компьютера, пользуясь оператором эквивалентности.

При составлении программы некоторые усложнения возникают из-за того, что суммирование в формулах (8.6) — (8.8) допустимо лишь тогда, когда верхний предел суммы не меньше нижнего.

Кроме того, вычисление решения х в (8.8) можно организовать так, чтобы переменная цикла i не убывала, а возрастала. Составленная с учетом этих особенностей блок-схема программы приведена на рис. 8.1.

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

В блоках 3-9 вычисляются матрица А и правая часть b системы уравнений, блоки 10-29 составляют прямой ход, блоки 31-38 — обратный ход, в блоке 39 вычисляется норма невязки.

59

1	Начало			Цикл			Цикл по		xn	= un, na	31
				по m= j,n		12	m= ia,na	20			32
	Начальные								Цикл
2	Данные								по i=1,n-1
	na=n+1			w = amj		13	w = aim	21			33
									n2	= n −i,
	Цикл
					14		22		k0	= na −i,
3 по i = 1,n					Да		i1 <1	Да	w = 0
				j	<1						34
				1
4	a , a			Нет	15		Нет	23	Цикл
	ii	ina							по k=k0,n
				Цикл			Цикл
	Цикл			по k=1,j1			по k=1,i1				35
5 по j = 1,n					16			24	w = w + un2k xk
6	aij (i ≠ j)			w = w −lmk ukj			w = w −lik ukm			Цикл
				Цикл			Цикл			по k
7	Цикл			по k			по k			36
	по j				17		25		xn2	= un2 na	− w
	Цикл							26			37
										Цикл
8
	по i						ulm = w / lii
				lmj = w		18				по i
											38
9	A,b,n			Цикл		19	Цикл	27		x, ξ c
				по m			по m				39
	i=1,j=1						28Да			Конец
	10								L,U	40
							i ≥ m
	i1 = i −1,ia		= i +1,				Нет		30
	j1 = j −1						i = i +1, j = j +1 29
	11
				Рис. 8.1. Блок-схема метода факторизации

60