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

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

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

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

$A={\begin{pmatrix}n&{\frac {1}{2}}&{\frac {1}{3}}&\ldots &{\frac {1}{n-1}}\\{\frac {1}{2}}&n&{\frac {1}{3}}&\ldots &{\frac {1}{n-1}}\\{\frac {1}{2}}&{\frac {1}{3}}&n&\ldots &{\frac {1}{n-1}}\\\vdots &\vdots &\vdots &\ddots &\vdots \\{\frac {1}{2}}&{\frac {1}{3}}&{\frac {1}{4}}&\ldots &n\\\end{pmatrix}}\ b={\begin{pmatrix}1\\2\\3\\\vdots \\n\end{pmatrix}}$

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

Аналіз

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

1) Так як $A$ симетрична, то ми можемо розкласти її на добуток матриць $A=S^{T}DS$, так, що елементи матриці $s_{ij}=0$, при $i>j$. А $D$ - взагалі діагональна матриця, тобто ненульові елементи містяться тільки не на діагоналі.

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

Спочатку $d_{11}={\mbox{sign }}a_{11}$, і $s_{11}={\sqrt {|a_{11}|}}$, $s_{1j}={\frac {a_{1j}}{d_{11}s_{11}}}\;j={\overline {2..n}}$.

Потім для для всіх i від 1 до n, $d_{ii}={\mbox{sign}}(a_{ii}-\sum _{k=1}^{i-1}s_{ki}^{2}\cdot d_{kk})$;

$s_{ii}={\sqrt {a_{ii}-\sum _{k=1}^{i-1}d_{kk}s_{ki}^{2}}}$;

$s_{ij}={\frac {a_{ij}-\sum _{k=1}^{i-1}s_{ki}s_{kj}d_{kk}}{s_{ii}d_{ii}}}\forall j<i$

2) З рівняння $S^{T}Dy=b$ знаходимо $y$. Здавалось би, який сенс, робити скільки дій (складність методу, до речі, $O(n^{3})$, щоб в кінці прийти до ще одного рівняння? Може тому, що матриця системи стає трикутна? Але поясніть мені, чим нам наприклад метод Гауса не підходить?

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

Метод Якобі

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

1. ${\mbox{det}}A\neq 0$
2. $a_{ii}\neq 0$
3. Діагональна перевага: $|a_{ii}|>\sum _{j\neq i}|a_{ij}|$ (кожне число на діагоналі більше по модулю за суму модулів всіх інших чисел в тому ж рядку).

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

Спочатку беремо будь-який ${\vec {X}}^{(0)}=(X_{1}^{(0)},X_{2}^{(0)},\ldots ,X_{n}^{(0)})$ - вектор початкового наближення.

А потім для кожного $i={\overline {1..n}}$ обчислюємо: $X_{i}^{k+1}=-\sum _{j\neq i}{\frac {a_{ij}}{a_{ii}}}X_{j}^{(k)}+{\frac {b_{i}}{a_{ii}}}$ Спочатку для 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.]

Висновок

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