实验3.1
 
1 实验目的
 
1.1 实验3.1:分别用复化梯形公式与复化Simpson公式计算下列积分，并与精确积分值相比较，探讨两类积分公式的精度。
 
2 实验内容
 
编写相应的M文件实现下列问题：
 
分别用复化梯形公式与复化Simpson公式计算下列积分，并与精确积分值相比较，探讨两类积分公式的精度。
 
（1），将区间8等分； 
 
 （2），将区间4等分；
 
 （3），将区间6等分；
 
 
 3 实验知识点
 
3.1 复化梯形公式与复化Simpson公式求积分。
 
3.2 求积函数计算定积分
 
4 算法思想
 
3.1 复化梯形公式
 
 
3.2 复化Simpson公式
 
 
5 实验代码及结果
 
（一）实验3.1
 
5.1 ，将区间8等分；
 
 5.1.1复化梯形求积分
 
代码
 
T_quad.m 
 
function [In,er]=T_quad(a,b,n)
 
h=(b-a)/n;
 
x=a:h:b;
 
y=x./(x+4.^2);
 
c=[1  2*ones(1,n-1),1];
 
In=h/2*sum(c.*y);
 
I=quad(f,a,b,1e-5)
 
er=abs(I-In);
 
In,er
 
 
f.m 
 
function f=f(x)
 
f=inline('x./(x+4.^2)')
 
%f=inline('sqrt(x)');
 
%f=inline('sqrt(4-(sinx).^2)')
 
运行结果
 
 
 
5.1.2复化Simpson公式求积分
 
代码
 
复合辛普森方法函数如下：
 
function s=simpsion(f,a,b,n)
 
%复化辛普森公式求积分.
 
h=(b-a)/n;  
 
x=linspace(a,b,2*n+1);
 
y=feval(f,x);
 
s=(h/6)*(y(1)+2*sum(y(3:2:2*n-1))+4*sum(y(2:2:2*n))+y(2*n+1));
 
end
 
运行结果
 
 
 
 
 
5.2 ，将区间4等分
 
5.2.1复化梯形公式求积分
 
代码
 
T_quad.m 
 
function [In,er]=T_quad(a,b,n)
 
h=(b-a)/n;
 
x=a:h:b;
 
y=sqrt(x);
 
c=[1  2*ones(1,n-1),1];
 
In=h/2*sum(c.*y);
 
I=quad(f,a,b,1e-5)
 
er=abs(I-In);
 
In,er
 
 
f.m 
 
function f=f(x)
 
%f=inline('x./(x+4.^2)')
 
f=inline('sqrt(x)');
 
%f=inline('sqrt(4-(sinx).^2)')
 
运行结果
 
 
 
5.2.2复化Simpson公式求积分
 
代码
 
复合辛普森方法函数如下：
 
function s=simpsion(f,a,b,n)
 
%复化辛普森公式求积分.
 
h=(b-a)/n;  
 
x=linspace(a,b,2*n+1);
 
y=feval(f,x);
 
s=(h/6)*(y(1)+2*sum(y(3:2:2*n-1))+4*sum(y(2:2:2*n))+y(2*n+1));
 
end
 
运行结果
 
 
 
5.3 ，将区间6等分
 
 5.3.1复化梯形公式求积分
 
代码
 
T_quad.m 
 
function [In,er]=T_quad(a,b,n)
 
h=(b-a)/n;
 
x=a:h:b;
 
y=sqrt(4-(sin(x)).^2);
 
c=[1  2*ones(1,n-1),1];
 
In=h/2*sum(c.*y);
 
I=quad(f,a,b,1e-5)
 
er=abs(I-In);
 
In,er
 
 
f.m 
 
function f=f(x)
 
%f=inline('x./(x+4.^2)')
 
%f=inline('sqrt(x)');
 
f=inline('sqrt(4-(sinx).^2)')
 
运行结果
 
 
 
 
5.3.2复化Simpson公式求积分
 
代码
 
function s=simpsion(f,a,b,n)
 
%复化辛普森公式求积分.
 
h=(b-a)/n;  
 
x=linspace(a,b,2*n+1);
 
y=feval(f,x);
 
