1.使用444（x）对数据进行插值，并写出误差分析理论。

建立脚本x1=[0.2 0.4 0.6 0.8 1.0];y1=[0.98 0.92 0.81 0.64 0.38];n=length(y1);c=y1(:);for j=2:n %求差商for i=n:-1:jc(i)=(c(i)-c(i-1))/(x1(i)-x1(i-j+1));endendsyms x df d;df(1)=1;d(1)=y1(1);for i=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-filesyms x l;x1=[0.2 0.4 0.6 0.8 1.0];y1=[0.98 0.92 0.81 0.64 0.38];n=length(x1);Ls=sym(0);for i=1:nl=sym(y1(i));for k=1:i-1l=l*(x-x1(k))/(x1(i)-x1(k));endfor k=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.试求3次、合曲线，用图示数据曲线及相应的三种拟合曲线。

建立脚本X=[0.0 0.1 0.2 0.3 0.5 0.8 1.0];Y=[1.0 0.41 0.50 0.61 0.91 2.02 2.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);for k=1:n+1V=1;for j=1:n+1if k~=jV=conv(V,poly(x(j)))/(x(k)-x(j));endendL(k,:)=V;endC=y*Lend%在命令窗口中输入以下的命令：x=[0.0 0.1 0.2 0.3 0.5 0.8 1.0];y=[1.0 0.41 0.50 0.61 0.91 2.02 2.46];cc=polyfit(x,y,4);xx=x(1):0.1:x(length(x));yy=polyval(cc,xx);plot(xx,yy,'r');hold on;plot(x,y,'x');xlabel('x');ylabel('y');x=[0.0 0.1 0.2 0.3 0.5 0.8 1.0];y=[0.94 0.58 0.47 0.52 1.00 2.00 2.46]; %y中的值是根据上面两种拟合曲线而得到的估计数据，不是真实数据function [C,L]=lagran(x,y);xx=0:0.01:1.0;yy=polyval(C,xx);hold on;plot(xx,yy,'b',x,y,'.');命令窗口运行p1 =-6.6221 12.8147 -4.6591 0.9266p2 =2.8853 -12.3348 16.2747 -5.2987 0.9427Y1 =0.9266 0.5822 0.4544 0.5034 0.9730 2.0103 2.4602 Y2 =0.9427 0.5635 0.4399 0.5082 1.0005 1.9860 2.4692 f1 =- (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.2183 110.3671 -57.3264 23.4994 -5.4764 0.94003.对于积分149xdx=-，取不同的步长h。

分别用复合梯形及复合辛普森求积计算积分，给出误差中关于h的函数，并与积分精确值比较两个公式的精度，是否存在一个最小的h，使得精度不能再被改善？syms xx=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.444444由结果（1度要高，且当步长取不同值时即n 越大、h 越小时，积分精度越高。

实验结果说明不存在一个最小的h ,使得精度不能再被改善。

又两个相应的关于h 的误差（余项）其中h 属于a 到b 。

可知h h 愈小，余项愈小，从而积分精度越高。

4.用LU 分解及列主元高斯消去法解线性方程组123410701832.09999962 5.9000015151521021x x x x ⎡⎤⎡⎤⎡⎤-⎢⎥⎢⎥⎢⎥-⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥--⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦输出Ax=b 中系数A=LU 分解的矩阵L 及U ，解向量x 及detA;列主元法的行交换次序，解向量x 及detA;比较两种方法所得的结果。

LU 分解法建立函数function h1=zhijieLU(A,b) %h1各阶主子式的行列式值[n n]=size(A);RA=rank(A);if RA~=ndisp('请注意：因为A 的n 阶行列式h1等于零，所以A 不能进行LU 分解。

A 的秩RA 如下：')RA,h1=det(A);returnendif RA==nfor p=1:nh(p)=det(A(1:p,1:p));endh1=h(1:n);for i=1:nif h(1,i)==0disp('请注意：因为A的r阶主子式等于零，所以A不能进行LU 分解。

A的秩RA和各阶顺序主子式h1依次如下：')h1;RAreturnendendif h(1,i)~=0disp('请注意：因为A的r阶主子式都不等于零，所以A能进行LU分解。

A的秩RA和各阶顺序主子式h1依次如下：')for j=1:nU(1,j)=A(1,j);endfor k=2:nfor i=2:nfor j=2:nL(1,1)=1;L(i,i)=1;if i>jL(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 -7 0 1;-3 2.099999 6 2;5 -1 5 -1;2 1 0 2];b=[8;5.900001;5;1];h1=zhijieLU(A,b)请注意：因为A的r阶主子式都不等于零，所以A能进行LU分解。

A的秩RA 和各阶顺序主子式h1依次如下：RA =4U = 10.0000 -7.0000 0 1.00000 2.1000 6.0000 2.30000 0 -2.1429 -4.23810 -0.0000 0 12.7333L = 1.0000 0 0 0-0.3000 1.0000 0 00.5000 1.1905 1.0000 -0.00000.2000 1.1429 3.2000 1.0000X = -0.2749-1.32981.29691.4398h1 = 10.0000 -0.0000 -150.0001 -762.0001列主元法建立函数function [RA,RB,n,X]=liezhu(A,b)B=[A b];n=length(b);RA=rank(A);RB=rank(B);zhicha=RB-RA;if zhicha>0disp('请注意：因为RA~=RB，所以方程组无解')returnwarning off MATLAB:return_outside_of_loopendif RA==RBif RA==ndisp('请注意：因为RA=RB，所以方程组有唯一解')X=zeros(n,1);C=zeros(1,n+1);for p=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;for k=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);for q=n-1:-1:1X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)))/A(q,q);endelsedisp('请注意：因为RA=RB<n，所以方程组有无穷多解') endendend命令窗口运行A=[10 -7 0 1;-3 2.099999 6 2;5 -1 5 -1;2 1 0 2];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,86109337591031A b ⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥==⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦则解T x (1,1,1,1)=.用MATLAB 内部函数求detA 及A 的所有特征值和cond(A)2。