数值微积分课件

MATLAB数学建模与仿真数值微积分2定积分计算的基本公式是牛顿－莱布尼兹公式。但当被积函数的原函数不知道时，如何计算？这时就需要利用近似计算。特别是在许多实际应用中，被积函数甚至没有解析表达式，而是一条实验记录曲线，或一组离散的采样值，此时只能用近似方法计算定积分。本实验主要研究定积分的三种近似计算算法：矩形法、梯形法和抛物线法。同时介绍Matlab

计算定积分的相关函数。

问题背景和实验目的定积分的近似计算1.极限和连续数列极限:>0,N>0，使当n>N时有xn-a<，则函数极限:如果当xx0时有f(x)A，则连续:如果当xx0时，有f(x)f(x0)

则称f(x)在x0连续。闭区间上连续函数必有最大值和最小值。预备知识：微积分当n=0得，微分中值定理

f(x)-f(x0)=f’()(x-x0)

其中是x0与x之间某个值Taylor公式:当f(x)在含有x0某个开区间内具有直到n+1阶的导数，3.多元函数微分学

设f(x,y)在点(x0,y0)附近有定义，当(x,y)以任何方式趋向于(x0,y0)时，f(x,y)趋向于一个确定的常数A，则若A=f(x0,y0),称f(x,y)在(x0,y0)点连续f(x,y)在点(x0,y0)的偏导数分别定义为4.积分

函数f(x)在区间[a,b]上的积分定义为其中a=x0<x1<…<xn=b,

xi=xi-xi-1,i(xi-1,xi),i=1,2,…,n若在[a,b]上,F’(x)=f(x),则二重积分定义为8矩形法梯形法抛物线法

数值积分的常见算法主要内容

Matlab

求积分函数数值积分函数：trapz、quad、dblquad符号积分函数：int10矩形法n

充分大，x

充分小

通常我们取左点法右点法中点法点可以任意选取，常见的取法有：

左端点,右端点和中点。定积分的近似：12矩形法举例例：用不同的矩形法计算下面的定积分(取n=100)，

并比较这三种方法的相对误差。左点法：右点法：中点法：解：h=1/n=0.01,xi=i*h,a=0,b=1,n=100(i=0,1,2,...,100)14定积分几何意义15

曲边小梯形的面积可以由直边小梯形的面积来近似整个曲边梯形的面积：梯形法16

如果我们n

等分区间[a,b]，即令：则==>梯形公式梯形法梯形公式与中点公式有什么区别

?

fuluB.m17解：==>例：用梯形法计算下面定积分(取n=100)，并计算相对误差梯形法举例a=0,b=1,n=100,f(x)=1/(1+x2)==>h=1/100=0.01,xi=i*h,yi=f(xi)

相对误差：182n

等分区间[a,b]，得用抛物线代替该直线，计算精度是否会更好？

计算每个节点上的函数值：抛物线法

在区间[x0,x2]上，用过以下三点的抛物线来近似原函数f(x)。20同理可得：相加即得：抛物线法21整理后可得：或辛卜生(Simpson)公式抛物线法公式抛物线法

fuluC.m23矩形法梯形法抛物线法

数值积分的常见算法Matlab函数

Matlab

求积分函数数值积分函数：trapz、quad、dblquad符号积分函数：int24矩形法总结

Matlab

数值积分函数：trapz、quad、dblquad梯形法抛物线法26前面的做法例：用梯形法计算下面定积分(取n=100)解：a=0,b=1,n=100,yi

=f(xi)=1/(1+xi2)>>

x=0:1/100:1;>>

y=1./(1+x.^2);>>

trapz(x,y)trapz函数trapz(x,1./(1+x.^2))trapz举例27quad(f,a,b,tol)f=f(x)为被积函数，[a,b]为积分区间，tol

为计算精度将自变量看成是向量不用自己分割积分区间可以指定计算精度，若不指定，缺省精度是10-6精度越高，函数运行的时间越长此处的函数

f是数值形式，应该使用数组运算，即：

.*

./

.\

.^

quad

quad抛物线法28解：>>

quad('1./(1+x.^2)',0,1)>>

quad('1./(1+x.^2)',0,1,1e-10)>>

quad('1./(1+x.^2)',0,1,1e-16)函数表达式一定要用单引号括起来！涉及的运算一定要用数组运算！例：用quad

计算定积分：quad举例30>>

f=inline('4*x*y+3*y^2');>>

I=dblquad(f,-1,1,0,2)

f

中关于第一自变量的运算是数组运算，即把x

看成是向量，y

看成是标量。也可以全部采用数组运算例：计算二重积分>>

dblquad(inline('4*x*y+3*x^2'),-1,1,0,2)>>

dblquad(inline('4*x*y+3*x.^2'),-1,1,0,2)X例：计算二重积分dblquad举例31例：计算二重积分>>

dblquad(@(x,y)4*x*y+3*x.^2,-1,1,0,2)指定x、y

分别是第一和第二积分变量>>

dblquad(inline('4*x*y+3*x.^2'),-1,1,0,2)被积函数f(x,y)

的另一种定义方法：匿名函数>>

dblquad(@(y,x)4*x*y+3*x.^2,-1,1,0,2)下面的命令运行结果和上面的一样吗？dblquad举例32int(f,a,b)

计算

f

关于默认自变量

的定积分，积分区间为[a,b]。int(f)

计算

f

关于默认自变量

的不定积分。int(f,v,a,b)

计算函数f

关于自变量v

的定积分，积分区间为[a,b]int(f,v)

计算函数

f

关于自变量

v

的不定积分findsym(f,1)int符号积分：int33例：用int

函数计算定积分：解：>>

symsx;>>

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

int(f,x,0,1)>>

f=sym('1/(1+x^2)');>>