s=(h/6)*(y(1)+2*sum(y(3:2:2*n-1))+4*sum(y(2:2:2*n))+y(2*n+1));
 
end
 
结果
 
 
 
 
6 实验结果分析
 
 
 由复化梯形公式和复化辛普森公式两种方法求解可以看出，两种方法得到的结果相差不是很大。但在一般情况下，当分开的区间数n相等时，复化辛普森得到的结果比复化梯形公式得到的结果更加准确。若想通过复化梯形公式求解得到复化辛普森求解的精确值，就需要选取更大的n，即划分成更多的区间进行求解。
 
 
 
实验3.2
 
1 实验目的
 
已知地球卫星飞行轨迹、部分距离及轨迹周长计算公式等信息，选用适当的求积函数计算定积分，求解卫星轨迹长度。
 
2 实验内容
 
地球卫星飞行轨道是一个椭圆，椭圆周长的计算公式是:
 
 
 式中，a是椭圆的半长轴，c是地球中心与轨道中心（椭圆中心）的距离。令h为近地点距离，H为远地点距离，R=6371（km)为地球半径，则
 
 
 
我国第一颗人造地球卫星近地点距离h=439(km),远地点距离H=3484km,试求卫星轨道长度。
 
3 实验知识点
 
在科学研究和工程技术中，经常遇到积分的计算，虽然有些函数的不定积分可以求出其初等函数表示式，但有更多的函数，它们的不定积分不是初等函数，这样就无法利用牛顿莱布尼兹公式求出其定积分，甚至经常遇到只知道函数在一些离散点的值，但函数表达式未知的情况，在上述情况下就必须以数值方法求定积分的近似值。用数值方法求定积分的近似值，通常称为数值积分。
 
4 算法思想
 
龙贝格求积公式也称为逐次分半加速法。它是在梯形公式、辛普森公式和柯特斯公式之间的关系的基础上，构造出一种加速计算积分的方法，在不增加计算量的前提下提高了误差的精度。在等距基点的情况下，用计算机计算积分值通常都采用把区间逐次分半的方法进行。前一次分割得到的函数值在分半以后仍可被利用,且易于编程。
 
 
5 实验代码
 
代码：
 
R=6371;
 
h=439;
 
H=3484;
 
a=(2*R+H+h)/2;
 
c=(H-h)/2;
 
syms theta
 
f=sqrt(1-(c/a)^2*sin(theta)^2);
 
s=a*int(f,theta,0,pi/2);
 
s=double(vpa(s))
 
 
6 实验结果
 

 
分别用复化梯形公式、复化Simpson 公式计算定积分dx e x ⎰+201，复化梯形公式求定积分:，function I=tquad(x，%复化梯形求积公式，%x 为向量，m=length(y)，%积分自变量的结点数应与它的函数值的个数相同 h=(x(n)-x(1))/(n-1)，a=[1 2*ones(1，I=h/2*sum(a，复化Simpson 公式求定积分：，function I=squad(x，%复化Simpson 求积公式，%x 为向量，m=length(y)，%积分自变量的结点数应与它的函数值的个数相同 if rem(n-1，I=tquad(x，return，end，N=(n-1)/2，h=(x(n)-x(1))/N，a=zeros(1，for k=1:N，a(2*k-1)=a(2*k-1)+1，a(2*k)=a(2*k)+4，a(2*k+1)=a(2*k+1)+1，end，I=h/6*sum(a。
 
取n=2，y)，其中，被积函数自变量的等距结点，n-2) 1]，*y)，y)，其中，被积函数自变量的等距结点，2)~=0，y)，n)，*y)。
 
4， %y 为向量， %y 为向量。
 
