Практикум j практическое примщенше численных методов

¥

К

ill!

0,

0.

liliii

0.

illlli

0.

vMOV-V

-VvM-

1

«1

liflllf

1.

illll

0.

llll

0

o

^:/М4\

J

тт; ^

^:'Лт

-1.

111111

M^:lf

mmm illll

0

ills

9ШШ

4 о

h:'%.

¥

II

-1

~1

-1

llll

Ml

ШШШ illll

1

¥io

M0‘^; '

V 16

)

-1

-1

“!

ill» lllllll

-1.

Illll

Mm-

VO;

■ 32

1:

-1.

Illll

да

111

Mil •

64

1 л

ШМ - Ш}-Ь:

-1:

ууМ

-1

llll

*1.

Ш'.'т-

^:fM'

128

и.. ■ Ш-

■'■i.S

1.

<1

W^:i

1

:;Й

Ml

ШИ

Рассматриваемая модельная задача имеет точное решение. В частности, тре­угольная матрица U есть
uij = < 2^J, j = n, г = 1,2,...,п, j = l,2,...,n.
I 0, г 7^ j, j ф n,
в LCZ-разложении матрицы A.
Упражнение 4.2 Напишите программу, реализующую разлоэюение Холец­кого А = LL* для симметричной положительно определенной матрицы А и вычисляющей определитель матрицы на основе этого разложения. Найди­те разлоэюение Холецкого и определитель матрицы Гильберта, для которой
^aij = т—4—т, г = 1,2,n, j = 1,2, ...,п г +j - 1
при п = 4.
Функция decCholO в модуле chol, реализует разложение Холецкого (вычис­ляется нижняя треугольная матрица L), функция detCholO вычисляет опре­делитель матрицы А.
import numpy as пр
def decChol(A):
м и и
Returns the decompositon Choleski for matrix A.
и и и
n = len(A)
L = np.copy(A) for j in range(n): try:
L[j,j] = np.sqrt(L[j,j] - np.dot(L[j ,0: j] , L[j,0:j])) except ValueError:
print ’Matrix is not positive definite’ for i in range(j+1,n):
L[i,j] = (L[i,j] - np.dot(L[i,0:j] , L[j,0:j])) / L[j,j] for j in range(l,n):
L[0: j , j] = 0. return L def detChol(A):
и и и
Returns the determinant of the matrix A.
и и и
n = len(A)
L = decChol(A) det = 1.
for i in range(n): det *= L[i,i]**2 return det
Разложение матрицы Гильберта и вычисление ее определителя проводится в следующей программе.
import numpy as np
from chol import decChol, detChol
n = 4
A = np.zeros((n, n), ’float’) for i in range(0,n):
for j in range(0,n):
A[i,j] = 1./(i+j+1) print ’A:\n’, A L = decChol(A) print ’L:\n’, L det = detChol(A) print ’det = ’, det
Определитель матрицы Гильберта при возрастании п стремиться к нулю.
Упражнение 4.3 Напишите программу, которая решает систему линей- них уравнений для трехдиагональной (ац = 0 при \i — j\ > I) матрицы пхп на основе LU-разлоясения. Найдите решение уравнения с
ап 2, ttj — Q-i.t+i — ~ 1
при постоянной правой части /* = 2/i², h = (n + l)^-1, i — 1,2, ...,n и срав­ните его с точным решением у_{ = ih( 1 — г/i), г — 1,2,..,, п. Рассчитайте определитель матрицы и сравните его значение с точным (п+ 1).
Трехдиагональная матрица А задается тремя диагоналями:
tti ttij, bi — &г,г+ 1j Cj — 0-г,г—l¹
В модуле Iu3 функция decLU3() предназначена для LU-разложение трехдиа- гоналы-юй матрицы А. Результат записан в виде трех, диагоналей, причем
di ⁼ 1ц, Щ ⁼ Щ,1-j_i, Ц — ^г,г— 1 •
Для решения системы уравнений используется функция solveLU3().
import numpy as np def decLU3(a, b, c):
* и и и
Input of tridiagonal matrix A: a[i] = A[i,i] , b[i] = A[i,i+i], c[i] = A[i,i-1]
Returns the decompositon LU for tridiagonal matrix: d[i] = L[i,i] u[i] = U[i,i+1]
1 [i] = L [i, i-1]
и и и
n = len(a) d = np.copy(a) u = np.copy(b)
1 = пр.сору(с) for i in ranged, n): al = 1 [i] / d[i-l] d[i] = d[i] - al*u[i-l]
1 [i] = al return d, u, 1 def solveLU3(a, b, c, f):
к и и
Solve the linear system Ax = b with tridiagonal matrix: a[i] = A [i, i], b[i] = A[i,i+1] , c [i] = A[i,i-1]
и и и
n = len(a)

• 
  LU decomposition
  

d, u, 1 = decLU3(a, b, c) x = np.copy(f)

• 
  forward substitution process for i in ranged, n):
  

x[i] = x[i] - 1 [i]*x[i-l]

• 
  back substitution process x[n-l] = x[n-l] / d[n-l] for i in range(n-2,-1,-1):
  

.x [i] = (x[i] - u[i]*x[i+l]) / d[i] return ^?x
Решение системы уравнений, вычисление определителя и сравнение с точны­ми данными проводится следующей программой.
import numpy as np
from lu3 import solveLU3, decLU3
n = 50
h = 1. / (n + 1) a = 2*np.ones((n), ’float’) b = - np.ones((n), ’float’) c = - np.ones((n), ’float’) f = h**2*2*np.ones((n), ’float’)

• 
  solution
  

x = solveLU3(a, b, c, f)

• 
  test function
  

у = np.zeros((n), ’float’)

• 
  LU decomposition
  

d, u, 1 = decLU3(a, b, c) er = 0.
det = 1.
for i in range(n):
у [i] = (i+l)*h*(l-(i+l)*h) er = er + (y [i] -x[i] )**2*h det = det*d[i] er = np.sqrt(er) print ’error = ’, er
print ’determinant calculated =’, det, ’\nexact =’, n+1
Ш®Я1ш1ШШШ1811^^Ш!®11!!1ВИИ11!1И1ШВВ!111Ш1Ш
Аналогичные данные мы получим и при других значениях п.

2. 
   Задачи
   

Задачи представляют собой задание для самостоятельного создания программ на языке Python. Как и в упражнениях, предусматривается написание моду­ля для решения поставленной задачи и программы, которая иллюстрирует работу этого модуля.
Задача 4.1 Напишите программу для вычисления норм вещественных век­торов (||х||_р, р = 1,2,оо) и норм вещественных матриц
п
IHIIi = max £|а_у|,
¹-³-^пы1
/ п \ ^1//2
№=[EKij -
*
п
HU = max EKI-