数值分析各种求积公式的应用.doc

数值分析第二次作业 班级：双控 姓名：王雪萍 学号：201102208第4章习题2请用复合辛普森公式计算下列积分：(1), n=8;(2) ,n=8;(3) ，n=8.分析：辛普森（Simpson）公式是牛顿-科特斯公式当n=2时的情形，也称为三点公式。辛普森公式是利用区间二等分的三个点来进行积分插值。其科特斯系数分别为1/6，4/6，1/6.相应的求积公式可表示为： （1），复化的辛普森公式讲区间a ,b分为n等份，在每个子区间xk,xk+1,上采用辛普森公式（1）,将每个小区间2等份记 叫做复合的辛普森公式。Rn（f）为余项。解：在MATLAB中用辛普森求积分的函数paud求积分，以及在编程用复合辛普森公式求解各题如下：（1） 用求定积分的函数int求积分的值：在命令窗口中键入如下的命令：sym xint(x./(4+x.2),0,1)I=vpa(int(x./(4+x.2),0,1)回车后输出积分结果为： I =0.1115717用变步长辛普森求积分函数paud求积分 在m文件中编写如下函数：function y=myfun(x)y=x./(4+x.2);在命令窗口中输入如下指令：S=quad(myfun,0,1)回车后输出积分结果：S= 0.1116 编写求复合辛普森公式的m文件如下： function Simpson1x=0:1/16:1;f= x./(4+x.2);a=0,b=1;%f(1)=1;h=(b-a)/8;s=f(1)+f(17);for k=1:1:8 s=s+4*f(2*k);endfor k=2:1:8 s=s+2*f(2*k-1);ends=(h/6)*s 运行上述程序：s = 0.1116 从上面的结果可以看出，两个结果的误差为0.0000283 （2） 用求定积分的函数int求积分的值：在命令窗口中键入如下的命令：sym xint(sqrt(x),1,9)I=vpa(int (sqrt(x),1,9)回车后输出积分结果为：I =17.3333用变步长辛普森求积分函数paud求积分在m文件中编写如下函数：function y=myfun1(x)y=sqrt(x);在命令窗口中输入如下指令：S=quad(myfun1,1,9)回车后输出积分结果：S= 17.3333 编写求复合辛普森公式的m文件如下： function Simpson2x=1:1/2:9;f= sqrt(x);a=1,b=9;%f(1)=1;h=(b-a)/8;s=f(1)+f(17);for k=1:1:8 s=s+4*f(2*k);endfor k=2:1:8 s=s+2*f(2*k-1);ends=(h/6)*s 运行上述程序：s=17.3332 （3） 用求定积分的函数int求积分的值：在命令窗口中键入如下的命令：sym xint(sqrt(4-sin(x).2),0,pi/6)I=vpa(int (sqrt(4-sin(x).2),0,pi/6)回车后输出积分结果为：I =1.03576用变步长辛普森求积分函数paud求积分在m文件中编写如下函数：function y=myfun2(t)y=sqrt(4-sin(t).2);在命令窗口中输入如下指令：S=quad(myfun2,0,pi/6)回车后输出积分结果：S= 1.0358 编写求复合辛普森公式的m文件如下： function Simpson3x=0:pi/96:pi/6;f=sqrt(4-sin(x).2);a=0,b=pi/6;%f(1)=1;h=(b-a)/8;s=f(1)+f(17);for k=1:1:8 s=s+4*f(2*k);endfor k=2:1:8 s=s+2*f(2*k-1);ends=(h/6)*s运行上述程序：s= 1.0358 误差为0.00004由上述三个小题可以知道，复合辛普森公式的编程只是积分上下限、被积函数以及等分点数的不同，所以我们可以编一个程序，只是改变上面的三个内容，编一个可以普遍使用的程序。程序代码如下：function s=Simpson()%辛普生法求积分clear;clc;options=积分下限a,积分上限b ,插入点相关的值M;topic=seting;lines=1;def=-5,5,1000;h=inputdlg(options,topic,lines,def);a=eval(h1);%积分下限b=eval(h2);%积分上限M=eval(h3);%子区间个数的一半 ，这些变量根据具体的积分而设置%*f=func;%用f来调用被积函数funch=(b-a)/(2*M);s1=0;s2=0;for k=1:M x=a+h*(2*k-1); s1=s1+feval(f,x);endfor k=1:(M-1) x=a+h*2*k; s2=s2+feval(f,x);ends=h*(feval(f,a)+feval(f,b)+4*s1+2*s2)/3;%s是辛普生规则的总计end%定义被积函数funcfunction y=func(x)y=sqrt(x); %不同的积分在程序中只需该这里。% y=x./(4+x.2);%y=sqrt(4-sin(t).2);end 运行程序后，会出现一个设置积分上下限，还有等分点的对话框，在相应的位置填上相应的信息即可。运行结果和上述结果是一样的，因为程序编写方法不同，但完成的任务是一样的。12. 地球卫星轨道是一个椭圆，椭圆的周长的计算公式是 ，这里a是椭圆的半长轴，c是地球中心与轨道中心的的距离，记h为近地点距离，H为远地点距离，R=6371（km）为地球半径，则 a =(2R+H+h)/2,c=(H-h)/2.我国第一课人造地球卫星近地点距离h=439km,远地点距离H=2384km,试求卫星轨道的周长。分析：本题是求解定积分的问题，可以采用复化的辛普森公式求解。根据题中的条件可以计算出c,a的值。解：用matlab求解才c,a如下： a =(2R+H+h)/2 ,c=(H-h)/2. c=972.5, a= 7782.5用求定积分的函数int求积分值：在命令窗口中键入如下的命令：sym xint(4*7782.5*sqrt(1-(972.5/7782.5)2*sin(x).2),0,pi/2)I=vpa(int (4*7782.5*sqrt(1-(972.5/7782.5)2*sin(x).2),0,pi/2)回车后输出积分结果为：I =48707.60111用变步长辛普森求积分函数paud求积分建立M文件：function y=myfun4(x)f=4*7782.5*sqrt(1-(972.5/7782.5)2*sin(x).2);在命令窗口中输入：S=quad(myfun4,0,pi/2)然后回车，可以看到周长为：S= 4.870743851190015e+004km下面编写M文件用复合的辛普森公式计算周长： function Simpson4 e=7782.5 ; %因为下面要用到a，为了避免变量冲突，所以这里将 c=972.5; %a改成e了。x=0:pi/16:pi/2;f=4*e*sqrt(1-(c/e)2*(sin(x).2);a=0,b=pi/2;h=(b-a)/4;s=f(1)+f(9);for k=1:1:4 s=s+4*f(2*k);endfor k=2:1:4 s=s+2*f(2*k-1);ends=(h/6)*s 运行上述程序得出周长为：4.870743851190015e+004km.用梯形求积分函数trapz求积分在命令行输入下列程序：d=pi/16;x=0:d:pi/2;T=d*trapz(4*7782.5*sqrt(1-(972.5/7782.5)2*sin(x).2)回车后，可见 T= 4.870743851190014e+004km用复化的梯形公式求解该题： Matlab程序如下： function t=tixing()%梯形法求积分clear;clc;options=积分下限a,积分上限b ,插入点相关的值n;topic=seting;lines=1;def=0,1,1000;h=inputdlg(options,topic,lines,def);a=eval(h1);%积分下限b=eval(h2);%积分上限n=eval(h3);%子区间个数%*f=func;%用f来调用被积函数funch=(b-a)/(n);t=0;for k=1:1:n-1 x=a+h*(k); t=t + feval(f,x);endformat longt=h*( feval(f,a) +feval(f,b)+2*t)/2;%s是