第5章定积分





第4章介绍了不定积分的概念及求解方式，本章将介绍定积分的定义、定积分的近似计算求解、反常积分的有关计算等内容。本章内容侧重理论与计算，第6章将会对定积分在几何与物理上的应用做更详尽的介绍说明。
5.1本章目标
本章将尝试用MATLAB编写有关计算函数来加深对定义的理解，并解决如下问题： 
(1) 运用定积分的几何意义编写求解函数； 
(2) 定积分的近似计算与数值积分； 
(3) 定积分的符号求解； 
(4) 积分上限函数及其求导； 
(5) 两种反常积分的符号求解与近似求解——Γ函数与Β函数； 
(6) 定积分方法的综合运用； 
(7) 积分的数值求解。
5.2相关命令
下面介绍涉及定积分的MATLAB命令。
(1) int： 定积分求解函数。用法如下： 
 int(expr,var,a,b)： 计算出区间［a,b］上关于var的表达式expr的定积分，如果未指定变量，int将使用symvar确定的默认变量。如果expr是常量，则默认变量为x。
 int(,Name,Value)： 使用一个或多个Name,Value对参数指定选项。例如，'IgnoreAnalyticConstraints',true指定int对积分器应用额外的简化。

(2) trapz： 梯形法求定积分。用法如下： 

 trapz(Y)： 通过梯形法计算Y的近似积分(采用单位间距)，Y的大小确定求积分所沿用的维度。如果Y为向量，则trapz(Y)是Y的近似积分； 如果Y为矩阵，则trapz(Y)对每列求积分并返回积分值的行向量； 如果Y为多维数组，则trapz(Y)对其大小不等于1的第一个维度求积分。
 trapz(X,Y)： 根据X指定的坐标或标量间距对Y进行积分。
 trapz(,dim)： 使用以前的任何语法沿维度dim求积分，必须指定Y，也可以指定X。如果指定X，则它可以是长度等 size(Y,dim)的标量或向量。例如，如果Y为矩阵，则trapz(X,Y,2)对Y的每行求积分。
(3) quad： 抛物线(simpson)法求定积分。用法如下： 
 quad(fun,a,b,tol)： 使用递归自适应simpson积分法求取函数fun从a到b的近似积分，误差为tol，默认为1e-6。
(4) integral： 计算数值积分。用法如下： 
 integral(fun,xmin,xmax)： 使用全局自适应积分和默认误差容限在xmin至xmax间以数值形式为函数 fun 求积分。
 integral(fun,xmin,xmax,Name,Value)： 指定具有一个或多个Name,Value对组参数的其他选项。例如，指定 'WayPoints'，后跟实数或复数向量，为要使用的积分器指示特定点。
(5) gamma： 求gamma积分。用法如下： 
 gamma(x)： 其中x为实数参数。
(6) beta： 求beta积分。用法如下： 
 beta(p,q)： 求beta积分，其中p,q为实数参数。
(7) feval： 将变量数值代入符号函数。用法如下： 
 feval(fun,x1,…,xm)： 将x1,…,xm分别代入fun方程求解。
(8) vpasolve： 求方程数值解。用法如下： 
 vpasolve(fun,var)： 用数值方法求解变量为var的方程fun的根。
