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

本文格式为Word版，下载可任意编辑数值分析各种求积公式的应用 数值分析其次次作业 班级：双控 姓名：王雪萍 学号：202202208 第4章习题 2．请用复合辛普森公式计算以下积分： 1x(1)?dx, n=8; 04?x2(2) (3) ?91?/6xdx,n=8; 4?sin2?d?，n=8. ?0分析：辛普森（Simpson）公式是牛顿-科特斯公式当n=2时的情形，也称为三点公式辛普森公式是利用区间二等分的三个点来举行积分插值其科特斯系数分别为1/6，4/6，1/6.相应的求积公式可表示为： s?b?a?a?b?，复化的辛普森公式讲区间[a ,b]f(a)?4f()?f(b)? （1）?6?2?分为n等份，在每个子区间[xk,xk+1],上采用辛普森公式（1）,将每个小区间2等份 I??f(x)dx???ak?0bn?ixk?1xkhn?1f(x)dx??[f(xk)?4f(xk?1/2)?f(xk?1)]?Rn(f)6k?0记 hn?1Sn??[f(xk)?4f(xk?1/2)?f(xk?1)]6k?0?h[f(a)?4?f(xk?1/2)?2?f(xk)?f(b)]6k?0k?1n?1n?1 叫做复合的辛普森公式。

Rn（f）为余项 解：在MATLAB中用辛普森求积分的函数paud求积分，以及在编程用复合辛普森公式求解各题如下： （1） 用求定积分的函数int求积分的值： 在命令窗口中键入如下的命令： sym x int(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文件如下： 数值分析其次次作业 班级：双控 姓名：王雪萍 学号：202202208 function Simpson1 x=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); end for k=2:1:8 s=s+2*f(2*k-1); end s=(h/6)*s 运行上述程序：s = 0.1116 从上面的结果可以看出，两个结果的误差为0.0000283 （2） 用求定积分的函数int求积分的值： 在命令窗口中键入如下的命令： sym x int(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 Simpson2 x=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); end for k=2:1:8 数值分析其次次作业 班级：双控 姓名：王雪萍 学号：202202208 s=s+2*f(2*k-1); end s=(h/6)*s 运行上述程序：s=17.3332 （3） 用求定积分的函数int求积分的值： 在命令窗口中键入如下的命令： sym x int(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 Simpson3 x=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); end for k=2:1:8 s=s+2*f(2*k-1); end s=(h/6)*s 运行上述程序：s= 1.0358 误差为0.00004 由上述三个小题可以知道，复合辛普森公式的编程只是积分上下限、被积函数以及等分点数的不同，所以我们可以编一个程序，只是变更上面的三个内容，编一个可以普遍使用的程序。

程序代码如下： function s=Simpson()%辛普生法求积分 clear; clc; options={'积分下限a','积分上限b' ,'插入点相关的值M'}; topic='seting'; 数值分析其次次作业 班级：双控 姓名：王雪萍 学号：202202208 lines=1; def={'-5','5','1000'}; h=inputdlg(options,topic,lines,def); a=eval(h{1});%积分下限 b=eval(h{2});%积分上限 M=eval(h{3});%子区间个数的一半 ，这些变量根据概括的积分而设置 %******************************************** f='func';%用f来调用被积函数func h=(b-a)/(2*M); s1=0; s2=0; for k=1:M x=a+h*(2*k-1); s1=s1+feval(f,x); end for k=1:(M-1) x=a+h*2*k; s2=s2+feval(f,x); end s=h*(feval(f,a)+feval(f,b)+4*s1+2*s2)/3;%s是辛普生规矩的总计 end %定义被积函数func function y=func(x) y=sqrt(x); %不同的积分在程序中只需该这里。

% y=x./(4+x.^2); %y=sqrt(4-sin(t).^2); end 运行程序后，会展现一个设置积分上下限，还有等分点的对话框，在相应的位置填上相应的信息即可运行结果和上述结果是一样的，由于程序编写方法不同，但完成的任务是一样的 12. 地球卫星轨道是一个椭圆，椭圆的周长的计算公式是 s?4a??/2021?(c/a)2sin?d?， 这里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如下： 数值分析其次次作业 班级：双控 姓名：王雪萍 学号：202202208 a =(2R+H+h)/2 ,c=(H-h)/2. c=972.5, a= 7782.5 用求定积分的函数int求积分值： 在命令窗口中键入如下的命令： sym x int(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); end for k=2:1:4 s=s+2*f(2*k-1); end s=(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; 数值分析其次次作业 班级：双控 姓名：王雪萍 学号：202202208 。