《数值分析作业(程序)》由会员分享，可在线阅读，更多相关《数值分析作业(程序)（8页珍藏版）》请在金锄头文库上搜索。
1、习题二习题二4 4（1 1） 方法一：Doolittle 分解法function my_LU(a, b) n = length(a); l = zeros(n, n); u = zeros(n, n); for i=1:n l(i,i) = 1; end u(1,1:n) = a(1,1:n); l(2:n, 1) = a(2:n, 1) ./ u(1,1); for r=2:n for i=r:n u(r, i) = a(r, i) - sum(l(r,1:r-1) .* (u(1:r-1,i); end for i=r+1:n if (r=n) l(i, r) = (a(i, r) - s
2、um(l(i,1:r-1) .* (u(1:r-1,r)./u(r,r); end end end L=l,U=u y(1) = b(1); for i=2:n y(i) = b(i) - sum(l(i, 1:i-1).*y(1:i-1); end x(n) = y(n)/u(n,n); for i=n-1:-1:1 x(i) = (y(i) - sum(u(i,i+1:n).*x(i+1:n)./u(i,i); end x=x end在 matlab 命令栏中输入矩阵 a 及 b，再输入函数名，即可运行出结果。方法二： Crout 分解A=4 3 2 1;3 4 3 2;2 3 4 3;1
3、 2 3 4;b=0 1 -1 0; format ratn=length(A); l=zeros(n); u=eye(n); for i=1:nl(i,1)=a(i,1); end for i=2:n u(1,i)=a(1,i)/l(1,1); end for i=2:nfor j=i:nll=0;for k=1:i-1ll=ll+l(j,k)*u(k,i);endl(j,i)=a(j,i)-ll;endfor j=i+1:nuu=0;for k=1:i-1uu=uu+l(i,k)*u(k,j);endu(i,j)=(a(i,j)-uu)/l(i,i);end end l u n=lengt
4、h(A); y(1)=b(1)/l(1,1); for i=2:nsum1=0;for k=1:i-1sum1=sum1+l(i,k)*y(k);end y(i)=(b(i)-sum1)/l(i,i); end x(n)=y(n); for i=n-1:-1:1sum2=0;for k=i+1:nsum2=sum2+u(i,k)*x(k);endx(i)=y(i)-sum2; endx方法三：平方根法A=4 3 2 1;3 4 3 2;2 3 4 3;1 2 3 4;b=0 1 -1 0; n=length(A); for i=1:n for j=1:n if A(i,j)=A(j,i) di
5、sp(输入的方阵不是对称矩阵，请重新输入); return; end end end d=eig(A); for i=1:n if d(i)=0 disp(输入的矩阵不是正定矩阵，请重新输入); return; else break; end end disp(输入的矩阵可以进行 Cholesky 分解); 矩阵 A 具体 Cholesky 分解程序如下 L=zeros(n,n); L(1,1)=sqrt(A(1,1); for i=2:1:n L(i,1)=A(i,1)/L(1,1); end for j=2:1:n sum1=0; for k=1:1:j-1 sum1=sum1+L(j,k
6、)2; end L(j,j)=sqrt(A(j,j)-sum1); for i=j+1:1:n; sum2=0; for k=1:1:j-1 sum2=sum2+L(i,k)*L(j,k); end L(i,j)=(A(i,j)-sum2)/L(j,j); end end L 求解线性方程组的解 L1=L;%L 的转置 L1 n=length(A); y(1)=b(1)/L(1,1); for i=2:nsum1=0;for k=1:i-1sum1=sum1+L(i,k)*y(k);end y(i)=(b(i)-sum1)/L(i,i); end x(n)=y(n)/L(n,n); for i
7、=n-1:-1:1sum2=0;for k=i+1:nsum2=sum2+L1(i,k)*x(k);endx(i)=(y(i)-sum2)/L1(i,i); end x4 4（4 4）方法一：Doolittle 分解A=2 1 1;1 3 2;1 2 2;b=4 6 5; format rat n=length(A); L=eye(n); U=zeros(n); for i=1:nU(1,i)=A(1,i); end for i=2:n L(i,1)=A(i,1)/U(1,1); endfor i=2:nfor j=i:nuu=0;for k=1:i-1uu=uu+L(i,k)*U(k,j);
8、endU(i,j)=A(i,j)-uu;endfor j=i+1:nll=0;for k=1:i-1ll=ll+L(j,k)*U(k,i);endL(j,i)=(A(j,i)-ll)/U(i,i);end end L U求解 x n=length(A); y(1)=b(1); for i=2:nsum1=0;for k=1:i-1sum1=sum1+L(i,k)*y(k);end y(i)=b(i)-sum1; end x(n)=y(n)/U(n,n); for i=n-1:-1:1sum2=0;for k=i+1:nsum2=sum2+U(i,k)*x(k);endx(i)=(y(i)-su
9、m2)/U(i,i); end x方法二： Crout 分解A=2 1 1;1 3 2;1 2 2;b=4 6 5; format rat n=length(A); l=zeros(n);u=eye(n); for i=1:nl(i,1)=a(i,1); endfor i=2:nu(1,i)=a(1,i)/l(1,1); endfor i=2:nfor j=i:nll=0;for k=1:i-1ll=ll+l(j,k)*u(k,i);endl(j,i)=a(j,i)-ll;endfor j=i+1:nuu=0;for k=1:i-1uu=uu+l(i,k)*u(k,j);endu(i,j)=(
10、a(i,j)-uu)/l(i,i);end end l u求解 x n=length(A); y(1)=b(1)/l(1,1); for i=2:nsum1=0;for k=1:i-1sum1=sum1+l(i,k)*y(k);end y(i)=(b(i)-sum1)/l(i,i); end x(n)=y(n); for i=n-1:-1:1sum2=0;for k=i+1:nsum2=sum2+u(i,k)*x(k);endx(i)=y(i)-sum2; end x方法三：平方根法A=2 1 1;1 3 2;1 2 2;b=4 6 5; n=length(A); for i=1:n for
11、j=1:n if A(i,j)=A(j,i) disp(输入的方阵不是对称矩阵，请重新输入); return; end end end d=eig(A); for i=1:n if d(i)=0 disp(输入的矩阵不是正定矩阵，请重新输入); return; else break; end end disp(输入的矩阵可以进行 Cholesky 分解); 矩阵 A 具体 Cholesky 分解程序如下 L=zeros(n,n); L(1,1)=sqrt(A(1,1); for i=2:1:n L(i,1)=A(i,1)/L(1,1); end for j=2:1:n sum1=0; for
12、k=1:1:j-1 sum1=sum1+L(j,k)2; end L(j,j)=sqrt(A(j,j)-sum1); for i=j+1:1:n; sum2=0; for k=1:1:j-1 sum2=sum2+L(i,k)*L(j,k); end L(i,j)=(A(i,j)-sum2)/L(j,j); end end L求解线性方程组的解 L1=L; L1 n=length(A); y(1)=b(1)/L(1,1); for i=2:nsum1=0;for k=1:i-1sum1=sum1+L(i,k)*y(k);end y(i)=(b(i)-sum1)/L(i,i); end x(n)=y(n)/L(n,n); for i=n-1:-1:1sum2=0;for k=i+1:nsum2=sum2+L1(i,k)*x(k);endx(i)=(y(i)-sum2)/L1(i,i); end x1111解：A=2 6;2 6.00001;b1=8 8.00001;b2=8 8.00002; C=cond(A,inf) X1=Ab1 X2=Ab2