5.3定积分的几何意义与近似计算
先回顾一下定积分的定义。
定义51设函数fx在a，b上有界，在a，b中插入若干个分点
a=x0<x1<…<xn－1<xn=b
把区间a，b分成n个小区间
x0，x1,x1，x2,…,xn－1，xn
各个小区间长度依次为
Δx1=x1－x0，Δx2=x2－x1,…,Δxn=xn－xn－1
在每个小区间xi-1，xi上任取一点ξixi－1≤ξi≤xi，作函数值fξi与小区间长度Δxi的乘积fξiΔxii=1，2，…,n，并求和
S=∑ni=1fξiΔxi
记λ=maxΔx1，Δx2,…,Δxn，如果当λ→0时，和S的极限总存在，且与闭区间a，b的分法及点ξi的取法无关，那么称这个极限I为函数fx在区间a，b上的定积分(简称积分)，记作∫bafxdx，即
∫bafxdx=I= limλ→0∑ni=1fξiΔxi
其中，fx叫作被积函数，fxdx叫作被积表达式，x叫作积分变量，a叫作积分下限，b叫作积分上限，a，b叫作积分区间。
下面再讨论定积分的几何意义。在a，b上fx≥0时，∫bafxdx表示由曲线y=fx、两条直线x=a，x=b与x轴所围成的曲边梯形的面积； 在a，b上fx≤0时，表示由曲线y=fx、两条直线x=a，x=b与x轴所围成的曲边梯形位于x轴的下方，∫bafxdx表示上述曲边梯形面积的负值； 在a，b上fx既取得正值又取得负值时，即函数fx的图形某些部分在x轴的上方，而其他部分在x轴下方，此时定积分∫bafxdx表示x轴上方图形面积减去x轴下方图形面积所得之差(如图51所示)。


图51定积分几何意义示意图


例51利用定积分的几何意义计算∫30x2－3x+1dx。
解： 先利用零点存在定理找到所有零点存在的区间，所用函数为root_int，函数代码如下： 



function r=root_int(fun,a,b,h)

% ROOT_INT 寻找函数的所有的零点所在长度为h的区间

% 输入参数:

%——fun:被搜索函数

%——a:搜索区间下界

%——b:搜索区间上界

%——h:区间长度，缺省为(b-a)/100

% 输出参数:

%——r:所有的零点所在的长度为h的区间

% 调用说明:

%   r=root_int(fun,a,b,h):求函数fun在［a,b］区间内所有的零点所在的长度为h的区间



if nargin==3

h=(b-a)/100;

end

a1=a;b1=a1+h;

r=［］;

fun=matlabFunction(fun);

while b1<b

if feval(fun,a1)*feval(fun,b1)<0

% 应用零点存在定理，搜索零点存在区间

r=［r;［a1,b1］］;

a1=b1;b1=a1+h;

else

a1=b1;b1=a1+h;

continue

end

end




然后利用零点存在的小区间以及内置函数fzero找出函数fun在被积区间内所有零点，并找出x轴上方部分及x轴下方部分，分别求面积作差，所用函数为defInt，实现代码如下： 



function I=defInt(fun,a,b,n)

% DEFINNT 应用定积分的几何性质求解定积分

% 输入参数:

% ——fun:被积函数

%——a:搜索区间下界

%——b:搜索区间上界

%——n:区间等分数

% 输出参数:

%——I:应用定积分几何性质求解得出结果

% 调用说明:

%I=defint(fun,a,b,n):求函数fun在［a,b］区间内的几何性质求解







fun=matlabFunction(fun);

r=root_int(fun,a,b);

if ~isempty(r)

N=size(r,1);

solution=zeros(1,N+2);

solution(1)=a;

solution(end)=b;

for i=1:N

solution(i+1)=fzero(fun,r(1,i));

end

Interval=zeros(N+1,0);

for j=1:N+1

Interval(j,1)=solution(j);

Interval(j,2)=solution(j+1);

end




这样，利用上述两个函数，可以求得： 



syms x

fun=x^2-3*x+1;

defInt(fun,0,3,100)




ans = -1.874728644728536
从上述方法可以看到，对于任一确定的正整数n，积分和为
∑ni=1fξiΔxi

都是定积分∫bafxdx的近似值，当n取不同值时，可得到定积分∫bafxdx精度不同的近似解。定义中，将窄条矩形的面积作为窄条曲边梯形面积的近似值。整体上用台阶形的面积作为曲边梯形面积的近似值。事实上，常用的求定积分近似值的方法都采用把区间a，b等分的分法，即用分点a=x0，x1,…,xn－1，xn=b将a，b分成n个长度相等的小区间，每个小区间的长度为Δx=b－an，记fxi=yii=0，1，2，…,n，然后再利用几种可求面积的方法进行求解。
1. 矩形法
矩形左上角与曲线刚好相交时(如图52所示)，有
∫bafxdx≈b－any0+y1+…+yn－1(取ξi=xi－1)



