常州大学数值分析作业—第二章

20.分别用Jacobi迭代法、Gauss-seidel迭代法解方程组.122,1,122321321321xxxxxxxxx解：Jacobi迭代法收敛133321xxx，Gauss-seidel迭代法不收敛。27.编写LU分解法、改进平方根法、追赶法的Matlab程序，并进行相关数值实验。3.将矩阵1100110211100201A进行Doolittle和Crout分解解：Doolittle分解：结果如下，程序见后面。2.100015001110020112.000010200100001LUACrout分解：结果如下，程序见后面。10002.0100111002012.1100050200100001LUA7.用改进平方根法解方程组000142002511013101144321xxxx解：结果如下，程序见后面。0256.00513.00769.02821.04321xxxx8（2）.用追赶法求解方程组00001210001210001210001210001254321xxxxx解：结果如下，程序见后面。1429.02413.02857.03571.04321xxxxfunction[x,iternum,flag]=jacobi(A,b,x0,delta,max1)%检验输入参数，初始化ifnargin2,error('moreaugmentsareneeded');endifnargin3,x0=zeros(size(b));endifnargin4,delta=1e-13;endifnargin5,max1=100;endifnargin5,error('incorrectnumberofinput');endn=length(b);x=0*b;flag=0;iternum=0;%用Jacobi迭代法解方程组fork=1:max1iternum=iternum+1;fori=1:nifabs(A(i,i))epserror('A(i,i)equaltozero,dividedbyzero');endx(i)=(b(i)-A(i,[1:i-1,i+1:n])*x0([1:i-1,i+1:n]))/A(i,i);x(i)=(b(i)-A(i,[1:i-1,i+1:n])*x0([1:i-1,i+1:n]))/A(i,i);enderr=norm(x-x0);relerr=err/(norm(x)+eps);x0=x;if(errdelta)||(relerrdelta)flag=1;break;endendifflag==1disp('TheJacobimethodconverges.')x=x;elsedisp(['TheJacobimethoddoesnotconvergewith'...,num2str(max1),'iterations'])endreturnA=[12-2;111;221]b=[1;1;1]A=12-2111221b=111jacobi(A,b)TheJacobimethodconverges.ans=-331function[x,iternum,flag]=gseid(A,b,x0,delta,max1)%检验输入参数，初始化ifnargin2,error('moreaugmentsareneeded');endifnargin3,x0=zeros(size(b));endifnargin4,delta=1e-13;endifnargin5,max1=100;endifnargin5,error('incorrectnumberofinput');endn=length(b);flag=0;iternum=0;%用Gauss-Seidel迭代法解方程组fork=1:max1iternum=iternum+1;x=x0;fori=1:nifabs(A(i,i))epserror('A(i,i)equaltozero,dividedbyzero');endx0(i)=(b(i)-A(i,[1:i-1,i+1:n])*x0([1:i-1,i+1:n]))/A(i,i);enderr=norm(x-x0);relerr=err/(norm(x0)+eps);if(errdelta)||(relerrdelta)flag=1;break;endendifflag==1disp('TheGauss-seidelmethodconverges.')x=x0;elsedisp(['TheGauss-seidelmethoddoesnotconvergewith'...,num2str(max1),'iterations'])endreturngseid(A,b)TheGauss-seidelmethoddoesnotconvergewith100iterationsans=1.0e+31*-9.19059.2222-0.0634A=[12-2;111;221]b=[1;1;1]A=12-2111221b=111function[L,U]=Crout(A)%检验输入参数，初始化b=size(A);n=b(1);ifb(1)~=b(2)error('MatrixAshouldbeaSquarematrix.');endifn~=rank(A)error('MatrixAshouldbeFULLRANK.');endL=zeros(n,n);U=zeros(n,n);fori=1:nU(i,i)=1;end%将矩阵A进行Crout分解fork=1:nfori=k:ntemp_sum=0;fort=1:k-1temp_sum=temp_sum+L(i,t)*U(t,k);endL(i,k)=A(i,k)-temp_sum;endforj=k+1:ntemp_sum=0;fort=1:k-1temp_sum=temp_sum+L(k,t)*U(t,j);endU(k,j)=(A(k,j)-temp_sum)/L(k,k);endendreturnA=[1020;0111;20-11;0011]A=1020011120-110011[L,U]=Crout(A)L=1.000000001.0000002.00000-5.00000001.00001.2000U=1.000002.0000001.00001.00001.0000001.0000-0.20000001.0000function[L,U]=Crout(A)%检验输入参数，初始化b=size(A);n=b(1);ifb(1)~=b(2)error('MatrixAshouldbeaSquarematrix.');endifn~=rank(A)error('MatrixAshouldbeFULLRANK.');endL=zeros(n,n);U=zeros(n,n);fori=1:nU(i,i)=1;end%将矩阵A进行Doolittle分解fork=1:nfori=k:ntemp_sum=0;fort=1:k-1temp_sum=temp_sum+L(i,t)*U(t,k);endL(i,k)=A(i,k)-temp_sum;endforj=k+1:ntemp_sum=0;fort=1:k-1temp_sum=temp_sum+L(k,t)*U(t,j);endU(k,j)=(A(k,j)-temp_sum)/L(k,k);endendreturnA=[1020;0111;20-11;0011]A=1020011120-110011[L,U]=Doolittle(A)L=1.000000001.0000002.000001.0000000-0.20001.0000U=1.000002.0000001.00001.00001.000000-5.00001.00000001.2000function[x]=improvecholesky(A,b,n)%检验输入参数，初始化L=zeros(n,n);D=diag(n,0);S=L*D;fori=1:nL(i,i)=1;endfori=1:nforj=1:nif(eig(A)=0)|(A(i,j)~=A(j,i))disp('MatrixAshouldbeaPositive-definitematrix.');break;endendendD(1,1)=A(1,1);%用改进平方根法解方程组fori=2:nforj=1:i-1S(i,j)=A(i,j)-sum(S(i,1:j-1)*L(j,1:j-1)');L(i,1:i-1)=S(i,1:i-1)/D(1:i-1,1:i-1);endD(i,i)=A(i,i)-sum(S(i,1:i-1)*L(i,1:i-1)');endy=zeros(n,1);x=zeros(n,1);fori=1:ny(i)=(b(i)-sum(L(i,1:i-1)*D(1:i-1,1:i-1)*y(1:i-1)))/D(i,i);endfori=n:-1:1x(i)=y(i)-sum(L(i+1:n,i)'*x(i+1:n));endreturnA=[41-10;13-10;-1-152;0024]b=[1;0;0;0]A=41-1013-10-1-1520024b=1000n=4[x]=improvecholesky(A,b,n)n=4x=0.2821-0.07690.0513-0.0256function[L,U,x]=pursue(a,b,c,f)n=length(b);u(1)=b(1);fori=2:nif(u(i-1)~=0)l(i-1)=a(i-1)/u(i-1);u(i)=b(i)-l(i-1)*c(i-1);elsebreak;endendL=eye(n)+diag(l,-1);U=diag(u)+diag(c,1);x=zeros(n,1);y=x;y(1)=f(1);fori=2:ny(i)=f(i)-l(i-1)*y(i-1);endif(u(n)~=0)x(n)=y(n)/u(n);endfori=n-1:-1:1x(i)=(y(i)-c(i)*x(i+1))/u(i);enda=[1-1-1-1];b=[22222];c=[-1-1-1-1];f=[10000];[L,U,x]=pursue(a,b,c,f)L=1.000000000.50001.00000000-0.40001.00000000-0.62501.00000000-0.72731.0000U=2.0000-1.000000002.5000-1.000000001.6000-1.000000001.3750-1.000000001.2727x=0.3571-0.2857-0.2143-0.1429