数值分析计算实习题.doc

插值法1. 下列数据点的插值x 0 1 4 9 16 25 36 49 64y 0 1 2 3 4 5 6 7 8可以得到平方根函数的近似，在区间0,64上作图.(1)用这9个点作8次多项式插值Ls(x).(2)用三次样条(第一边界条件)程序求S(x).从得到结果看在0,64上，哪个插值更精确；在区间0,1上，两种插值哪个更精确？解：(1)拉格朗日插值多项式，求解程序如下syms x l;x1=0 1 4 9 16 25 36 49 64;y1=0 1 2 3 4 5 6 7 8;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).输出结果为Ls = -24221063/63504000*x2+95549/72072*x-1/3048192000*x8-2168879/435456000*x4+19/283046400*x7+657859/10886400*x3+33983/152409600*x5-13003/2395008000*x6(2)三次样条插值，程序如下x1=0 1 4 9 16 25 36 49 64;y1=0 1 2 3 4 5 6 7 8;x2=0:1:64;y3=spline(x1,y1,x2);p=polyfit(x2,y3,3); %得到三次样条拟合函数S=p(1)+p(2)*x+p(3)*x2+p(4)*x3 %得到S(x)输出结果为：S = 2288075067923491/73786976294838206464-2399112304472833/576460752303423488*x+4552380473376713/18014398509481984*x2+999337332656867/1125899906842624*x3(3)在区间0,64上，分别对这两种插值和标准函数作图，plot(x2,sqrt(x2),b,x2,y2,r,x2,y3,y)蓝色曲线为y=x函数曲线,红色曲线为拉格朗日插值函数曲线，黄色曲线为三次样条插值曲线可以看到蓝色曲线与黄色曲线几乎重合，因此在区间0,64上三次样条插值更精确。在0,1区间上由上图看不出差别，不妨代入几组数据进行比较 ，取x4=0:0.2:1x4=0:0.2:1;sqrt(x4) %准确值subs(Ls,x,x4) %拉格朗日插值spline(x1,y1,x4) %三次样条插值运行结果为ans = 0 0.4472 0.6325 0.7746 0.8944 1.0000ans = 0 0.2504 0.4730 0.6706 0.8455 1.0000ans = 0 0.2429 0.4630 0.6617 0.8403 1.0000从这几组数值上可以看出在0,1区间上，拉格朗日插值更精确。数据拟合和最佳平方逼近2. 有实验给出数据表x 0.0 0.1 0.2 0.3 0.5 0.8 1.0y 1.0 0.41 0.50 0.61 0.91 2.02 2.46试求3次、4次多项式的曲线拟合，再根据数据曲线形状，求一个另外函数的拟合曲线，用图示数据曲线及相应的三种拟合曲线。解：(1)三次拟合，程序如下sym x;x0=0.0 0.1 0.2 0.3 0.5 0.8 1.0;y0=1.0 0.41 0.50 0.61 0.91 2.02 2.46;cc=polyfit(x0,y0,3);S3=cc(1)+cc(2)*x+cc(3)*x2+cc(4)*x3 %三次拟合多项式xx=x0(1):0.1:x0(length(x0);yy=polyval(cc,xx);plot(xx,yy,-);hold on;plot(x0,y0,x);xlabel(x);ylabel(y);运行结果S3 =-7455778416425083/1125899906842624+1803512222945437/140737488355328*x-655705280524945/140737488355328*x2+4172976892199509/4503599627370496*x3图像如下(2)4次多项式拟合sym x;x0=0.0 0.1 0.2 0.3 0.5 0.8 1.0;y0=1.0 0.41 0.50 0.61 0.91 2.02 2.46;cc=polyfit(x0,y0,4);S3=cc(1)+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);hold on;plot(x0,y0,x);xlabel(x);ylabel(y);运行结果S3 = 3248542900396215/1125899906842624-3471944732519151/281474976710656*x+4580931990070637/281474976710656*x2-5965836931688425/1125899906842624*x3+8491275464650307/9007199254740992*x4图像如下(3)另一个拟合曲线,新建一个M-file程序如下：function 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)/(x(k)-x(j); end end L(k,:)=V;endC=y*L在命令窗口中输入以下的命令：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中的值是根据上面两种拟合曲线而得到的估计数据，不是真实数据C,L=lagran(x,y);xx=0:0.01:1.0;yy=polyval(C,xx);hold on;plot(xx,yy,b,x,y,.);图像如下解线性方程组的直接解法3. 线性方程组Ax=b的A及b为A=10 7 8 77 5 6 58 6 10 97 5 9 10，b=32233331，则解x=1111.用MATLAB内部函数求detA及A的所有特征值和cond(A)2.若令A+A=10 7 8.1 7.27.08 5.04 6 58 5.98 9.89 96.99 5 9 9.98，求解(A+A)(x+x)=b,输出向量x和|x|2.从理论结果和实际计算两方面分析线性方程组Ax=b解得相对误差|x|2/|x|2及A的相对误差|A|2/|A|2的关系.解：(1)程序如下clear;A=10 7 8 7;7 5 6 5;8 6 10 9;7 5 9 10;det(A)cond(A,2)eig(A)输出结果为ans = 1ans = 2.9841e+003ans = 0.0102 0.8431 3.8581 30.2887(2)程序如下A=10 7 8.1 7.2;7.08 5.04 6 5;8 5.98 9.89 9;6.99 5 9 9.98;b=32 23 33 31;x0=1 1 1 1;x=Ab %扰动后方程组的解x1=x-x0 %x的值norm(x1,2) %x的2-范数运行结果为x = -9.5863 18.3741 -3.2258 3.5240x1 = -10.5863 17.3741 -4.2258 2.5240ans = 20.9322程序如下A=10 7 8.1 7.2;7.08 5.04 6 5;8 5.98 9.89 9;6.99 5 9 9.98;A0=10 7 8 7;7 5 6 5;8 6 10 9;7 5 9 10;b=32 23 33 31;x0=1 1 1 1;x=Ab;x1=x-x0;norm(x1,2);A1=A-A0 ; %A的值norm(x1,2)/norm(x0,2) % |x|2/|x|2的值norm(A1,2)/norm(A0,2) %|A|2/|A|2的值输出结果为ans = 10.4661ans =0.0076可见A相对误差只为0.0076，而得到的结果x的相对误差就达到了10.4661，该方程组是病态的，A的条件数为2984.1远远大于1，当A只有很小的误差就会给结果带来很大的影响。非线性方程数值解法4. 求下列方程的实根(1) x2-3x+2-ex=0;(2) x3+2x2+10x-20=0.要求：(1)设计一种不动点迭代法，要使迭代序列收敛，然后再用斯特芬森加速迭代，计算到|x(k)-x(k-1)|10(-8)为止。(2)用牛顿迭代，同样计算到|x(k)-x(k-1)|10(-8)。输出迭代初值及各次迭代值和迭代次数k，比较方法的优劣。解:(1)先用画图的方法估计根的范围ezplot(x2-3*x+2-exp(x);grid on;可以估计到方程的根在区间（0,1）;选取迭代初值为x0=0.5；构造不动点迭代公式x(k+1)=( x(k)2+2-ex(k)/3;(x)= ( x2+2-ex)/3;当0x1时，0(x)1; (x)1;因此迭代公式收敛。程序如下：format long;f=inline(x2+2-exp(x)/3)disp(x=);x=feval(f,0.5);disp(x);Eps=1E-8;i=1;while 1 x0=x; i=i+1; x=feval(f,x); disp(x); if abs(x-x0)Eps break; endendi,x运行结果为f = Inline function: f(x) = (x2+2-exp(x)/3x= 0.20042624309996 0.27274906509837 0.25360715658413 0.25855037626494 0.25726563633509 0.25759898516219 0.25751245451483 0.25753491361525 0.25752908416796 0.25753059723833 0.25753020451046 0.25753030644564 0.25753027998767 0.25753028685501i = 14x = 0.25753028685501斯特芬森加速法，程序如下:format long;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+2-exp(x)/3+x);disp(x=);x=feval(f,0.5);disp(x);Eps=1E-8;i=1;while 1 x0=x; i=i+1; x=feval(f,x); disp(x); if abs(x-x0)Eps break; endendi,x运行结果为x= 0.25868442756579 0.25753031771981 0.25753028543986 0.25753028543986i = 4x = 0.25753028543986牛顿迭代法，程序如下：format long;x=sym(x);f=sym(x2-3*x+2-exp(x);df=diff(f,x);FX=x-f/df;Fx=inline(FX);disp(x=);x1=0.5;disp(x1);Eps=1E-8;i=0;while 1 x0=x1; i=i+1; x1=feval(Fx,x1); disp(x1); if abs(x1-x0)1E10 break; end if abs(x-x0)Eps break; endendi,x运行结果：f = Inline function: f(x) = (-2*x2-10*x+20)1/3x= 0.16666666666667 6.09259259259259 -38.38843164151806 -8.478196837919431e+002 -4.763660785374071e+005 -1.512815059604763e+011i = 6x = -1.512815059604763e+011迭代6次后x的值大得令人吃惊,表明构造的式子并不收敛.也无法构造出收敛的不动点公式牛顿迭代法，程序如下：format long;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;while 1 x0=x1; i=i+1; x1=feval(Fx,x1); disp(x1); if abs(x1-x0)Eps break; endendi,x1运行结果：x= 1.50000000000000 1.37362637362637 1.36881481962396 1.36880810783441 1.36880810782137i = 4x1 = 1.36880810782137比较三种方法，牛顿法的收敛性比较好，相比不动点迭代法要构造出收敛的公式比较难，牛顿法迭代次数也较少，收敛速度快，只是对初值的要求很高，几种方法各有利弊，具体采用哪种也需因题而异。常微分方程初值问题数值解法5. 给定初值问题y=-50y+50x2+2x,0x1;y(0)=1/3;用经典的四阶R-K方法解该问题，步长分别取h=0.1,0.025,0.01,计算并打印x=0.1i(i=0,1,10)各点的值，与准确值y(x)=1/3e(-50x)+x2比较。解：取步长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;for i=1:n x=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+2*K2+2*K3+K4)/6;end%准确值temp=;f=dsolve(Dy=-50*y+50*x2+2*x,y(0)=1/3,x);df=zeros(1,n+1);for i=1:n+1 temp=subs(f,x,X(i); df(i)=double(vpa(temp);enddisp( 步长 四阶经典R-K法 准确值);disp(X,Y,df);运行结果： 步长 四阶经典R-K法 准确值 1.0e+010 * 0 0.00000000003333 0.00000000003333 0.00000000001000 0.00000000046055 0.00000000000122 0.00000000002000 0.00000000630625 0.00000000000400 0.00000000003000 0.00000008640494 0.00000000000900 0.00000000004000 0.00000118436300 0.00000000001600 0.00000000005000 0.00001623545110 0.00000000002500 0.00000000006000 0.00022256067134 0.00000000003600 0.00000000007000 0.00305093542778 0.00000000004900 0.00000000008000 0.04182323921740 0.00000000006400 0.00000000009000 0.57332690347809 0.00000000008100 0.00000000010000 7.85935630083771 0.00000000010000%画图观察结果figure;plot(X,df,k*,X,Y,-r);grid on;title(四阶经典R-K法解常微分方程);legend(准确值,四阶经典R-K法);当x值接近1的时候，偏离准确值太大。当步长h=0.025时，将上面程序中的h改为0.025即可，运行结果： 步长 四阶经典R-K法 准确值 0 0.33333333333333 0.33333333333333 0.10000000000000 0.10313524034288 0.01224598233303 0.20000000000000 0.04428527327599 0.04001513330992 0.30000000000000 0.05196795755507 0.09000010196744 0.40000000000000 0.09395731149439 0.16000000068705 0.50000000000000 0.16034531435757 0.25000000000463 0.60000000000000 0.248