图52矩形法求定积分示意图(矩形左上角与曲线刚好相交)



矩形右上角与曲线刚好相交时(如图53所示)，有
∫bafxdx≈b－any1+y2+…+yn(取ξi=xi)



图53矩形法求定积分示意图(矩形右上角与曲线刚好相交)


2. 梯形法
根据梯形面积公式，有
∫bafxdx≈b－any0+yn2+y1+y2+…+yn-1
梯形法求定积分示意图如图54所示。


图54梯形法求定积分示意图


3. 抛物线法(simpson法)
取n为偶数，有
∫bafxdx≈b－a3ny0+yn+4y1+y3+…+yn-1+2y2+y4+…+yn-2
抛物线法求定积分示意图如图55所示。


图55抛物线法求定积分示意图


数值积分的必要性是由于计算函数的原函数比较困难。原函数可由初等函数表示的为数不多，大部分可积函数的积分无法用初等函数表示，有些甚至无法用解析表达式表示。同时，在实际应用中，只能知道某些积分函数在特定点的取值，此时无法用求原函数方法计算定积分。数值积分为我们提供了不依赖原函数求得函数积分的方法，在数学应用中十分重要，MATLAB提供了几种计算数值积分的方法： trapz函数基于梯形法设计； quad函数基于simpson法编写设计； integral法则可以自动选择最佳的算法进行数值积分。
下面通过几个例子介绍相关函数的用法及具体应用过程。
例52应用矩形法、梯形法及simpson法分别求定积分∫1041+x2dx的近似值。
解： 对于矩形法，以左取法(矩形左上角与曲线刚好相交，见图52)为例，可以编写如下函数： 



function Q=rectangular_left(fun,a,b,n)

% RECTANGULAR_LEFT 根据矩形法求定积分近似值

% 输入参数:

%——fun:函数的MATLAB表述;

%——a,b:积分的下限与上限;

%——n:区间等分数;

% 输出参数:

%——Q:矩形法定积分近似值

% 调用说明:

%Q=rectangular_left(fun,a,b,n):根据矩形法求定积分近似值，等分数为 n


if nargin<4

n=100;

end

X=linspace(a,b,n+1)

x=a

y=eval(fun);

s=0;

for i=1:n

s=s+(b-a)/n*y;

x=X(i);

y=eval(fun);

end

Q=s;




利用该函数以及quad、trapz，取n=100，有： 



syms x

f=4/(1+x^2);

x=［0:1/100:1］;

y=4*(1+x.^2).^(-1);

format long

%用矩形法求出积分近似值，等分长度为1/100

Q_rec=rectangular_left(f,0,1,100)




Q_rec = 3.171374986973629




%用抛物线法求出积分近似值

Q_quad=quad(matlabFunction(f),0,1)




Q_quad = 3.141592653589793



%用梯形法求出积分近似值

Q_trapz=trapz(x,y)




Q_trapz = 3.141575986923129
与本例函数积分的精确值pi相比较可知，精确度排序为： simpson法>梯形法>矩形法。
计算定积分近似值的方法还有很多，如： NewtonCotes公式法、Gauss法、复化梯形法、复化simpson法等，它们精确度更高，依赖的函数更为复杂，被更多地应用于追求高精度数值解的情形，这部分内容将在本章末尾做简要说明。
5.4定积分的符号计算
通常用求原函数的方式计算定积分的值，这依赖于一个重要公式，即： 
定理51(牛顿莱布尼茨公式) 如果函数Fx是连续函数fx在区间a，b上的一个原函数，那么
∫bafxdx=Fb－Fa
在MATLAB中，与不定积分的求解相同，仍然使用int函数求解符号函数的定积分，调用方式为： 

