Перейти до вмісту

Чисельні методи. Лабораторний практикум/Розв'язування систем лінійних алгебраїчних рівнянь

Матеріал з Вікіпідручника

Один з варіантів звіту до лаболаторної. Бажано доповнити іншими методами.

Постановка задачі

[ред.]

Розв'язати СЛАР методом квадратного кореня і Якобі.

Знайти число обумовленості та нев'язку.

Аналіз

[ред.]

Метод квадратного кореня

[ред.]

1) Так як симетрична, то ми можемо розкласти її на добуток матриць , так, що елементи матриці , при . А - взагалі діагональна матриця, тобто ненульові елементи містяться тільки не на діагоналі.

Рокзлад на такі матриці відбувається за такими формулами:

Спочатку , і , .

Потім для для всіх i від 1 до n, ;

;

2) З рівняння знаходимо . Здавалось би, який сенс, робити скільки дій (складність методу, до речі, , щоб в кінці прийти до ще одного рівняння? Може тому, що матриця системи стає трикутна? Але поясніть мені, чим нам наприклад метод Гауса не підходить?

3) З рівняння знаходимо x.

Метод Якобі

[ред.]

Щоб воно розв'язувалось ітераційними методами потрібно щоб виконувались певні умови:

  1. Діагональна перевага: (кожне число на діагоналі більше по модулю за суму модулів всіх інших чисел в тому ж рядку).

Метод полягає в поступовому уточненні розв'язку.

Спочатку беремо будь-який - вектор початкового наближення.

А потім для кожного обчислюємо: Спочатку для 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.]

Висновок

[ред.]

Метод Якобі, хоча й працює за час залежний від потрібної точності, але може давати точніші результати.