数值分析第二次作业
班级：双控 姓名：王雪萍 学号：201102208
第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文件如下：
数值分析第二次作业
班级：双控 姓名：王雪萍 学号：201102208
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
数值分析第二次作业
班级：双控 姓名：王雪萍 学号：201102208
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';
数值分析第二次作业
班级：双控 姓名：王雪萍 学号：201102208
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如下：
数值分析第二次作业
班级：双控 姓名：王雪萍 学号：201102208
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;
数值分析第二次作业
班级：双控 姓名：王雪萍 学号：201102208
options={'积分下限a','积分上限b' ,'插入点相关的值n'}; topic='seting'; lines=1;
def={'0','1','1000'};
h=inputdlg(options,topic,lines,def); a=eval(h{1});%积分下限 b=eval(h{2});%积分上限 n=eval(h{3});%子区间个数
%********************************************
f='func';%用f来调用被积函数func h=(b-a)/(n); t=0;
for k=1:1:n-1 x=a+h*(k); t=t + feval(f,x);
end
format long
t=h*( feval(f,a) +feval(f,b)+2*t)/2;%s是梯形规则的总计 end
%定义被积函数func function y=func(x)
y=4*7782.5*sqrt(1-(972.5/7782.5)^2*sin(x).^2); end
运行上面的程序得周长为：4.870743851190014e+004。
由上述计算结果可知编程计算与用积分函数计算的结果是一样的。与精确值之间的误差也在允许的范围内。
终上所述，卫星轨道的周长约为48707km.