Q=int(fun,x,a,b):fun为待积分符号函数表达式(可为单个方程或矩阵)，x为符号变量，当fun为单一变量时，可以默认。a、b分别为定积分的下限与上限	。

例53计算下列表达式： 
(1) ∫3-1dx1+x2； 
(2) ∫-1-2dxx； 
(3) ∫π0xcosxx2+1

xcoskx3+sinx2+cosxdx。

解： 



syms x k

f1=1/(1+x^2);

f2=1/x;

f3=［x,cos(x)/((x^2+1)^0.5); x*cos(k*x),(3+sin(x))/(2+cos(x))］;

int(f1,3^0.5,-1)



ans = -7π12




int(f2,-2,-1)




ans = -log(2)



int(f3,0,pi)




ans = 
π22∫π0cos(x)x2+1dx


-2sinπk22-πksin(πk)k2log(3)+π3

题(3)中出现了不可积的定积分。对于该类情况，可以使用vpa函数，得到其任意精度的数值解： 



F=int(cos(x)/(x^2 + 1)^(1/2), x, 0, pi);

vpa(F,5)




ans = 0.48827
5.5积分上限函数及其性质
Φx=∫xaftdt(a≤x≤b)
定理52如果函数fx在积分区间a，b上连续，那么积分上限函数
Φx=∫xaftdt
在a，b上可导，并且它的导数
Φ′x=ddx∫xaftdt=fx
定理53如果函数fx在积分区间a，b上连续，那么函数
Φx=∫xaftdt
就是fx在a，b上的一个原函数。
例54计算ddx∫cosxsinxcosπt2dt的导数。
解： 



syms t x

F1=cos(pi*t^2);

F11=int(F1,sin(x),cos(x))




