MATLAB解决数值分析问题(总9页

1、1.已知函数在下列各点的值为：xi0.20.40.60.81.0f(xi)0.980.920.810.640.38使用4次牛顿插值多项式P4（x）及4次拉格朗日插值多项式L4（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:j c(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 %求牛顿差

2、值多项式 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*x4 + 0.83333*x3 - 1.1042*x2 + 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

3、.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:n l=sym(y1(i); for k=1:i-1 l=l*(x-x1(k)/(x1(i)-x1(k); end for k=i+1:n l=l*(x-x1(k)/(x1(i)-x1(k); end Ls=Ls+l;endLs=simplify(Ls) %为所求插值多项式Ls(x).figureezplot(Ls,0.2,1.08);输出结果为Ls =- (25*x4)/48 + (5*x3)/6 - (53*x2)/48 + (23*x

4、)/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.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

5、,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+1 V=1; for j=1:n+1 if k=j V=conv(V,poly(x(j)

6、/(x(k)-x(j); end end L(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.

7、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.4602Y2 =0.9427 0.5635 0.4399 0.5082 1.0005 1.9860 2

8、.4692f1 =- (3727889208212549*x3)/562949953421312 + (3607024445890881*x2)/281474976710656 - (1311410561049893*x)/281474976710656 + 4172976892199517/4503599627370496f2 =(406067862549531*x4)/140737488355328- (6943889465038337*x3)/562949953421312 + (4580931990070649*x2)/281474976710656 - (29829184658442

9、19*x)/562949953421312 + 8491275464650319/9007199254740992C = 40.6746 -110.2183 110.3671 -57.3264 23.4994 -5.4764 0.94003.对于积分 ，取不同的步长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

10、,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）可知复合辛普森法求积分精度明显比复合梯形法求积的精度要高，且当步长取不同值时即越大、越小时，积分精度越高。实验结果说明不存在一个最小的,使得精度不能再被改善。又两个相应的关于的误差（余项）其中属于到。可知h愈小

11、，余项愈小，从而积分精度越高。4.用LU分解及列主元高斯消去法解线性方程组输出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=n disp(请注意：因为A的n阶行列式h1等于零，所以A不能进行LU分解。A的秩RA如下：) RA,h1=det(A); returnendif RA=n for p=1:n h(p)=det(A(1:p,1:p); end h1=h(1

12、:n); for i=1:n if h(1,i)=0 disp(请注意：因为A的r阶主子式等于零，所以A不能进行LU分解。A的秩RA和各阶顺序主子式h1依次如下：) h1;RA return end end if h(1,i)=0 disp(请注意：因为A的r阶主子式都不等于零，所以A能进行LU分解。A的秩RA和各阶顺序主子式h1依次如下：) for j=1:n U(1,j)=A(1,j); end for k=2:n for i=2:n for j=2:n L(1,1)=1;L(i,i)=1; if ij L(1,1)=1;L(2,1)=A(2,1)/U(1,1);L(i,1)=A(i,1)

13、/U(1,1); L(i,k)=(A(i,k)-L(i,1:k-1)*U(1:k-1,k)/U(k,k); else U(k,j)=A(k,j)-L(k,1:k-1)*U(1:k-1,j); end end end end h1;RA,U,L, X=inv(U)*inv(L)*b 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;h1=zhijieLU(A,b)请注意：因为A的r阶主子式都不等于零，所以A能进行LU分解。A的秩RA和各阶顺序主子式h1依次如下：RA =4U = 10.000

14、0 -7.0000 0 1.0000 0 2.1000 6.0000 2.3000 0 0 -2.1429 -4.2381 0 -0.0000 0 12.7333L = 1.0000 0 0 0 -0.3000 1.0000 0 0 0.5000 1.1905 1.0000 -0.0000 0.2000 1.1429 3.2000 1.0000X = -0.2749 -1.3298 1.2969 1.4398h1 = 10.0000 -0.0000 -150.0001 -762.0001列主元法建立函数function RA,RB,n,X=liezhu(A,b)B=A b;n=length(b

15、);RA=rank(A);RB=rank(B);zhicha=RB-RA;if zhicha0 disp(请注意：因为RA=RB，所以方程组无解) return warning off MATLAB:return_outside_of_loopendif RA=RB if RA=n disp(请注意：因为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:n m=B(k,p)/B(p,p

16、); B(k,p:n+1)=B(k,p:n+1)-m*B(p,p:n+1); end end b=B(1:n,n+1);A=B(1:n,1:n);X(n)=b(n)/A(n,n); for q=n-1:-1:1 X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)/A(q,q); end else disp(请注意：因为RA=RB0 jX=Ajb;Pnb=norm(b,p);pnjx=norm(jX,p);deltaX=jX-X; pnjdX=norm(deltaX,p);jxX=pnjdX/pnjX; pnX=norm(X,p);xX=pnjdX/pnX; pndb=norm(deltab,p);xAb=pndb/pnb;pnbj=norm(jb,p);xAbj=pndb/pnbj; Xgxx=Acp*xAb;endif pndA0 jX=jAb;deltaX=jX-X;pnX=norm(X,p); pnjdX=norm(deltaX,p); pnjX=norm(jX,p);jxX=pnjdX/pnjX;xX=pnjdX/pnX; pnjA=norm(jA,p);pnA=norm(A,p); pndA=norm(deltaA,p);xAbj