现在时间是:
当前位置:首 页 >> 矩阵计算>> 教学区>> 文章列表

130213027王者上机作业九QR分解

作者:王者   发布时间:2015-06-23 16:30:21   浏览次数:150
 
例3.3.1
(1) 用Cholesky分解法求解正则化方程组
function x=cholesky1(A,b)
[m,n]=size(A);
C=A'*A;
d=A'*b;
for k=1:n
   C(k,k)=sqrt(C(k,k)-C(k,1:k-1)*C(k,1:k-1)');
   C(k+1:n,k)=(C(k+1:n,k)-C(k+1:n,1:k-1)*C(k,1:k-1)')/C(k,k);
end   
L=tril(C);
U=L';
 
%Ly=b
d(1)=d(1)/L(1,1);
for j=2:n
    d(j)=(d(j)-L(j,1:j-1)*d(1:j-1))/L(j,j);
end
 
%Ux=y
x=zeros(n,1);
x(n)=d(n)/U(n,n);
for j=n-1:-1:1
    x(j)=(d(j)-U(j,j+1:n)*x(j+1:n))/U(j,j);
end
 
运行结果:
clc
U=[1/sqrt(3),1/sqrt(2);-1/sqrt(3),0;1/sqrt(3),-1/sqrt(2)];
V=[sqrt(3)/2,1/2;-1/2,sqrt(3)/2];
b=[3;2;-1];
x1=[1;sqrt(3)];
for k=[10^5,10^7,10^9]
    S=[sqrt(2)*k,0;0,sqrt(2)];
    A=U*S*V';
   x=cholesky1(A,b)
   r=norm(A*x-b)
   
   r2=norm(x-x1)
   r1=norm(A*x1-b)
end
 
 
x =
    1.0000
    1.7321
r =
    2.4495
r2 =
  4.7677e-007
r1 =
    2.4495
 
 
x =
    0.9922
    1.7186
r =
    2.4496
r2 =
    0.0155
r1 =
    2.4495
 
 
x =
   Inf
   Inf
r =
   NaN
r2 =
   Inf
r1 =
2.4495
(2)列主元Gauss消去法求解正则化方程组
function x=gauss2(A,b)
C=A'*A;
b=A'*b;
[m,n]=size(C);
tol=1e-14;
p=[];
if n==m
   for k=1:n-1
       [r,c]=find(C==max(abs(C(k:n,k))));
       p(k)=r(1);
       A([k,p(k)],:)=C([p(k),k],:);     
       if   abs(C(k,k))>tol
            C(k+1:n,k)=C(k+1:n,k)/C(k,k);
            C(k+1:n,k+1:n)=C(k+1:n,k+1:n)-C(k+1:n,k)*C(k,k+1:n);
       
       else
           break;
       end
   end
   L=eye(n)+tril(C,-1);
   U=triu(C);
   % the interchanges of b by P
   for k=1:length(p)
       b([k,p(k)])=b([p(k),k]);
   end
   %Ly=b
   b(1)=b(1)/L(1,1);
   for j=2:n
       b(j)=(b(j)-L(j,1:j-1)*b(1:j-1))/L(j,j);
   end
   %Ux=y
   x=zeros(n,1);
   x(n)=b(n)/U(n,n);
   for j=n-1:-1:1
       x(j)=(b(j)-U(j,j+1:n)*x(j+1:n))/U(j,j);
   end
   
else
    'A is not a square'
end
运行结果:
clc
U=[1/sqrt(3),1/sqrt(2);-1/sqrt(3),0;1/sqrt(3),-1/sqrt(2)];
V=[sqrt(3)/2,1/2;-1/2,sqrt(3)/2];
b=[3;2;-1];
for k=[10^5,10^7,10^9]
    S=[sqrt(2)*k,0;0,sqrt(2)];
    A=U*S*V';
   x=gauss2(A,b)
   r=norm(A*x-b)
   x1=[1;sqrt(3)];
   r2=norm(x-x1)
   r1=norm(A*x1-b)
end
 
x =
    1.0000
    1.7321
r =
    2.4495
r2 =
  2.3849e-007
r1 =
    2.4495
 
 
x =
    0.9981
    1.7287
r =
    2.4495
r2 =
    0.0039
r1 =
    2.4495
 
 
x =
   Inf
   Inf
r =
   NaN
r2 =
   Inf
r1 =
    2.4495
 






上一篇:没有了    下一篇:没有了

Copyright ©2019    计算数学达人 All Right Reserved.

技术支持:自助建站 | 领地网站建设 |短信接口 |燕窝 版权所有 © 2005-2019 lingw.net.粤ICP备16125321号 -5