试验10 数值微积分

一、实验目的与要求

1．观察从定义出发和用微分方程法计算定积分的原理及两者的联系。

2．体会微积分基本定理New-Leibniz公式在沟通不定积分和定积分中所起的作用。

3．熟练掌握和使用Matlab提供的几个数值积分函数进行计算。

二、问题描述

选取一个特殊函数，如f(x)sinx,0x,计算定积分sinxdx。

0 这个问题可以直接算出sinx的原函数，从而计算出定积分，易知sinxdx＝

02。对于更一般情况，如果原函数得不到，或者计算过程太复杂，如何去计算定积分呢。

三、问题分析

上例中，虽然可以直接算出sinx的原函数，从而计算出定积分，但这样做失去了一般意义，本实验希望用解微分方程的欧拉折线法和Matlab提供的数值积分函数来近似计算定积分，在实验过程中利用Matlab图形处理功能将结果可视化以便观察。

四、背景知识介绍

在工程问题和科学实验中，常常需要计算定积分，例如，力学和电力学中功和功率的计算，电流和电压平均值的计算以及一些几何图形的面积，体积和弧长的计算等，另外微积分方程的求解也是以积分计算为基础的。

（1） 直接以定义出发，取近似和的极限，这种方法原理简单，但计算量大，在

高等数学中不是重点内容，然而数值积分的各种算法却是基于定积分定义建立起来的。

（2） 用不定积分计算定积分，不定积分是求导的逆运算，而定积分是连续变量

的求和（如计算曲边梯形的面积）。表面上看是两个完全不同的概念，然而却通过牛顿－莱布尼兹公式联系在一起了，甚至两者的符号也是类似的，牛顿－莱布尼兹公式被称为微积分学基本定理。 微积分基本定理（牛顿－莱布尼兹公式）： 如果F(x)是f(x)的原函数，即

F(x)f(x)，则有f(x)dxF(b)F(a)成立

a'b因此理论上为了计算定积分f(x)dx，只需要求出f(x)的一个原函数，这可

ab通过求f(x)不定积分来完成。反之，如果将定积分的积分上限看出变量，便得

到了上限积分函数F(x)xaf(t)dt，它又是f(x)的一个原函数。这意味着，对于

满足某些条件的函数f(x)，不论其不定积分f(x)dx是否存在解析解，我们总可以算出（至少在数值上）它的不定积分f(x)dx=f(x)dxC，这在理论分

ax析和数值计算中都是很便利的。在高等数学中求不定积分是重点学习内容之一，我们总是期望求出不定积分的封闭解来，但在实践中能够算出封闭解的积分并不多，所以数值积分是非常有用的工具。牛顿-莱布尼兹公式不愧为微积分的“基本定理”。 （3）利用解微积分方程的方法计算定积分。设函数yy(x)，axb是连续可微的，它表述了平面上的一条曲线。如果y的导数满足

dydxf(x)

（1.1）

且曲线yy(x)过a,y0点，我们的问题是如何由式（1.1）式确定曲线

yy(x)？方程（1.1）又和定积分f(t)dtab有什么联系呢？

由牛顿-莱布尼兹公式有yy(x)y0ybxaf(t)dt，令初始值y00，则

baf(t)dt就是所需的积分。这就是说，如果算出了过初始点（a，0）的曲线

byy(x)，则定积分f(t)dta就是曲线在右端点处的值y(b)。方程（1.1）是一个

最简单的微分方程，把它右边的f(x)改为f(x,y)，就得到一般的一阶常微分方程的初值问题：

dyf(x,y)dxy(a)y0 （1.2）

相对于方程（1.1），求方程（1.2）的解析解增加了很大的难度，但是从数值

计算的角度来看并没有本质的困难。 五、实验过程 1．微分方程法

（1） 考虑最简单的微分方程

yy(x)，0xdydxsin(x),0x，它确定了一族解曲线

，观察这族解曲线。

运行观察程序exp1-1.m,得到如图1－1。

图1－1

（2）确定方程

dydxsin(x),0x，过(0,0)的解曲线yy(x)。观察步骤如下：