int(f,'x',0,1)>>

int('1/(1+x^2)','x',0,1)或>>

int('1/(1+x^2)',0,1)或或int举例34double(a)将

a

转化为双精度型，若

a

是字符，则取对应的

ASCII码>>

a=3;>>

double(a)>>

double('a')例：ans=3ans=97相关函数35梯形法：抛物线法：符号积分法：用Matlab函数近似计算定积分作业36抛物线法：符号积分法：作业用Matlab函数近似计算二重积分Matlab中几个数值积分函数的比较和优缺点

一、Z=trapz(X,Y,dim)

梯形数值积分，通过已知参数x,y按dim维使用梯形公式进行积分。例

计算int(sin(x),0,pi)>>x=0:pi/100:2*pi;

>>y=sin(x);

>>z=trapz(x,y)%或者说使用z=pi/100*trapz(y)

z=1.0300e-017>>z=pi/100*trapz(y)二、[q,fcnt]=quad(fun,a,b,tol,trace,p1,p2...)

自适应simpson公式数值积分，适用于精度要求低，被积函数平滑性较差的数值积分。注意事项：

1.被积函数fun必须是函数句柄；

2.积分限[a,b]必须是有限的，因此不能为inf；

3.p1为其他需要传递的参数，一般是数值。可能警告：

1.'Minimumstepsizereached'

意味着子区间的长度与计算机舍入误差相当，无法继续计算了。原因可能是有不可积的奇点；

2.'Maximumfunctioncountexceeded'

意味着积分递归计算超过了10000次。原因可能是有不可积的奇点；

3.'InfiniteorNot-a-Numberfunctionvalueencountered'

意味着在积分计算时，区间内出现了浮点数溢出或者被零除。例2计算积分1/(x^3-2*x-p)，其中参数p=5，积分区间为[0,2]。>>F=@(x,n)1./(x.^3-2*x-n);

>>Q=quad(@(x)F(x,5),0,2)%或者使用quad(F,0,2,[],[],5)效果是一样的，%只是前者使用的函数嵌套。Q=

-0.4605>>quad(F,0,2,[],[],5)ans=

-0.4605三、[q,fcnt]=quadl(fun,a,b,tol,trace,p1,p2...)

自适应Lobatto数值积分，适用于精度要求高，被积函数曲线比较平滑的数值积分。注意事项：

同quad可能警告：

同quad例

计算积分1/(x^3-2*x-p)，其中参数p=5，积分区间为[0,2]。

>>F=@(x,p)1./(x.^3-2*x-p);

>>Q=quadl(F,0,2,[],[],5)%或者Q=quadl(@(x)F(x,5),0,2)。Q=

-0.4605四、[q,errbnd]=quadgk(fun,a,b,param1,val1,param2,val2,...)

自适应Gauss-Kronrod数值积分，适用于高精度和震荡数值积分，支持无穷区间，并且能够处理端点包含奇点的情况，同时还支持沿着不连续函数积分，复数域线性路径的围道积分法。注意事项：

1.积分限[a,b]可以是[-inf,inf]，但必须快速衰减；

2.被积函数在端点可以有奇点，如果区间内部有奇点，将以奇点区间划分成多个，也就是说奇点只能出现在端点上；

3.被积函数可以剧烈震荡；

4.可以计算不连续积分，此时需要用到'Waypoints'参数，'Waypoints'中的点必须严格单调；

5.可以计算围道积分，此时需要用到'Waypoints'参数，并且为复数，各点之间使用直线连接；

6.param,val为函数的其它控制参数，比如上面的'waypoints'就是，具体看帮助。例

计算有奇点积分int(exp(x)*log(x),0,1)。>>F=@(x)exp(x).*log(x);%奇点必须在端点上，否则请先进行区间划分。

>>Q=quadgk(F,0,1)Q=

-1.3179例

计算半无限震荡积分int(x^5*exp(-x)*sin(x),0,inf)。

>>F=@(x)x.^5.*exp(-x).*sin(x);

>>fplot(F,[0,100])%绘图，看看函数的图形

>>[q,errbnd]=quadgk(F,0,inf,'RelTol',1e-8,'AbsTol',1e-12)%积分限中可以有inf，但必须快速收敛q=-15.0000

errbnd=9.4386e-009五、q=dblquad(fun,xmin,xmax,ymin,ymax,tol,method)

矩形区域二重数值积分，一般区域二重积分参见NIT(数值积分工具箱)的quad2dggen函数。例9计算下面二重积分

>>F=@(x,y)y*sin(x)+x*cos(y);

>Q=dblquad(F,pi,2*pi,0,pi)Q=

-9.8696六q=triplequad(fun,xmin,xmax,ymin,ymax,zmin,zmax,tol,method)

长方体区域三重数值积分，注意此时没有一般区域的三重积分。例10计算下面三重积分。>>F=@(x,y,z)y*sin(x)+z*cos(x);

>>Q=triplequad(F,0,pi,0,1,-1,1)Q=

2.0000总结(1)quad：采用自适应变步长simpson方法，速度和精度都是最差的，建议不要使用；

(2)quad8：使用8阶Newton-Cotes算法，精度和速度均优于quad，但在目前版本下已被取消；

(3)quadl：采用lobbato算法，精度和速度均较好，建议全部使用该函数；

(4)quadg：NIT(数值积分)工具箱函数，效率最高，但该工具箱需要另外下载；

(5)quadgk：很有用的函数，功能在Matlab中最强大；

(6)quad2dggen：一般区域二重积分，效率很好，需要NIT支持；

(7)dblquad：长方形区域二重积分；

(8)triplequad：长方体区域三重积分；

1.导数、单调性与极值

当f’(x0)>0,函数在x0点附近是上升的，

f’(x0