首先设出要求解的矩阵
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);
end
b(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);
end
y(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
end
end
%这一算法将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
Cholesky分解
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;
tic
v = 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);
end
A(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
实际计算通过计算时间得到同样结果!
于是就结束了!