步骤1. 将[0,Pi]区间n=4等分（共有5个分点）

【 n=4;n1=n+1;x=linspace(0,pi,n1);f=myfun2_2(x);[0:n;x;f] 】

步骤2. 在第一个子区间[x(1),x(2)]上，令初始点y(1)=y0=0, 按（4.4）式画出相应的折线段，折线段的起点为[x(1),y(1)],终点记为[x(2),y(2)],其中 【 y(1)=0; y(2)=y(1)+f(1)*(x(2)-x(1));

P_intial=[x(1),y(1)],P_final=[x(2),y(2)], 】

步骤3. 观察在第二子区间[x(2),x(3)]中的情况，按任意键后观察程序继续运算。 【 y(3)=y(2)+f(2)*(x(3)-x(2));

P_intial=[x(2),y(3)], P_final=[x(3),y(3)] 】 步骤4。观察增量Dy的意义 【 Dy=y(3)-y(2) 】 【 DS=f(2)*(x(3)-x(2)) 】

步骤5 . 按回车键将程序运行完毕。

运行观察程序exp1-2.m，在程序中设置了Pause语句，运行后只显示第一个子区间计算结果，按回车键则显示下一子区间的结果，如此下去直至完成计算。图1－2显示了所有区间上的情况。

图1－2

以下是观察程序及其说明：

观察程序exp1_1.m （观察方程的解曲线族，图1-1） 【 clf,clear,a=0;b=pi;c=0;d=2.2;n=20;

[X,Y]=meshgrid(linspace(a,b,n),linspace(c,d,n));

z=X.*Y; Fx=cos(atan(sin(X)));Fy=sin(atan(sin(X))); quiver(X,Y,Fx,Fy,0.5),hold on,axis([a,b,c,d]) [x,y]=ode45('exp1_1f',[0,pi],0);

[x1,y1]=ode45(' exp1_1f ',[0,pi],0.2); [x2,y2]=ode45(' exp1_1f ',[0,pi],0.4); [x3,y3]=ode45(' exp1_1f ',[0,pi],0.6);

plot(x,y,'r.-',x1,y1,'r.-',x2,y2,'r.-',x3,y3,'r.-'), xlabel('x') ,ylabel('y') 】 exp1_1f.m

【 function dy=exp1_1f(x,y)

dy=sin(x); 】 . 观察程序exp1_2.m （图1-2）

【 clf, a=0;b=4*pi;n=31;

x=linspace(a,b,n); y=zeros(1,n); y(1)=0; s(1)=0;

for i=1:n-1 dy=myfun1_2(x(i)); y(i+1)=y(i)+dy*(x(i+1)-x(i)); hh(i)=myfun1_2(x(i));

s(i+1)=s(i)+hh(i)*(x(i+1)-x(i)); end

a1=0.9*a;b1=1.1*b; % 设置坐标范围。 ymin=min(y');ymax=max(y'); ymin1=ymin*0.9;ymax1=ymax*1.1;

subplot(2,1,1), %在第一幅图中画垂线,和原函数。

fplot('1-cos(x)',[0,pi],'r:');axis([a1 b1 ymin1 ymax1]),hold on, plot([x(1) x(1)],[ymin ymax]); subplot(2,1,2), %在第二幅图中画被积函数图象。

fplot('myfun2_2',[a b],'-k');hold on,axis([a1 b1 ymin1 ymax1])

for i=1:n-1 %在第I个子区间上计算。

subplot(2,1,1),

plot([x(i+1) x(i+1)],[ymin ymax]); %在各分点画竖线。 plot([x(i) x(i+1)],[y(i),y(i)]); %画水平线。 h=plot([x(i+1) x(i+1)],[y(i),y(i+1)]);%画表示增量的垂线。 set(h,'linewidth',3,'color','b'); %设置线宽和颜色。

h1=plot([x(i) x(i+1)],[y(i) y(i+1)],'.-'); %画折线，设置图形属性 。

set(h1,'linewidth',2,'markersize',18); subplot(2,1,2),

xx=[x(i) x(i) x(i+1) x(i+1)]; yy=[0 hh(i) hh(i) 0]; %计算矩形顶点坐标。

