常州大学数值分析课后习题答案第二章第三章第四章节

数值分析作业第二章1、用Gauss消元法求解下列方程组：2x1-x2+3x3=1,(1)4x1+2x2+5x3=4,x1+2x2=7;(2)解：A=[2-131;4254;1207]n=size(A,1);x=zeros(n,1);flag=1;%消元过程fork=1:n-1fori=k+1:nifabs(A(k,k))epsA(i,k+1:n+1)=A(i,k+1:n+1)-A(k,k+1:n+1)*A(i,k)/A(k,k);elseflag=0;returnendendend%回代过程ifabs(A(n,n))epsx(n)=A(n,n+1)/A(n,n);elseflag=0;returnendfori=n-1:-1:1x(i)=(A(i,n+1)-A(i,i+1:n)*x(i+1:n))/A(i,i);endreturnxA=2-13142541207x=9-1-611x1-3x2-2x3=3,(2)-23x1+11x2+1x3=0,x1+2x2+2x3=-1;(2)解：A=[11-3-23;-231110;122-1]n=size(A,1);x=zeros(n,1);flag=1;%消元过程fork=1:n-1fori=k+1:nifabs(A(k,k))epsA(i,k+1:n+1)=A(i,k+1:n+1)-A(k,k+1:n+1)*A(i,k)/A(k,k);elseflag=0;returnendendend%回代过程ifabs(A(n,n))epsx(n)=A(n,n+1)/A(n,n);elseflag=0;returnendfori=n-1:-1:1x(i)=(A(i,n+1)-A(i,i+1:n)*x(i+1:n))/A(i,i);endreturnxA=11-3-23-231110122-1x=0.21240.5492-1.15544、用Cholesky分解法解方程组323x15220x233012x37解：.A=[323;220;3012];b=[537];lambda=eig(A);iflambdaeps&isequal(A,A')[n,n]=size(A);R=chol(A);%解R'y=by(1)=b(1)/R(1,1);ifn1fori=2:ny(i)=(b(i)-R(1:i-1,i)'*y(1:i-1)')/R(i,i);endend%解Rx=yx(n)=y(n)/R(n,n);ifn1fori=n-1:-1:1x(i)=(y(i)-R(i,i+1:n)*x(i+1:n)')/R(i,i);endendx=x';elsex=[];disp('该方法只适用于对称正定的系数矩阵!');endR=1.73211.15471.732100.8165-2.4495001.7321y=2.8868-0.40820.5774x=1.00000.50000.33335.用列主元Doolittle分解法解方程组解：A=[345;-134;-23-5;];345X12b=[2,-26]';-134X2-2[L,U,pv]=luex(A);-23-5X36y=L\b(pv);x=U\y结果如下：x=11-114.已知，计算.解：A=[10099;9998];cond(A,inf)ans=3.9601e+04cond(A,2)ans=3.9206e+0427.编写LU分解法，改进平方根法，追赶法的Matlab程序，并进行相关数值试验。解：LU分解法程序Function[L,U]=lup(A)%lup:LUfactorization%Synopsis:[L,U]=lup(A)%Input:A=coefficientmatrix%Output:L:lowertriangularmatrix%UuppertriangularmatrixFormatshort[m,n]=size(A);Ifm~=n,error(`Amatrixneedstobesquare`);EndPv=(1:n)`;%LUfactorizationFori=1:n-1Pivot=A（i,i）;Fork=i+1;nA(k,i)=A(k,i)/pivot;A(k,i+1;n)=A(k,i+1;n)-A(k,i)*A(i,i+1;n);EndEndL=eye(size(A))+tril(A,-1);%extractLandUU=triu(A)改进平方根法程序Function[x]=ave(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(`wrong`);Break;EndEndEndD(1,1)=A(1,1)；Fori=2;nForj=1;i-1S(i,j)=A(i,j)-sum(S(i,1;i-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(L(i,1;i-1)*L(i,1;i-1)`);EndD(i,i)=A(i,i)-sum(S(i,1;i-1)*D(1;i-1,1;i-1)*y(1;i-1)))/D(i,i);EndY=zeros(n,1);X=zero(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));End追赶法程序Function[x,L,U]=Thomas(a,b,c,f)N=length(b);%对A进行分解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(1,-1);U=diag(u)+diag(c,1);X=zeros(n,1);Y=x;%?求解ly=b?Y(1)=f(1);Fori=2;nY(i)=f(i)-l(i-1)*y(i-1);End%?求解Ux=y?If(u(n)`=0)X(n)=y(n)/u(n);EndFori=n-1;-1;1X(i)=(y(i)-c(i)*x(i+1))/u（i）；End第三章1、设节点x0=0,x1=π/8,x2=π/4,x3=3π/8,x4=π/2,适当选取上述节点用Lagrange插值法分别构造cosx在区间[0,π/2]上的一次，二次和四次插值多项式P1(x)P2(x)和P4(x),并分别计算P1(x),P2(x),P4(x)其中X取π/3。A=fliplr(A);Returnx=[π/8,3π/8];y=cos(x);x0=π/3;[A,Y]=lagrange(x,y,x0);P1=vpa(poly2sym(A),3)YP1=1.19x-0.689Y=0.4729x0=π/3;[A,Y]=lagrange(x,y,x0);P2=vpa(poly2sym(A),3)YP2=x2-0.109x-0.336Y=0.5174x=[0,π/8,π/4,3π/8,π/2];y=cos(x);x0=π/3;[A,Y]=lagrange(x,y,x0);P4=vpa(poly2sym(A),3)YP4=x4+0.00282x3-0.514x2+0.0232x+0.0287Y=0.50017.根据列表函数x0.01.02.03.04.0f(x)0.501.252.753.502.75选取适当的节点，用逐次线性插值法给出三次多项式在2.8处的值。答：Matlab程序function[T,y0]=aitken(x,y,x0,T0)ifnargin==3T0=[];endn0=size(T0,1);m=max(size(x));n=n0+m;T=zeros(n,n+1);T(1:n0,1:n0+1)=T0;T(n0+1:n,1)=x;T(n0+1:n,2)=y;ifn0==0i0=2;elsei0=n0+1;endfori=i0:nforj=3:i+1T(i,j)=fun(T(j-2,1),T(i,1),T(j-2,j-1),T(i,j-1),x0);endy0=T(n,n+1);returnfunction[y]=fun(x1,x2,y1,y2,x)y=y1+(y2-y1)*(x-x1)/(x2-x1);return%选取0、1、3、4四个节点，求三次插值多项式x=[0,1,3,4];y=[0.5,1.25,3.5,2.75];x0=2.8;[T,y0]=aitken(x,y,x0)y0=3.4190000000000008.根据上题中的列表函数，写出差商表，并写出Newton插值多项式N2（x）和N4（x）。答：差商表：xf(x)0.00.51.01.250.752.02.751.1250.3753.03.5010.125-0.254.02.750.5625-0.0625-0.218750.03125由Nn(x)=a0+a1(x-x0)+a2(x-x0)(x-x1)+…+an(x-x0)(x-x1)…(x-xn-1)得N2(x)=0.5+0.75(x-0.00)+0.375(x-0.0)(x-1.0)=0.375x2+0.375x+0.5N4(x)=0.5+0.75(x-0.00)+0.375(x-0.0)(x-1.0)-0.25(x-0.0)(x-1.0)(x-2.0)+0.03125(x-0.0)(x-2.0)(x-1.0)(x-3.0)=0.03125x4-0.4375x3+1.46875x2-0.3125x+0.516、选取适当的函数y=f(x)和插值节点，编写Matlab程序,分别利用Lagrange插值方法,Newton插值方法确定的插值多项式，并将函数y=f(x)的插值多项式和插值余项的图形画在同一坐标系中，观测节点变化对插值余项的影响。答：Matlab程序function[C,D,Y]=newpoly(x0,y0,x)%检验输入参数ifnargin2|nargin3error('IncorrectNumberofInputs');endiflength(x0)~=length(y0)error('Thelengthofx0mustbeequaltoitofy0');endn=length(x0);D=zeros(n,n);D(:,1)=y0';%计算差商表forj=2:nfork=j:nifabs(x0(k)-x0(k-j+1))epserror('DividedbyZero,therearetwonodesarethesame');endD(k,j)=(D(k,j-1)-D(k-1,j-1))/(x0(k)-x0(k-j+1));endend%计算Newton插值多项式的系数C=D(n,n);fork=(n-1):-1:1C=conv(C,poly(x0(k)));m=length(C);C(m)=C(m)+D(k,k);endifnargin==3Y=polyval(C,x);endx=[01234];y=[0.5,1.25,2.75,3.5,2.75];x0=[01234];y0=[0.5,1.25,2.75,3.5,2.75];%用lagrang插值法计算[A,Y]=lagrange(x,y,x0)Lx=vpa(poly2sym(A),4)%用newton插值法计算[C,D,X]=newpoly(x0,y0,x)Nx=vpa(poly2sym(C),4)%绘制两者图像plot(x,Y,'b*',x0,X,'r-')A=0.5000-0.31251.4687-0.43750.0313Y=0.50001.25002.75003.500