数值分析计算实习题

1、插值法下列数据点的插值x01491625364964y012345678可以得到平方根函数的近似，在区间0,64上作图.用这9个点作8次多项式插值Ls(x).用三次样条(第一边界条件)程序求S(x).从得到结果看在0,64上，哪个插值更精确；在区间0,1上，两种插值哪个更精确？解：(1)拉格朗日插值多项式，求解程序如下symsxl;x1=01491625364964;y1=012345678;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(

2、k)/(x1(i)-x1(k);endLs=Ls+l;endLs=simplify(Ls)%为所求插值多项式Ls(x).输出结果为Ls=-24221063/63504000*x2+95549/72072*x-1/3048192000*x8-2168879/435456000*x4+19/283046400*x7+657859/10886400*x3+33983/152409600*x5T3003/2395008000*x飞(2)三次样条插值，程序如下x1=01491625364964;y1=012345678;x2=0:1:64;y3=spline(x1,y1,x2);p=polyfit(x2

3、,y3,3);%得到三次样条拟合函数S=p(l)+p(2)*x+p(3)*x2+p(4)*x3%得到S(x)输出结果为：S=2288075067923491/73786976294838206464-2399112304472833/576460752303423488*x+4552380473376713/18014398509481984*x2+999337332656867/1125899906842624*x3在区间0,64上，分别对这两种插值和标准函数作图，plot(x2,sqrt(x2),b,x2,y2,r,x2,y3,y)蓝色曲线为y=Vx函数曲线，红色曲线为拉格朗日插值函数曲线

4、，黄色曲线为三次样条插值曲线100806040200100806040200可以看到蓝色曲线与黄色曲线几乎重合，因此在区间0，64三次样条插值更精确。在0，1区间上由上图看不出差别，不妨代入几组数据进行比较，取x4二0:0.2:1x4=0:0.2:1;sqrt(x4)subs(Ls,x,x4)spline(x1,y1,x4)运行结果为%准确值%拉格朗日插值%三次样条插值ans=00.44720.63250.77460.89441.0000ans=00.25040.47300.67060.84551.0000ans=00.24290.46300.66170.84031.0000从这几组数值上可以

5、看出在0,1区间上，拉格朗日插值更精确。数据拟合和最佳平方逼近有实验给出数据表x0.00.10.20.30.50.81.0y1.00.410.500.610.912.022.46试求3次、4次多项式的曲线拟合，再根据数据曲线形状，求一个另外函数的拟合曲线，用图示数据曲线及相应的三种拟合曲线。解：(1)三次拟合，程序如下symx;x0=0.00.10.20.30.50.81.0;y0=1.00.410.500.610.912.022.46;cc=polyfit(x0,y0,3);S3二cc(l)+cc(2)*x+cc(3)*x2+cc(4)*x3%三次拟合多项式xx=x0(1):0.1:x0(l

6、ength(x0);yy=polyval(cc,xx);plot(xx,yy,-);holdon;plot(x0,y0,x);xlabel(x);ylabel(y);运行结果S3=-7455778416425083/1125899906842624+1803512222945437/140737488355328*x-655705280524945/140737488355328*x2+4172976892199509/4503599627370496*x3图像如下(2)4次多项式拟合symx;x0=0.00.10.20.30.50.81.0;y0=1.00.410.500.610.912.0

7、22.46;cc=polyfit(x0,y0,4);S3二cc(l)+cc(2)*x+cc(3)*x2+cc(4)*x3+cc(5)*x4xx=x0(1):0.1:x0(length(x0);yy=polyval(cc,xx);plot(xx,yy,r);holdon;plot(x0,y0,x);xlabel(x);ylabel(y);运行结果S3=3248542900396215/1125899906842624-3471944732519151/281474976710656*x+4580931990070637/281474976710656*x2-5965836931688425/n2

8、5899906842624*x3+8491275464650307/9007199254740992*x4图像如下2.521.5y2.521.5y10.5x(3)另一个拟合曲线,新建一个M-file程序如下：functionC,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*L在命令窗口中输入以下的命令：x=0.00.10.20.30.50.81.0;y=1.00.410.500.610.912

9、.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中的值是根据上面两种拟合曲线而得到的估计数据，不是真实数据C,L=lagran(x,y);xx=0:0.01:1.0;yy=polyval(C,xx);holdon;plot(xx,yy,b,x,y,.);图像如下x解线性方程组的直接解法线性方

10、程组Ax=b的A及b为1078732A=75651078732A=7565，b=2386109337591031则解1111.用MATLAB内部函数求detA及A的所有特征值和cond(A)2若令1078.17.2A+6A=708q504()65，求解(A+SA)(x+6x)=b，输出向量x和85.989.8996.99599.98|6x|2从理论结果和实际计算两方面分析线性方程组Ax=b解得相对误差l|6x|2/|x|2及a的相对误差|6a|2/|a|2的关系.解：(1)程序如下clear;A=10787;7565;86109;75910;det(A)cond(A，2)eig(A)输出结果为

11、ans=1ans=2.9841e+003ans=0.01020.84313.858130.2887(2)程序如下A=1078.17.2;7.085.0465;85.989.899;6.99599.98;b=32233331;x0=1111;x=Ab%扰动后方程组的解x1=x-x0norm(x1,2)x1=x-x0norm(x1,2)%8x的值%5x的2-范数运行结果为x=-9.586318.3741-3.22583.5240 x1=-10.586317.3741-4.22582.5240ans=20.9322程序如下A=1078.17.2;7.085.0465;85.989.899;6.995

12、99.98;A0=10787;7565;86109;75910;b=32233331;x0=1111;x=Ab;x1=x-x0;norm(x1,2);A1=A-AO;%SA的值norm(xl,2)/norm(x0,2)%|x|2/|x|2的值norm(Al,2)/norm(A0,2)%|A|2/|A|2的值输出结果为ans=10.4661ans=0.0076可见A相对误差只为0.0076，而得到的结果x的相对误差就达到了10.4661,该方程组是病态的，A的条件数为2984.1远远大于1,当A只有很小的误差就会给结果带来很大的影响。非线性方程数值解法求下列方程的实根x八2-3x+2-e八x=0

13、;x八3+2x八2+10 x-20=0.要求：(1)设计一种不动点迭代法，要使迭代序列收敛，然后再用斯特芬森加速迭代，计算到|x(k)-x(k-1)|10,-8)为止。(2)用牛顿迭代，同样计算到|x(k)-x(k-1)|10,-8)。输出迭代初值及各次迭代值和迭代次数k,比较方法的优劣。解:(1)先用画图的方法估计根的范围ezplot(，x“2_3*x+2_exp(x)，);gridon;可以估计到方程的根在区间(0,1)；选取迭代初值为x0=0.5；构造不动点迭代公式x(k+1)=(x(k厂2+2-e八x(k)/3;W(x)二(x八2+2-e八x)/3;当0 x1时，0W(x)1;旷(x)

14、、1;因此迭代公式收敛。0.257534913615250.25753491361525程序如下：formatlong;f二inline(x2+2-exp(x)/3)disp(x=);x=feval(f,0.5);disp(x);Eps=1E-8;i=1;while1x0=x;i=i+1;x=feval(f,x);disp(x);ifabs(x-x0)Epsbreak;endendi,x运行结果为f=Inlinefunction:f(x)=(x2+2-exp(x)/3x=0.200426243099960.272749065098370.253607156584130.258550376264

15、940.257265636335090.257598985162190.257512454514830.257529084167960.257530285439860.257530285439860.257530597238330.257530204510460.257530306445640.257530279987670.25753028685501i=14x=0.25753028685501斯特芬森加速法，程序如下:formatlong;f二inline(x-(x2+2-exp(x)/3-x)2/(x2+2-exp(x)/3)2+2-exp(x2+2-exp(x)/3)/3-2*(x2+

16、2-exp(x)/3+x);disp(x=);x=feval(f,0.5);disp(x);Eps=1E-8;i=1;while1x0=x;i=i+1;x=feval(f,x);disp(x);ifabs(x-x0)Epsbreak;endendi,x运行结果为x=0.258684427565790.257530317719810.25753028543986i=4i=40.25753028543986i=40.25753028543986i=4x=0.25753028543986牛顿迭代法，程序如下：formatlong;x=sym(x);f二sym(，x“2-3*x+2-exp(x);df

17、=diff(f,x);FX=x-f/df;Fx=inline(FX);disp(，x=，);x1=0.5;disp(x1);Eps=1E-8;i=0;while1x0=x1;i=i+1;x1=feval(Fx,x1);disp(x1);ifabs(x1-x0)1E10break;endifabs(x-x0)Epsbreak;endendi,x运行结果：f=Inlinefunction:f(x)=(-2*x2-10*x+20厂1/3x=0.166666666666676.09259259259259-38.38843164151806-8.478196837919431e+002-4.76366

18、0785374071e+005-1.512815059604763e+011i=6x=-1.512815059604763e+011迭代6次后x的值大得令人吃惊,表明构造的式子并不收敛.也无法构造出收敛的不动点公式牛顿迭代法，程序如下：formatlong;x=sym(x);f=sym(x3+2*x2+10*x-20);df=diff(f,x);FX=x-f/df;Fx=inline(FX);disp(x=);x1=0.5;disp(x1);Eps=1E-8;i=0;while1x0=x1;i=i+1;x1=feval(Fx,x1);disp(x1);ifabs(x1-x0)Epsbreak;

19、endendi,x1运行结果：x=1.500000000000001.373626373626371.368814819623961.368808107834411.36880810782137i=4x1=1.36880810782137比较三种方法，牛顿法的收敛性比较好，相比不动点迭代法要构造出收敛的公式比较难，牛顿法迭代次数也较少，收敛速度快，只是对初值的要求很高，几种方法各有利弊，具体采用哪种也需因题而异。常微分方程初值问题数值解法给定初值问题y=-50y+50 x八2+2x,0 x1;y(0)=1/3;用经典的四阶R-K方法解该问题，步长分别取h=0.1,0.025,0.01,计算并打

20、印x=0.1i(i=0,l,，10)各点的值，与准确值y(x)=l/3e八(-50 x)+x八2比较。解：取步长h=0.1,程序如下：%经典的四阶R-K方法clear;F二-50*y+50*x2+2*x;a=0;b=1;h=0.1;n=(b-a)/0.1;X=a:0.1:b;Y=zeros(1,n+1);Y(1)=1/3;fori=1:nx=X(i);y=Y(i);K1=h*eval(F);x=x+h/2;y=y+K1/2;K2=h*eval(F);y=Y(i)+K2/2;K3=h*eval(F);x=X(i)+h;y=Y(i)+K3;K4=h*eval(F);Y(i+1)=Y(i)+(K1+

21、2*K2+2*K3+K4)/6;end%准确值temp=;f二dsolve(Dy=-50*y+50*x2+2*x,y(0)=l/3,x);df=zeros(1,n+1);fori=1:n+1temp=subs(f,x,X(i);df(i)=double(vpa(temp);enddisp（步长四阶经典R-K法准确值）disp（X,Y,df）;运行结果：步长四阶经典R-K法准确值1.0e+010*00.000000000010000.0000000000200000.000000000010000.000000000020000.000000000030000.000000000040000.0

22、00000000050000.000000000060000.000000000070000.000000000080000.000000000090000.00000000010000%画图观察结果figure;plot(X,df,k*,X,Y,gridon;0.000000000033330.000000000460550.000000006306250.000000086404940.000001184363000.000016235451100.000222560671340.003050935427780.041823239217400.573326903478097.8593563

23、00837710.000000000033330.000000000001220.000000000004000.000000000009000.000000000016000.000000000025000.000000000036000.000000000049000.000000000064000.000000000081000.00000000010000titletitle（四阶经典R-K法解常微分方程）;legend（准确值，四阶经典R-K法）;当X值接近1的时候，偏离准确值太大。当步长h=0.025时，将上面程序中的h改为0.025即可，运行结果:步长00.10000000000

24、000步长00.100000000000000.200000000000000.300000000000000.400000000000000.500000000000000.600000000000000.700000000000000.800000000000000.900000000000001.00000000000000图像如下：四阶经典R-K法0.333333333333330.103135240342880.044285273275990.051967957555070.093957311494390.160345314357570.248085701305570.356241884727580.484525906616270.6