通過前代法後代法列主元Cholesky求解Ax b

2021-10-10 00:20:45 字數 2725 閱讀 8754

矩陣分解

function a = matrix_builder(n)

%對於乙個n,我們按課本120頁要求生成(n-1)^2*(n-1)^2的矩陣

a = diag(repmat([4], 1, (n-1)^2))+diag(repmat([-1], 1, n*(n-2)), 1)+diag(repmat([-1], 1, n*(n-2)), -1)-diag(repmat([1],1,(n-1)*(n-2)),1-n)-diag(repmat([1],1,(n-1)*(n-2)),n-1);

for i = 1:n-2

a((n-1)*i,(n-1)*i+1)=0;

a((n-1)*i+1,(n-1)*i)=0;

end

前代法後代法主要是用來求解a為上三角或者下三角的情況,我們首先通過列主元消去或者是cholesky分解將乙個一般的矩陣轉化為兩個上下三角矩陣的乘積,然後通過上下三角線性方程組的求解演算法,即可解出解

function b = predecessor(l,b)

n = length(b);

for j=1:n

b(j) = b(j)/l(j,j);

b(j+1:n) = b(j+1:n) - b(j)*l(j+1:n,j);

endb(n) = b(n)/l(n,n);

function y = descendant(u,y)

n = length(y);

for j =n:-1:2

y(j) = y(j)/u(j,j);

y(1:j-1) = y(1:j-1)-y(j)*u(1:j-1,j);

endy(1) = y(1)/u(1,1);

我們給出的是通過列主元消去和cholesky分解方法對矩陣進行分解:

function  column_principal_elimination(n)

%首先我們定義乙個指令碼matrix_builder生成(n-1)^2維的方陣

%我們首先採用列主元消去進行求解方程ax=b,計算量為2/3*n^3

a = matrix_builder(n);

b = randn([(n-1)*(n-1),1]);tic

p = eye((n-1)^2);

for i = 1:(n-1)^2

%首先確定最大值的位置,row_max為此時的列最大值,line_indicator為最大值的行指標

a_remaining = a(i:(n-1)^2,i:(n-1)^2);

a_row_i = a(:,i);

row_max = max(abs(a_row_i(i:(n-1)^2)));

line_indicator = find(a_row_i==row_max);

%然後我們設定容器變數temp,並利用line_indicator進行交換行

temp = a(line_indicator,:);

a(line_indicator,:)=a(i,:);

a(i,:)=temp;

temp = p(line_indicator,:);

p(line_indicator,:)=p(i,:);

p(i,:)=temp;

if abs(a(i,i))

a(i:(n-1),i) = a(i+1:n,i)/a(i,i);

a(i+1:n,i+1:n) = a(i+1:n,i+1:n)-a(i+1:n,i)*a(i,i+1:n);

else

break

endend%這一演算法將lu分別儲存在a的上三角和下三角部分,最後p即為置換矩陣.

%計算ly = pb,ux = y;

l = tril(a,-1)+eye((n-1)^2);

u = triu(a);

y = descendant(l,p*b);

x = predecessor(u,y);

toc

function cholesky(n)

%首先我們定義乙個指令碼matrix_builder生成(n-1)^2維的方陣

%我們首先採用cholesky進行求解方程ax=b,計算量為1/3*n^3

a = matrix_builder(n);

b = randn([(n-1)*(n-1),1]);

temp=a\b;

ticv = zeros([1,(n-1)^2]) ;

n=(n-1)^2;

for j = 1:n

for i =1:j-1

v(i) = a(j,i)*a(i,i);

enda(j,j) = a(j,j)-a(j,1:j-1)*v(1:j-1)';

a(j+1:n,j) = (a(j+1:n,j)-a(j+1:n,1:j-1)*v(1:j-1)')/a(j,j);

end%最後我們解ly=b,dl'x=y

l = tril(a,-1)+eye(n);

d = diag(diag(a));

y = predecessor(l,b);

x = descendant(d*l',y);norm(x-temp);

toc

理論分析結果為:

列主元消去計算量為2/3n^3

cholesky分解計算量為1/3n^3

實際計算通過計算時間得到同樣結果!

於是就結束了!

高斯消去法與列主元消去法

兩種消去法的實現主要是,通過函式的實現,傳入引數來實現的。如有其他需要,請另行修改 function time gauss n,a,b b a b tic for k 1 n 1 if a k,k 0 disp the matrix has too many answers,please chang...

列主元的高斯消元法(FORTRAN)

program guass1 real,dimension allocatable arr real,dimension allocatable x real a integer i,j,k,n 輸入 輸出數列 write 請輸入需要計算的係數矩陣的大小n read n allocate arr n...

列主元高斯消去法(C語言)

高斯消元法是將方程組中的一方程的未知數用含有另一未知數的代數式表示,並將其代人到另一方程中,這就消去了一未知數,得到一解 或將方程組中的一方程倍乘某個常數加到另外一方程中去,也可達到消去一未知數的目的。消元法主要用於二元一次方程組的求解。核心 1 兩方程互換,解不變 2 一方程乘以非零數k,解不變 ...