MATLAB解决数值分析问题

1.已知函数在下列各点的值为：xi0.20.40.60.81.0f(xi)0.980.920.810.640.38使用4次牛顿插值多项式P4（x）及4次拉格朗日插值多项式L4（x）对数据进行插值，并写出误差分析理论。建立脚本x1=[0.20.40.60.81.0];y1=[0.980.920.810.640.38];n=length(y1);c=y1(:);forj=2:n%求差商fori=n:-1:jc(i)=(c(i)-c(i-1))/(x1(i)-x1(i-j+1));endendsymsxdfd;df(1)=1;d(1)=y1(1);fori=2:n%求牛顿差值多项式df(i)=df(i-1)*(x-x1(i-1));d(i)=c(i)*df(i);enddisp('4次牛顿插值多项式');P4=vpa(collect((sum(d))),5)%P4即为4次牛顿插值多项式,保留小数点后5位数figureezplot(P4,[0.2,1.08]);输出结果为P4=-0.52083*x^4+0.83333*x^3-1.1042*x^2+0.19167*x+0.98插值余项：R4(x)=f(5)(ξ)/(5!)*(x-0.6)*(x-0.4)*(x-0.8)*(x-1)*(x-0.2)新建一个M-filesymsxl;x1=[0.20.40.60.81.0];y1=[0.980.920.810.640.38];n=length(x1);Ls=sym(0);fori=1:nl=sym(y1(i));fork=1:i-1l=l*(x-x1(k))/(x1(i)-x1(k));endfork=i+1:nl=l*(x-x1(k))/(x1(i)-x1(k));endLs=Ls+l;endLs=simplify(Ls)%为所求插值多项式Ls(x).figureezplot(Ls,[0.2,1.08]);输出结果为Ls=-(25*x^4)/48+(5*x^3)/6-(53*x^2)/48+(23*x)/120+49/50插值余项：R4(x)=f(5)(ξ)/(5!)*(x-0.6)*(x-0.4)*(x-0.8)*(x-1)*(x-0.2)2.由实验给出数据表x0.00.10.20.30.50.81.0y1.00.410.500.610.912.022.46试求3次、4次多项式的曲线拟合，再根据数据曲线形状，求一个另外函数的拟合曲线，用图示数据曲线及相应的三种拟合曲线。建立脚本X=[0.00.10.20.30.50.81.0];Y=[1.00.410.500.610.912.022.46];p1=polyfit(X,Y,3)p2=polyfit(X,Y,4)Y1=polyval(p1,X)Y2=polyval(p2,X)plot(X,Y,'b*',X,Y1,'g-.',X,Y2,'r--')f1=poly2sym(p1)f2=poly2sym(p2)plot(X,Y,'b*',X,Y1,'m-.',X,Y2,'c--')legend('数据点','3次多项式拟合','4次多项式拟合')xlabel('X轴'),ylabel('Y轴')另一个拟合曲线,新建一个M-filefunction[C,L]=lagran(x,y)w=length(x);n=w-1;L=zeros(w,w);fork=1:n+1V=1;forj=1:n+1ifk~=jV=conv(V,poly(x(j)))/(x(k)-x(j));endendL(k,:)=V;endC=y*Lend%在命令窗口中输入以下的命令：x=[0.00.10.20.30.50.81.0];y=[1.00.410.500.610.912.022.46];cc=polyfit(x,y,4);xx=x(1):0.1:x(length(x));yy=polyval(cc,xx);plot(xx,yy,'r');holdon;plot(x,y,'x');xlabel('x');ylabel('y');x=[0.00.10.20.30.50.81.0];y=[0.940.580.470.521.002.002.46];%y中的值是根据上面两种拟合曲线而得到的估计数据，不是真实数据function[C,L]=lagran(x,y);xx=0:0.01:1.0;yy=polyval(C,xx);holdon;plot(xx,yy,'b',x,y,'.');命令窗口运行p1=-6.622112.8147-4.65910.9266p2=2.8853-12.334816.2747-5.29870.9427Y1=0.92660.58220.45440.50340.97302.01032.4602Y2=0.94270.56350.43990.50821.00051.98602.4692f1=-(3727889208212549*x^3)/562949953421312+(3607024445890881*x^2)/281474976710656-(1311410561049893*x)/281474976710656+4172976892199517/4503599627370496f2=(406067862549531*x^4)/140737488355328-(6943889465038337*x^3)/562949953421312+(4580931990070649*x^2)/281474976710656-(2982918465844219*x)/562949953421312+8491275464650319/9007199254740992C=40.6746-110.2183110.3671-57.326423.4994-5.47640.94003.对于积分104ln9xxdx，取不同的步长h。分别用复合梯形及复合辛普森求积计算积分，给出误差中关于h的函数，并与积分精确值比较两个公式的精度，是否存在一个最小的h，使得精度不能再被改善？symsxx=0:0.001:1;%h为步长,可分别令h=0.1,0.05,0.001,0.005y=sqrt(x).*log(x+eps);T=trapz(x,y);T=vpa(T,7)f=inline('sqrt(x).*log(x)',x);S=quadl(f,0,1);S=vpa(S,7)运行结果：T=-0.4443875S=-0.444444h复合梯形复合辛普森0.2-0.378909-0.4444440.3-0.3315311-0.4444440.001-0.4443875-0.4444440.1-0.4170628-0.444444由结果（1）、（2）可知复合辛普森法求积分精度明显比复合梯形法求积的精度要高，且当步长取不同值时即n越大、h越小时，积分精度越高。实验结果说明不存在一个最小的h,使得精度不能再被改善。又两个相应的关于h的误差（余项）其中h属于a到b。可知hh愈小，余项愈小，从而积分精度越高。4.用LU分解及列主元高斯消去法解线性方程组123410701832.099999625.9000015151521021xxxx输出Ax=b中系数A=LU分解的矩阵L及U，解向量x及detA;列主元法的行交换次序，解向量x及detA;比较两种方法所得的结果。LU分解法建立函数functionh1=zhijieLU(A,b)%h1各阶主子式的行列式值[nn]=size(A);RA=rank(A);ifRA~=ndisp('请注意：因为A的n阶行列式h1等于零，所以A不能进行LU分解。A的秩RA如下：')RA,h1=det(A);returnendifRA==nforp=1:nh(p)=det(A(1:p,1:p));endh1=h(1:n);fori=1:nifh(1,i)==0disp('请注意：因为A的r阶主子式等于零，所以A不能进行LU分解。A的秩RA和各阶顺序主子式h1依次如下：')h1;RAreturnendendifh(1,i)~=0disp('请注意：因为A的r阶主子式都不等于零，所以A能进行LU分解。A的秩RA和各阶顺序主子式h1依次如下：')forj=1:nU(1,j)=A(1,j);endfork=2:nfori=2:nforj=2:nL(1,1)=1;L(i,i)=1;ifijL(1,1)=1;L(2,1)=A(2,1)/U(1,1);L(i,1)=A(i,1)/U(1,1);L(i,k)=(A(i,k)-L(i,1:k-1)*U(1:k-1,k))/U(k,k);elseU(k,j)=A(k,j)-L(k,1:k-1)*U(1:k-1,j);endendendendh1;RA,U,L,X=inv(U)*inv(L)*bendendend命令窗口运行A=[10-701;-32.09999962;5-15-1;2102];b=[8;5.900001;5;1];h1=zhijieLU(A,b)请注意：因为A的r阶主子式都不等于零，所以A能进行LU分解。A的秩RA和各阶顺序主子式h1依次如下：RA=4U=10.0000-7.000001.000002.10006.00002.300000-2.1429-4.23810-0.0000012.7333L=1.0000000-0.30001.0000000.50001.19051.0000-0.00000.20001.14293.20001.0000X=-0.2749-1.32981.29691.4398h1=10.0000-0.0000-150.0001-762.0001列主元法建立函数function[RA,RB,n,X]=liezhu(A,b)B=[Ab];n=length(b);RA=rank(A);RB=rank(B);zhicha=RB-RA;ifzhicha0disp('请注意：因为RA~=RB，所以方程组无解')returnwarningoffMATLAB:return_outside_of_loopendifRA==RBifRA==ndisp('请注意：因为RA=RB，所以方程组有唯一解')X=zeros(n,1);C=zeros(1,n+1);forp=1:n-1[Y,j]=max(abs(B(p:n,p)));C=B(p,:);B(p,:)=B(j+p-1,:);B(j+p-1,:)=C;fork=p+1:nm=B(k,p)/B(p,p);B(k,p:n+1)=B(k,p:n+1)-m*B(p,p:n+1);endendb=B(1:n,n+1);A=B(1:n,1:n);X(n)=b(n)/A(n,n);forq=n-1:-1:1X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)))/A(q,q);endelsedisp('请注意：因为RA=RBn，所以方程组有无穷多解')endendend命令窗口运行A=[10-701;-32.09999962;5-15-1;2102];b=[8;5.900001;5;1];[RA,RB,n,X]=liezhu(A,b),H=det(A)请注意：因为RA=RB，所以方程组有唯一解RA=4RB=4n=4X=0.0000-1.00001.00001.0000H=-762.00015.线性方程组Ax=b的A及b为1078732756523,86109337591031Ab则解Tx(1,1,1,1).用MATLAB内部函数求detA及A的所有特征值和cond(A)2。若令1078.17.27.085.046585.989.8996.99599.98AA求解(A+A)(xx)b，输出向量x和2x.从理论结果和实际计算两方面分析线性