patch(xx,yy,[0.7 0.7 0.7]); %在第二幅图中画矩形块并填充颜色。

plot([x(i+1) x(i+1)],[ymin ymax]);

plot([x(1) x(1)],[s(i),s(i+1)],'r','linewidth',5); %在y轴上画面积增量线段。

plot(x(1),y(i+1),'r.','Markersize',18);

subplot(2,1,1),plot([x(i+1) x(i+1)],[ymin ymax]);

h=plot([x(1) x(1)],[s(i),s(i+1)]); %画相应的辅助线段。

set(h,'linestyle','-','linewidth',5,'color','red'); plot(x(1),y(i+1),'r.','Markersize',18); plot([x(1) x(i+1)],[s(i+1) y(i+1)],'r--')

pause %暂停，n大时应该去掉。 end 】 myfun1_2.m (选择其它的函数进行实验，可修改此程序) 【 function y=myfun1_2(x)

sin(x); %y=sin(x)./(x); 】

2．利用Matlab积分命令数值求解定积分

（1） 矩形求积指令 Cumsum 观察cumsum指令的功能。

【 x=[1 1 1 1 1]；I=cumsum(x) 】 观察：用累积和命令cumsum计算积分。 【 x=linspace(0,pi,50); y=sin(x);

T=cumsum(y)*pi/(50-1);I=T(50)

（2） 梯形公式指令trapz

观察梯形公式计算的效果。

【 x=linspace(0,pi,50); y=sin(x);T=trapz(x,y) 】 （3） 辛普森公式指令quad

观察辛普森积分公式的计算效果。

【 I1=quad('sin',0,pi), I2=quad8('sin',0,pi), 】

（4） 解一阶微分方程组指令 ode23, ode45(龙格－库塔法) 观察：用解常微分方程的方法计算积分0sin(x)dx。

【 y0=0;[x,y]=ode23('myfun2_2',[0,pi],y0);s=length(y);y(s) 】

【 y0=0;[x,y]=ode45('myfun2_2',[0,pi],y0);s=length(y);y(s) 】

六、结论与应用

1．观察图1-1最下面的解曲线y(x)，它过（0，0）点，而且y()2恰等于积分的值，这是正确的结果，由式（1.2）式知这条曲线为

y0sinxdx1coxs,0x，所以定积分

dydxsin(x),0x0sinxdx恰等于

y()1cos2，则确定方程的过（0，0）点的解曲线y(x)

2．实验中采用了适合微分方程数值计算的欧拉折线法，其原理如下： 1将[0,Pi]n等分：x0x1xixi1xn，计算分点和相应的○

函数值xi,fi,i0,,n1

2在每个子区间xi,xi1上，以折线段y(x)yifi(xxi)代替解曲线y(x)。 ○

其中x00,y00，折线段的起点为xi,fi，终点为xi1,fi1。这样就近似地确定了这条解曲线。

运行观察程序exp1-2.m可对上面的原理进行观察。积分的近似值是图1-2下图的小矩形面积之和，也是上图中各增量Dy之和。

进一步可以进行下面的观察：

1) 增加分割数，例如取n＝50，观察运行结果。 2) 在0,4中观察。

3) 选择其他的函数进行观察。

3．用积分定义计算数值积分是十分自然的，当分割较为细密时，用离散和式代替连续和（积分）可得到定积分的近似值。Matlab数值积分函数就是基于这样的原理进行计算的。

本实验所采用的数值方法适合所有的微分、积分模型，其应用范围十分广泛。

七、练习

我辑私雷达发现，距离d处有一走私船正以匀速a沿直线行驶，缉私舰立即以最大速度（匀速v）追赶。若用雷达进行跟踪，保持船的瞬时速度方向始终指向走私船，则辑私舰的运动轨迹是怎么的？是否能够追上走私船？如果能追上，需要多长时间？ y

M0

M(x,y) d  S0Sx

图1－3

1、试用Eular法进行数值模拟，并画出相应的追击轨线，进行适当的结果分析。 2、追线问题化为一个常微分方程,用ode45指令解方程并画出追击曲线。