F11 = 2(C(2cos(x))-C(2sin(x)）)2



diff(F11)




ans = -2(2cos(πcos(x)2)sin(x)+2cos(πsin(x)2)cos(x))2
5.6无穷限区间的反常积分
定义52设函数fx在区间a，+∞上连续，任取t>a，作定积分∫tafxdx，再求极限： 
limt→+∞∫tafxdx(51)
若极限存在，则上式称为函数fx在无穷区间a，+∞上的反常积分，记为∫+∞afxdx，即
∫+∞afxdx=limt→+∞∫tafxdx
此时也称反常积分收敛； 若上述极限不存在，则∫+∞afxdx无意义，称为反常积分∫+∞afxdx发散。
类似地，设函数fx在区间-∞，b上连续，任取t<b，作定积分∫btfxdx，再求极限： 
limt→-∞∫btfxdx
若极限存在，则上式称为函数fx在无穷区间-∞，b上的反常积分，记为∫b-∞fxdx，即
∫b-∞fxdx=limt→-∞∫btfxdx
此时也称反常积分收敛； 若上述极限不存在，则∫b-∞fxdx无意义，称为反常积分∫b-∞fxdx发散。
上述反常积分统称为无穷限区间的反常积分。由定义及牛顿莱布尼茨公式可得： 
设Fx为fx在a，+∞上一原函数，若limx→∞Fx存在，则反常积分
∫+∞afxdx=limx→∞Fx-Fa
若limx→∞Fx不存在，则反常积分发散。类似地，在区间-∞，b同理。
由上述可知，对于收敛的无穷限反常积分，依然可以使用int函数求其值。
例55计算下列无穷限反常积分值： 
(1) ∫+∞-∞dx1+x2(2) ∫+∞0te－ptdt(p>0)
解： 



syms x t

f1=1/(1+x^2);

int(f1,-inf,inf)




ans = π



syms a p positive

% 输入变量及变量范围

f2=t*exp(-p*t);

int(f2,t,0,inf)




ans = 1p2
5.7无界函数的反常积分
定义53若函数fx在点a的任意邻域无界，则称点a为fx的瑕点，该类积分又称为瑕积分。

设函数fx在区间a，b上连续，点a为fx的瑕点，任取t>a，作定积分∫btfxdx，再求极限： 
limt→a∫btfxdx(53)
若极限存在，则上式称为函数fx在区间a，b上的反常积分，记为∫bafxdx，即
∫bafxdx=limt→a∫btfxdx
此时也称反常积分收敛； 若上述极限不存在，则∫bafxdx无意义，称为反常积分∫bafxdx发散。
当b为瑕点时亦然，这里不再过多陈述。
例56证明下列无界反常积分，当0<q<1时收敛，当p≥1时发散。
∫badxx－aq
证明： 



syms x a b

int(1/(x-a)^(p+1),x,a,b)




ans = ∞



syms p positive

int(1/(x-a)^(1-p),x,a,b)




ans = （b-a）pp
5.8反常积分的近似计算
5.8.1无界函数的反常积分

无界函数的反常积分的近似计算比较简单，只需去掉相应的瑕点即可用定积分近似计算时的方法计算。
例57计算无界函数∫10xdx1－x2。

解： 



syms x

f1=x/sqrt(1-x^2);

%用抛物线法近似计算反常积分

integral(matlabFunction(f1),0,1)




ans = 1.000000000000077

5.8.2无穷限区间的反常积分
对于无穷限区间的反常积分，由于区间长度无穷，不能直接采用等间分割的方法。下面介绍无穷区间逼近的方法，以方程

∫+∞0fxdx的求解为例。
∫+∞0fxdx=limbn→+∞∫bn0fxdx

取单调递增数列bn，且bn→+∞(n→+∞)，则有
∫+∞0fxdx=∫b00fxdx+∫b1b0fxdx+…+∫bnbn－1fxdx+…
将无穷限函数分为无数个小积分区间，使上面的每个小区间的积分都是正常积分，可用正常方法计算，当∫bnbn－1fxdx≤ε时，停止迭代。
根据上述步骤，利用定积分中quad函数，编写如下求解程序quadInf，解决上述问题。



function Q=quadInf(fun,a,b,tol,range)

% QUANINF 无穷限函数的反常积分的近似计算

% 输入参数:

%——fun:被积函数

%——a,b:积分下限与上限,可为inf

%——tol:容差值，默认为1e-6

%——range:精度要求，迭代终止准则，默认为1e-5

% 输出参数:

%——Q:所求近似积分值

% 调用说明:

%Q=quadInf(fun,a,b):求函数fun在［a,b］上的积分近似计算，a与b可为无穷，下同

%Q=quadInf(fun,a,b,tol):求函数fun在［a,b］上的积分近似计算，容差为tol

%Q=quadInf(fun,a,b,tol,range):求函数fun在［a,b］上的积分近似计算，容差为tol，精度

% 为range



if nargin<5

range=1e-5;

end

if nargin<4






tol=1e-6;

end

if isinf(a) && isinf(b)

Q=quadInf(fun,-inf,0)+quadInf(fun,0,inf);

% 表示若积分范围两端都为无穷，则分两部分递归调用该函数。



elseif isinf(b)

Q=0;I=1;

while I>range

b=a+1;

I=quad(fun,a,b,tol);

Q=Q+I;

a=b;

end

elseif isinf(a)

Q=0;I=1;

while I>range

a=b-1;

I=quad(fun,a,b,tol);

Q=Q+I;

b=a;

end

else

Q=quad(fun,a,b,tol);

end



5.9Γ函数与B函数
下面介绍在理论和实际应用中都非常重要的反常积分Γ函数，并对与它联系紧密同样重要的反常积分Β函数做课本的补充说明。
定义54Γ函数的定义为
Γs=∫+∞0xs－1e－xdx
为查看其定义域，将Γ函数写为
Γs=∫10xs－1e－xdx+∫+∞1xs－1e－xdx
可见，该函数既为无穷限区间的反常积分，当s－1<0时，又为无界函数的反常积分。
s≤0时，等式右边第一个反常积分发散； s>0时，等式右边两个反常积分都收敛，因此该函数定义域为0，+∞。
定义55B函数的定义为
Bp，q=∫10xp－11－xq－1dx
为查看其定义域，将B函数写为
Bp，q=∫120xp－11－xq－1dx+∫112xp－11－xq－1dx
p>0时等式右边第一个反常积分收敛，q>0时等式右边第二个反常积分收敛，因此该函数的定义域为p，q∈0，+∞×0，+∞。
MATLAB提供了求Γ函数与B函数函数值的内置函数gamma，beta，其调用形式为

Y1=gamma(x); Y2=beta(p,q)

参数说明： x,p,q为参数，它必须为实数，Y1，Y2分别为对应积分值。
利用MATLAB绘制Γ函数与简单的B函数的图形(如图56、图57所示)： 



fplot(@gamma)

%运用内置fplot函数直接画出gamma函数的一种情况

legend('Gamma(x)')

%用内置legend函数做图例标注

grid on

xlabel('x'); ylabel('y');

p=linspace(0,5);q=linspace(0,5);

Y2=beta(p,q);

%画出beta函数的一种情况

plot(Y2)

legend('Beta(p,q)')

grid on

xlabel('x'); ylabel('y');






图56Γ函数图形




图57B函数图形


由于函数gamma的特殊性，也有自己的性质。对于整数n： 
(1) gamma(n+1)=factorial(n)=prod(1： n)。
(2) gamma函数的域通过解析延拓延伸到负实数，在负整数处有简单的极点。这种扩展源于以下递归关系的重复应用： 
Γs+1=sΓs
Γ函数与B函数的关系： 
Bp，q=ΓpΓqΓp+q(p>0，q>0)
5.10拓展实例
下面将综合运用以上方法，在具体实例中体会MATLAB在处理定积分问题上提供的便利条件。
例58求函数Fa=∫a－asinaxsinxadx(a>0)的最大值。
解： 



syms a x;

assume(a >=0);%设置参数范围(a>=0)

%a=1与a~=1时分别有两种情况

F = int(sin(a*x)*sin(x/a),x,-a,a)




F = 1-sin(2)2a=1

2a(sin(a2)cos(1)-a2cos(a2)sin(1))a4-1a≠1



F1 = piecewise(a == 1, 1 - sin(2)/2, a ~= 1, (2*a*(sin(a^2)*cos(1) - a^2*cos(a^2)*sin(1)))/(a^4 - 1))




F1 =4912087702119901
9007199254740992a=1

2a1216652631687587sin(a2)2251799813685248-3789648413623927a2cos(a2)450359927370496a4-1
a≠1



%为了后续计算便利，先只考虑a~=1的情况

assumeAlso(a ~= 1);

F= int(sin(a*x)*sin(x/a),x,-a,a)




F = 2a(sin(a2)cos(1)-a2cos(a2)sin(1))a4-1



fplot(F,［0,10］, 'linewidth',3)

xlabel('x'); ylabel('y');

hold on

Fa = diff(F,a)




Fa = 2σ1a4-1+2a(2acos(a2)cos(1)-2acos(a2)sin(1)+2a3sin(a2)sin(1))a4-1-8σ4σ1(a4-1)2
where
σ1=sin(a2)cos(1)-a2cos(a2)sin(1)
本节程序还可以得到F的图形(如图58所示)，根据F的图形，可以很清晰地看到F的最值分布情况及收敛情况，该函数的最大值点为区间［1,2］的极大值点。为了便于进一步确定，在该图上还可以绘制F一阶导数Fa的图形： 



fplot(Fa,［0,10］, '--','linewidth',2)

legend('F','Fa')

grid on

hold off




可以看出，F的极值点恰为Fa的零点，这样在区间［1,2］内求Fa的根，即为区间［1,2］F函数取极大值时a的取值： 



a_max = vpasolve(Fa,a,［1,2］)




a_max = 1.5782881585233198075558845180583



%将a的取值代入F函数，可得最大值点

F_max = vpa(subs(F,a,a_max))




F_max = 1.2099496860938456039155811226054



vpa(int(sin(x)*sin(x),x,-1,1))




ans = 0.54535128658715915230199006704413


图58函数F与Fa的图形


最后对a=1的情况进行验证，结果表明F(1)不是最大值，从而确定在区间［1,2］内的最大值就是函数在整个定义域内的最大值。
5.11定积分的数值求解的补充说明
除了在正文部分提到的定积分的近似计算的方法，还有很多方法对定积分进行数值求解。下面介绍依赖Lagrange插值多项式的插值法求解定积分的数值解。
定义56设函数fx在a，b上的m+1个互异点x0，x1,…,xm上的函数值和若干阶导数值f(j)xi(i=0，1，2,…,m； j=0，1,…,ni－1)是已知的，这里∑mi=0ni=n+1
若存在一个n次多项式pnx，满足如下插值条件
p(j)nxn=f(j)xi(i=0，1，2,…,m； j=0，1，2,…,ni－1)
则称pnx是fx在a，b上关于插值节点(一般简称节点)x0，x1,…,xm的n次插值多项式，而
rnx=fx－pnx
称为插值余项。
根据以上定义，介绍Lagrange插值法。
n0=n1=…=nm=1，m=n
这时n+1个插值条件均为函数值而不包括导数值，即pnx满足
pnxi=fxi，i=0，1，2,…,n
如果能找到一组n次多项式qkx，k=0，1，2,…,n，满足
qkxi=σik
这里
σik=0,i≠k


1,i=k
取
qkx=∏ni=0i≠kx－xi
∏ni=0i≠k
xk－xi，k=0，1，2,…,n
于是就得到了fx的n次插值多项式
pnx=∑nk=0fxk∏ni=0i≠kx－xixk－xi
这就是Lagrange插值多项式。
则
∫bafxdx=∫bapnxdx=∫ba∑nk=0fxk∏ni=0i≠kx－xixk－xidx=∑nk=0fxk∫ba∏ni=0i≠kx－xixk－xidx
取等分步长h=b－an，则
xk=x0+khk=0，1，2,…,n
作变换
t=x－x0h
则
∫ba∏ni=0i≠kx－xixk－xidx=h∫n0∏nj=0j≠kt－jk－jdt=-1n－khn－k!k!∫n0∏nj=0j≠kt－jdtk=0，1，2,…,n
这种方法通常称为NewtonCotes公式。
根据上述公式，可编写函数程序interpoly，用插值法求解定积分的数值解： 




function I=interpoly(fun,a,b,n)

% INTERPOLY 插值法求解定积分的数值解

% 输入参数:

%——fun:被积函数

%——a:积分下限

%——b:积分上限

%——h:等分步长

% 输出参数:

%——I:求解定积分的数值解

% 调用说明:

% I=interpoly(fun,a,b,h):插值法求解定积分的数值解，待积分函数为fun,积分区间为［a,b］，等

% 分步长为h



fun=matlabFunction(fun);

X=linspace(a,b,n+1);

Y=feval(fun,X);

I=0;

for k=1:n

A=(-1)^(n-k)*h/(factorial(n-k)*factorial(k));

T=1;

syms t

for j=0:n

T=T*(t-j);

end

T=T/(t-k);

Ak=A*int(T,0,n);

I=I+eval(Ak*Y(k));

end




例59用插值法求定积分∫1041+x2dx的近似值。
解： 



syms s x

fun=4/(1+x^2)




fun = 4x2+1



format long

I=interpoly(fun,0,1,1/50)




I = 3.165111293937786
5.12动手实践
请用MATLAB求解下列问题。
1. 依据右取法(矩形右上角与曲线刚好相交)编写矩形法求积分函数。
2. 计算下面函数的导数： 

limx→0∫1cosxe－t2dtx2。
3. 计算无穷限区间的反常积分值∫+∞adxxp(a>0)。
4. 求解函数∫+∞11x2dx。
5. 对函数f(x)=11+25x2进行拉格朗日插值，观察插值函数是否收敛于原函数。