8，16分别验证结果(精确值I=4，被积函数在结点处的函数值， n=length(x)，被积函数在结点处的函数值， n=length(x)。
 
导读：文档《数值分析复化梯形公式,复化Simpson公式MATLAB程序》共1页，当前为第1页，大小为20kb，是专业资料、自然科学、数学相关类别的资料，并提供若干种数值分析复化梯形公式,复化Simpson公式MATLAB程序的文本文档下载，如word文档下载、wps文档下载等。以下便是第1页的正文：
 
分别用复化梯形公式、复化Simpson 公式计算定积分dx e x ⎰+201，取n=2,4,8,16分别验证结果(精确值I=4.006994)。
 
复化梯形公式求定积分:
 
function I=tquad(x,y)
 
%复化梯形求积公式，其中，
 
%x 为向量，被积函数自变量的等距结点； %y 为向量，被积函数在结点处的函数值； n=length(x);
 
m=length(y);
 
%积分自变量的结点数应与它的函数值的个数相同 h=(x(n)-x(1))/(n-1);
 
a=[1 2*ones(1,n-2) 1];
 
I=h/2*sum(a.*y);
 
复化Simpson 公式求定积分：
 
function I=squad(x,y)
 
%复化Simpson 求积公式，其中，
 
%x 为向量，被积函数自变量的等距结点； %y 为向量，被积函数在结点处的函数值； n=length(x);
 
m=length(y);
 
%积分自变量的结点数应与它的函数值的个数相同 if rem(n-1,2)~=0
 
I=tquad(x,y);
 
return;
 
end
 
N=(n-1)/2;
 
h=(x(n)-x(1))/N;
 
a=zeros(1,n);
 
for k=1:N
 
a(2*k-1)=a(2*k-1)+1;
 
a(2*k)=a(2*k)+4;
 
a(2*k+1)=a(2*k+1)+1;
 
end
 
I=h/6*sum(a.*y);

 
分别用复化梯形公式、复化Simpson 公式计算定积分dx e x ?+201，取
 
n=2,4,8,16分别验证结果(精确值I=)。
 
复化梯形公式求定积分:
 
function I=tquad(x,y)
 
%复化梯形求积公式，其中，
 
%x 为向量，被积函数自变量的等距结点； %y 为向量，被积函数在结点处的函数值； n=length(x);
 
m=length(y);
 
%积分自变量的结点数应与它的函数值的个数相同 h=(x(n)-x(1))/(n-1);
 
a=[1 2*ones(1,n-2) 1];
 
I=h/2*sum(a.*y);
 
复化Simpson 公式求定积分：
 
function I=squad(x,y)
 
%复化Simpson 求积公式，其中，
 
%x 为向量，被积函数自变量的等距结点； %y 为向量，被积函数在结点处的函数值； n=length(x);
 
m=length(y);
 
%积分自变量的结点数应与它的函数值的个数相同 if rem(n-1,2)~=0
 
I=tquad(x,y);
 
return;
 
end
 
N=(n-1)/2;
 
h=(x(n)-x(1))/N;
 
a=zeros(1,n);
 
for k=1:N
 
a(2*k-1)=a(2*k-1)+1;
 
a(2*k)=a(2*k)+4;
 
a(2*k+1)=a(2*k+1)+1;
 
end
 
I=h/6*sum(a.*y);

从Newton-Cotes的截断误差公式可以看出，当积分区间$[a,b]$较大时，低阶的Newton-Cotes求积公式截断误差都比较大。由于高阶Newton-Cotes求积公式是数值不稳定的，因此通过不断增加阶数来提高求积公式的精度是不可行的。但是，如果将积分区间$[a,b]$分成几个小区间（任意的），在每个小区间上应用Newton-Cotes求积公式，其截断误差必然会减小，然后再把每个小区间上的积分值累加起来，这样却能大大提高整个积分的精度。这种方法称为复化求积方法。
 
常用的复化求积方法采用等分区间的做法，具体如下：
 
将区间$[a,b]$划分为n等分，步长为$H=\frac{(b-a)}{n}$，分点为$x_k=a+kH,k=0,1,2,\cdots ,n$。先用低阶Newton-Cotes求积公式求得每个子区间$[x_k,x_{k+1}]$上的积分值$I_k$，然后将它们累加起来求和，用${\sum }_{k=0}^{n-1}{I}_k$作为所求积分$I={\int }_a^{b}f(x)dx$的近似值。
 
1. 复化梯形公式
 
在区间$[a,b]$上采用复化求积方法，具体使用梯形求积公式进行计算，就得到复化梯形求积公式。用$T_k$表示$f(x)$在子区间$[x_k,x_{k+1}]$上的积分值，$T_n$表示$f(x)$在区间$[a,b]$上的积分值，有：
 $$T_k=\frac{1}{2}H[f(x_k)+f(x_{k+1})]$$
 其中
 $$H=(b-a)/n,x_k=a+kH(k=0,1,2,\cdots ,n)$$
 故
 $$T_n=\sum _{k=0}^{n-1}{T}_k=\frac{1}{2}H\sum _{k=0}^{n-1}[f(x_k)+f(x_{k+1})]$$
 即：
 $$T_n=\frac{1}{2}H[f(x_{k=0})+f(x_{k=n})+2\sum _{k=1}^{n-1}f(x_k)]$$
 或
 $$T_n=\frac{b-a}{2n}[f(a)+f(b)+2\sum _{k=1}^{n-1}f(x_k)]$$
 截断误差用$R_T$表示，由于$f(x)$在区间$[a,b]$上有连续的二阶导数，故有：
 $$R_T=\sum _{k=0}^{n-1}-\frac{1}{12}{H}^3{f}^{(2)}({\eta }_k)=-\frac{1}{12}{H}^3\sum _{k=0}^{n-1}{f}^{(2)}({\eta }_k){\eta }_k\in [x_k,x_{k+1}]$$
 即
 $$R_T=-\frac{1}{12}{H}^3\cdot n\cdot f^{(2)}(\eta )=-\frac{(b-a)}{12}{H}^2{f}^n(\eta )\eta \in [a,b]$$
 
2. 复化Simpson公式
 
在区间$[a,b]$上采用复化求积方法，具体使用Simpson求积公式进行计算，就得到复化Simpson公式。用$S_k$表示$f(x)$在子区间$[x_k,x_{k+1}]$上的积分值，$S_n$表示$f(x)$在区间$[a,b]$上的积分值，有：
 $$S_k=\frac{1}{6}H[f(x_k)+4f(x_{k+\frac{1}{2}})+f(x_{k+1})]$$
 式中，$x_{k+\frac{1}{2}}$为子区间$[x_k,x_{k+\frac{1}{2}}]$的中点，$H=(b-a)/n$。
 $$S_n=\sum _{k=0}^{n-1}=\frac{1}{6}H\sum _{k=0}^{n-1}[f(x_k)+4f(x_{k+\frac{1}{2}})+f(x_{k+1})]$$
 故
 $$S_n=\sum _{k=0}^{n-1}=\frac{1}{6}H[f(a)+4\sum _{k=0}^{n-1}f(x_{k+\frac{1}{2}})+2\sum _{k=1}^{n-1}f(x_k)+f(b)]$$
 截断误差用$R_s$表示，由于$f(x)$在区间$[a,b]$上有连续的四阶导数，故有：
 $$R_S=\sum _{k=0}^{n-1}-\frac{1}{2880}{H}^5{f}^{(4)}({\eta }_k)=-\frac{1}{2880}{H}^5\sum _{k=0}^{n-1}{f}^{(4)}({\eta }_k){\eta }_k\in [x_k,x_{k+1}]$$
 
3. 复化Cotes公式
 
在区间$[a,b]$上采用复化求积方法，具体使用Cotes求积公式进行计算，就得到复化Cotes公式为：
 $$C_n=\frac{1}{90}H[7f(a)+32\sum _{k=0}^{n-1}f(x_{k+\frac{1}{4}})+12\sum _{k=0}^{n-1}f(x_{k+\frac{1}{2}})+32\sum _{k=0}^{n-1}f(x_{k+\frac{3}{4}})+14\sum _{k=1}^{n-1}f(x_k)+7f(b)]$$
 截断误差为：
 $$R_c=-\frac{2(b-a)}{945}(\frac{H}{4}{)}^6{f}^{(6)}(\eta )\eta \in [a,b]$$
 从复化求积的余项公式中可以看出，复化梯形公式、复化Simpson公式、复化Cotes公式的余项和步长的关系为$R_T=O(h^2),R_s=O(h^4),R_c=O(h^6)$。因此，当$H\to 0$或$n\to \infty$时，$T_n,S_n,C_n\to I$。