Чисельні методи. Лабораторний практикум/Розв'язування систем лінійних алгебраїчних рівнянь
Один з варіантів звіту до лаболаторної. Бажано доповнити іншими методами.
Постановка задачі
ред.Розв'язати СЛАР методом квадратного кореня і Якобі.
Знайти число обумовленості та нев'язку.
Аналіз
ред.Метод квадратного кореня
ред.1) Так як симетрична, то ми можемо розкласти її на добуток матриць , так, що елементи матриці , при . А - взагалі діагональна матриця, тобто ненульові елементи містяться тільки не на діагоналі.
Рокзлад на такі матриці відбувається за такими формулами:
Спочатку , і , .
Потім для для всіх i від 1 до n, ;
;
2) З рівняння знаходимо . Здавалось би, який сенс, робити скільки дій (складність методу, до речі, , щоб в кінці прийти до ще одного рівняння? Може тому, що матриця системи стає трикутна? Але поясніть мені, чим нам наприклад метод Гауса не підходить?
3) З рівняння знаходимо x.
Метод Якобі
ред.Щоб воно розв'язувалось ітераційними методами потрібно щоб виконувались певні умови:
- Діагональна перевага: (кожне число на діагоналі більше по модулю за суму модулів всіх інших чисел в тому ж рядку).
Метод полягає в поступовому уточненні розв'язку.
Спочатку беремо будь-який - вектор початкового наближення.
А потім для кожного обчислюємо: Спочатку для k=0, потім для 1, і так до тих пір, поки наш розв'язок не буде достатньо точним. Він має збігатись до правильного.
Код
ред.# coding=utf-8
from numpy import *
n=input("n=")
A=[]
for i in range(0,n):#Генеруємо матрицю а
row=[]
div=1.0
for j in range(0,n):
if i==j:
row.append(n)
else:
div+=1
row.append(1.0/div)
A.append(row)
A=array(A) #Страшний код, але який є
b=arange(n)+1
print "A=",A
print "b=",b
def sign(x): #Знак числа
if x<0:
return -1.0
if x>0:
return 1.0
return 0.0
def abs(x): # Модуль
return x*sign(x)
from math import sqrt
S=zeros((n,n))
D=zeros((n,n))
S[0,0]=abs(A[0,0])# Далі починаємо метод квадратного кореня
D[0,0]=sign(A[0,0])
for j in range(1,n):
S[0,j]=A[0,j]/(D[0,0]*S[0,0])
for i in range(0,n):
j=i
s=0
for k in range(0,i):
s+=S[k,i]**2*D[k,k]
D[i,i]=sign(A[i,i] - s)
S[i,i]=sqrt(abs(A[i,i] - s))
for j in range(i+1,n):
s=0
for k in range(0,i):
s+=S[k,i]*S[k,j]*D[k,k]
S[i,j]=(A[i,j]-s)/(S[i,i]*D[i,i])
print "S=",S
print "D=",D
def solve(A,b,n):
"""Знаходить x з рівняння Ax=b, де A - нижньотрикутна матриця"""
x=empty(n)
for j in range(n):
sum=0
for i in range(j):
sum+=x[i]*A[j,i]
x[j]=(b[j]-sum)/A[j,j]
return x
y=solve(dot(S.transpose(),D),b,n)
def solve2(A,b,n):
"""Знаходить x з рівняння Ax=b, де A - верхньотрикутна матриця"""
x=empty(n)
j=n
while j>0:
j-=1
sum=0
for i in range(j+1,n):
sum+=x[i]*A[j,i]
x[j]=(b[j]-sum)/A[j,j]
return x
x=solve2(S,y,n)
print "Метод квадратного кореня"
print "x=",x
print "Вектор нев'язок:", dot(A,x)-b
print "Метод Якобі"
def diagonal_prevalence(A,n):
for i in range(n):
s=0.0
for j in range(n):
s+=abs(A[i,j])
if 2*A[i,i]<=s:
return False
return True
if diagonal_prevalence(A,n):
print "Діагональна перевага виконується"
else:
exit(0)
x=zeros(n)
print "Вектор початкового наближення"
print x
def next_vector(x):
y=zeros(n)
for i in range(n):
s=0
for j in range(n):
if j!=i:
s+=A[i,j]/A[i,i]*x[j]
y[i]=-s+b[i]/A[i,i]
return y
for i in range(100):
x=next_vector(x)
print x
print "Вектор нев'язок:",dot(A,x)-b
Вивід програми
ред.n=4 A= [[ 4. 0.5 0.33333333 0.25 ] [ 0.5 4. 0.33333333 0.25 ] [ 0.5 0.33333333 4. 0.25 ] [ 0.5 0.33333333 0.25 4. ]] b= [1 2 3 4] S= [[ 2. 0.25 0.16666667 0.125 ] [ 0. 1.98431348 0.14698618 0.11023964] [ 0. 0. 1.98761598 0.10714492] [ 0. 0. 0. 1.99016135]] D= [[ 1. 0. 0. 0.] [ 0. 1. 0. 0.] [ 0. 0. 1. 0.] [ 0. 0. 0. 1.]] Метод квадратного кореня x= [ 0.09043838 0.37615267 0.65299078 0.93002614] Вектор нев'язок: [ -1.11022302e-16 -2.22044605e-16 1.50730641e-02 5.39556520e-02] Метод Якобі Діагональна перевага виконується Вектор початкового наближення [ 0. 0. 0. 0.] [ 0.09142006 0.37713434 0.64986161 0.91652828] Вектор нев'язок: [ 0. 0. 0. 0.]
Висновок
ред.Метод Якобі, хоча й працює за час залежний від потрібної точності, але може давати точніші результати.