数值分析计算实习题.docx
《数值分析计算实习题.docx》由会员分享,可在线阅读,更多相关《数值分析计算实习题.docx(19页珍藏版)》请在沃文网上搜索。
1、插值法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
2、 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
3、 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,6
4、4上,分别对这两种插值和标准函数作图,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
5、.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
6、=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/1407374
7、88355328*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);
8、运行结果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
9、 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
10、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
11、 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
12、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);
13、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)设计一种不动点迭代法,要使迭代序列收敛,然后再用斯特芬森加速迭代,计算到
14、|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(
15、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.25752908
16、416796 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=
17、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
18、+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+0
19、11迭代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.
20、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比
- 1.请仔细阅读文档,确保文档完整性,对于不预览、不比对内容而直接下载带来的问题本站不予受理。
- 2.下载的文档,不会出现我们的网址水印。
- 3、该文档所得收入(下载+内容+预览)归上传者、原创作者;如果您是本文档原作者,请点此认领!既往收益都归您。
下载文档到电脑,查找使用更方便
10 积分
下载 | 加入VIP,下载更划算! |
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值 分析 计算 实习