数值分析算法 MATLAB 实践 线性方程组 分解法
Lu分解法
% LUmethod分解矩阵
function [L,U]=LUmethod(A)
[rows,~]=size(A);
temp_mat=A;
L=zeros(rows);
for i=1:rows
coefficient=temp_mat(:,i);
coefficient=coefficient./coefficient(i);
coefficient(1:i)=0;
L(:,i)=coefficient;
temp_mat=-coefficient*temp_mat(i,:)+temp_mat;
end
U=temp_mat;
L(eye(rows)==1)=1;
end
% LUmethod分解矩阵 求解线性方程组
function x=LUsolve(L,U,b)
[rows,~]=size(L);
aug_mat=[L,b];
for i=1:rows
aug_mat(i,:)=aug_mat(i,:)./aug_mat(i,i);
coefficient=-aug_mat(:,i);
coefficient(1:i)=0;
aug_mat=coefficient*aug_mat(i,:)+aug_mat;
end
aug_mat=[U,aug_mat(:,rows+1:end)];
for i=rows:-1:1
aug_mat(i,:)=aug_mat(i,:)./aug_mat(i,i);
coefficient=-aug_mat(:,i);
coefficient(i:end)=0;
aug_mat=coefficient*aug_mat(i,:)+aug_mat;
end
x=aug_mat(:,rows+1:end);
x=x';
end
function solution=LuFunmethon(M, Presion)
% LU分解 M为用户输入的增广矩阵
% Precision为用户所输入的精度要求
if nargin==2
try
digits(Precision);
cath
disp('你输入的精度有误');digits(10);
end
else
digits(10);
end
A=vpa(M)
row=size(A,1);
col=size(A,2);
if ndims(A)~=2|(col-row)~=1
disp('矩阵的大小有误');
return
end
if det(M(:,1:row))==0
disp('该方程的系数矩阵行列式为零');
return
end
%% 调用系统的LU命令
[L,U,P]=lu(double(A));
%% 回代求解过程
for i=row:-1:1
temp=U(i,col);
for k=i+1:row
temp=vpa(temp-t_solution(k)*U(i,k));
end
t_solution(i)=vpa(temp/U(i,i));
end
for i=1:row
temp=t_solution(i);
for k=1: i-1
temp=vpa(temp-t_solution(k)*U(i,k));
end
solution(i)=temp;
end
end
Cholesky分解
%% 0.平方根法解线性方程组，输出L矩阵和根
%% 1.对称正定矩阵的Cholesky分解
%对称正定矩阵A存在唯一的对角元素均为正数的下三角矩阵L，使得A=L*L'
%这种分解叫做Cholesky分解
A=[3,3,5;3,5,9;5,9,17];
b=[0;-2;-4];
%L=chol(A,'lower')基于矩阵A的对角线和下三角形生成下三角矩阵L，满足方程L*L'=A
L=chol(A,'lower')
%% 2.由Ly=b得到y
y=L\b;
%% 3.由L_转置*x=y得到方程组的解x
x=L'\y%输出线性方程组的根
function x = Cholesky_method(A,b)
%Cholesky平方根法解方程组
%A为方程组的系数矩阵 b为方程组的右端项;
n = length(A);
L = zeros(n);
for k = 1:n
delta = A(k,k);
for j = 1:k-1
delta = delta-L(k,j)^2;
end
L(k,k) = sqrt(delta);
for i = k+1:n
L(i,k) = A(i,k);
for j = 1:n-1
L(i,k) = L(i,k)-L(i,j)*L(k,j);
end
L(i,k) = L(i,k)/L(k,k);
end
end
L
x =zeros(n,1);
y =zeros(n,1);
y(1) = b(1)/L(1,1);
for i = 2:n
ly = 0;
for j = 1:i-1
ly = ly+L(i,j)*y(i);
end
y(i) = (b(i)-ly)/L(i,i);
end
x(n) = y(n)/L(n,n);
for i = n-1:-1:1
lx = 0;
for j = i+1:n
lx = lx+L(j,i)*x(j);
end
x(i) = (y(i)-lx)/L(i,i);
end
end
%Cholesky平方根法解方程组
%A为方程组的系数矩阵 b为方程组的右端项;
A6= [1 2 -1;2 5 1; -1 1 14];
b6 = [3;4;3];
x6 = Cholesky_method(A6,b6);
disp(['方程组的解：x6= ']);
disp(x6)
Cholesky分解--改进平方根法
function x=chol_ldlt_method(A,b)
%function x=chol_ldlt_method(A,b)
%Cholesky改进平方根法解方程组
%A为方程组的系数矩阵 b为方程组的右端项;
n = length(A);
L = eye(n);
D = zeros(n);
d = zeros(1,n);
T = zeros(n);
for k =1:n
d(k) = A(k,k);
for j = 1:k-1
d(k)=d(k)-L(k,j)*T(k,j);
end
for i=k+1:n
T(i,k) = A(i,k);
for j = 1:k-1
T(i,k) =T(i,k) -T(i,j)*L(k,j);
end
L(i,k) = T(i,k)/d(k);
end
end
D = diag(d);
L
D
x =zeros(n,1);
y =zeros(n,1);
d1 = zeros(n,1);
d1 = diag(D);
y(1) = b(1);
for i =2:n
ly = 0;
for k=1:i-1
ly = ly+L(i,k)*y(k);
end
y(i) = b(i)-ly;
end
x(n) = y(n)/d1(n);
for i = n-1:-1:1
lx = 0;
for k=i+1:n
lx = lx+L(k,i)*x(k);
end
x(i) = y(i)/d1(i)-lx;
end
%function x=chol_ldlt_method(A,b)
%Cholesky改进平方根法解方程组
%A为方程组的系数矩阵 b为方程组的右端项;
A7= [1 2 -1;2 5 1; -1 1 14];
b7 = [3;4;3];
x7=chol_ldlt_method(A7,b7);
disp(['方程组的解：x7= ']);
disp(x7)
奇异值分解法
%奇异值分解计算线性方程组
a=[6.5 -1 -1 3.6
6.2 7 -5 4
3 2.1 -6 4.8
1 5.6 3.7 2.1];
b=[12.3 21.4 -7.8 21]';
[u,s,v]=svd(a)
x=v*inv(s)*u'*b
