四、线性方程组直接解法解的性态分析

一、我们知道，选主元的LU分解为PA=LU，则解线性方程组Ax=b时，在其两边左乘P，得PAx=Pb，从而有LUx=Pb，最后解得x=U-1L-1Pb。

试利用[L,U,P]=lu(A)设计程序解之。

>> A=[0 3 4;1 -1 1;2 1 2];
>> b=[1;2;3];
>> [L,U,P]=lu(A)
L =
1.0000 0 0
0 1.0000 0
0.5000 -0.5000 1.0000
U =
2 1 2
0 3 4
0 0 2
P =
0 0 1
1 0 0
0 1 0
>> x=inv(U)*inv(L)*P*b
x =
1.1667
-0.3333
0.5000
二、已知A为对称正定矩阵，试编写函数，实现A=RRT的分解，其中R为上三角矩阵。

function R=rfj(A)
[n n]=size(A);
for i=1:n
h(i)=det(A(1:i,1:i));
if h(i)<=0
disp('no')
end
if A-A'==0
disp('yes')
end
end
R(n,n)=sqrt(A(n,n));
for i=n-1:-1:1
R(i,n)=A(i,n)/R(n,n);
end
for j=n-1:-1:1
R(j,j)=sqrt(A(j,j)-(R(j,j+1:n)*R(j,j+1:n)'));
end
for i=j-1:-1:1
R(i,j)=(A(i,j)-sum((R(i,j+1:n))^2))/R(i,j);
end
R=chol(A)
for i=n-1:-1:1
r(i,n)=A(i,n)/r(n,n);
end
for k=j+1:n
for j=n-1:-1:1
r(j,j)=sqrt(A(j,j)-sum((r(j,k))^2));
end
for i=j-1:-1:1
r(i,j)=(A(i,j)-sum((r(i,k)*r(j,k))))/r(j,j);
end
end
三、用矩阵A的P条件数讨论方程组的解的性态函数。

function Acp=zpjwc(A,b,jb,p)
Acp=cond(A,p);dA=det(A);X=A\b;
dertaA=A-jA;
pndA=norm(dertaA,p);dertab=b-jb;
pndb=norm(dertab,p);
if pndb>0
jx=A\jb;pnb=norm(b,p);
pnjx=norm(jx,p);dertax=x-jx;
pnjdx=norm(dertax,p);jxx=pnjdx/pnjx;
pnjx=norm(jx,p);
pnx=norm(x,p);jxx=pnjdx/pnjx;xx=pnjdx/pnx;
pndb=norm(dertab,p);
xAb=pndb/pnb;pnbj=norm(jb,p);xabj=pndb/pnbj;
xgxx=acp*xab;
end
if PndA>0
jX=jA\b;dertaX=X-jX;PnX=norm(X,P);
PnjdX=norm(dertaX,P);
PnjX=norm(jX,P);jxX=PnjdX/PnjX;PnjdX/PnX;
PnjA=norm(jA,P);PnA=norm(A,P);
PndA=norm(dertaA,P);xAbj=PndA/PnjA;
xAb=PndA/PnA;
Xgxx=Acp*xAb;
end
if (Acp>50)&(dA<0.1)
disp('dddd')
Acp,dA,X',jX',xX',jxX',Xgxx',xAb',xAbj'。