matlab积分例子

matlab积分例子

【篇一：matlab积分例子】

符号积分由函数int来实现。该函数的一般调用格式为：

int(s)：没有指定积分变量和积分阶数时，系统按findsym函数指示的默认变量对被积函数或符号表达式s求不定积分；

int(s,v)：以v为自变量，对被积函数或符号表达式s求不定积分；int(s,v,a,b)：求定积分运算。a,b分别表示定积分的下限和上限。该函数求被积函数在区间[a,b]上的定积分。a和b可以是两个具体的数，也可以是一个符号表达式，还可以是无穷(inf)。当函数f关于变量x在闭区间[a,b]上可积时，函数返回一个定积分结果。当a,b中有一个是inf时，函数返回一个广义积分。当a,b中有一个符号表达式时，函数返回一个符号函数。

例：

求函数x^2+y^2+z^2的三重积分。内积分上下限都是函数，对z积分下限是sqrt(x*y)，积分上限是x^2*y；对y积分下限是sqrt(x)，积分上限是x^2；对x的积分下限1，上限是2，求解如下：

syms x y z %定义符号变量

f2=int(int(int(x^2+y^2+z^2,z,sqrt(x*y),x^2*y),y,sqrt(x),x^2),x,1,2) %注意定积分的书写格式

f2 =

1610027357/6563700-

6072064/348075*2^(1/2)+14912/4641*2^(1/4)+64/225*2^(3/4) %给出有理数解

vf2=vpa(f2) %给出默认精度的数值解

vf2 =

224.92153573331143159790710032805

二、数值积分

1.数值积分基本原理

求解定积分的数值方法多种多样，如简单的梯形法、辛普生(simpson)?法、牛顿－柯特斯(newton-cotes)法等都是经常采用的方法。它们的基本思想都是将整个积分区间[a,b]分成n个子区间[xi,xi+1]，i=1,2, ,n，其中x1=a，xn+1=b。这样求定积分问题就分解为求和问题。

2.数值积分的实现方法

基于变步长辛普生法，matlab给出了quad函数来求定积分。该函

数的调用格式为：

[i,n]=quad(fname,a,b,tol,trace)

基于变步长、牛顿－柯特斯(newton-cotes)法，matlab给出了

quadl函数来求定积分。该函数的调用格式为：

[i,n]=quadl(fname,a,b,tol,trace)

其中fname是被积函数名。a和b分别是定积分的下限和上限。tol

用来控制积分精度，缺省时取tol=0.001。trace控制是否展现积分

过程，若取非0则展现积分过程，取0则不展现，缺省时取trace=0。返回参数i即定积分值，n为被积函数的调用次数。

例：

求函数exp(-x*x)的定积分，积分下限为0，积分上限为1。

fun=inline(exp(-x.*x),x); %用内联函数定义被积函数fname

isim=quad(fun,0,1) %辛普生法

isim =

0.746824180726425

il=quadl(fun,0,1) %牛顿－柯特斯法

il =

0.746824133988447

三、梯形法求向量积分

trapz(x,y) 梯形法沿列方向求函数y关于自变量x的积分(向量形式，数值方法)。

d=0.001;

x=0:d:1;

s=d*trapz(exp(-x.^2))

s=

0.7468

或：

format long g

x=0:0.001:1; %x向量，也可以是不等间距

y=exp(-x.^2); %y向量，也可以不是由已知函数生成的向量

s=trapz(x,y); %求向量积分

s =

0.746824071499185

int的积分可以是定积分，也可以是不定积分（即有没有积分上下限

都可以积）可以得到解析的解，比如你对x^2积分，得到的结果是

1/3*x^3，这是通过解析的方法来解的。如果int(x^2,x,1,2)得到的结

果是7/3

quad是数值积分，它只能是定积分（就是有积分上下限的积分），

它是通过simpson数值积分来求得的（并不是通过解析的方法得到

解析解，再将上下限代入，而是用小梯形的面积求和得到的）。如

果f=inline(x.^2);quad(f,1,2)得到的结果是2.333333，这个数并不

是7/3

int是符号解，无任何误差，唯一问题是计算速度；quad是数值解，有计算精度限制，优点是总是能有一定的速度，即总能在一定时间

内给出一个一定精度的解。

[from: 58.192.116.*]

对于y=exp(-(x.^2+x+1)/(1+x))，被积函数之原函数无封闭解析表达式 ,符号计算无法解题，这是符号计算有限性，结果如下：

syms x

y=exp(-(x.^2+x+1)/(1+x))

s=int(y,x,0,inf)

y =

exp((-x^2-x-1)/(1+x))

warning: explicit integral could not be found.

in sym.int at 58

s =

int(exp((-x^2-x-1)/(1+x)),x = 0 .. inf)

只有通过数值计算解法

dx=0.05; %采样间隔

x=0:dx:1000; %数值计算适合于有限区间上,取有限个采样点，只要

终值足够大，精度不受影响

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

s=dx*cumtrapz(y); %计算区间内曲线下图形面积,为小矩形面积累

加得

s(end)

ans =

0.5641 %所求定积分值

或进行编程，积分上限人工输入，程序如下：

%表达式保存为函数文件

function y=fxy(x)

y=exp(-(x.^2+x+1)./(1+x)); % save fxy.m

% main --------主程序

clear,clc

h=.001;p=0;a=0;

r=input(请输入积分上限，r=)

while a r

p=p+(fxy(a)+fxy(a+h))*h/2;

a=a+h;

end

p=vpa(p,10)

运行主程序后得到结果：

请输入积分上限，r=1000

r =

1000

p =

.5641346055

其它结果如下：

0-1: int=.3067601686

0-2: int=.4599633159

0-5: int=.5583068217

0-10: int=.5640928975

0-100: int=.5641346055

0-1000: int=.5641346055

[from: 211.65.33.*]

在积分函数

中,sqrt(e1*e2*e3)*cos(n1*pi*x/12).*cos(n2*pi*y/11).*cos(n3*pi*z/ 9);已知变量e1,e2,e3,n1,n2,n3通过函数参数输入,如果直接用

inline或字符串的形式，则表达式中

的未知数有9个，分别是e1,e2,e3,n1,n2,n3,x,y,z。而用匿名函数时，已知变量e1,e2,e3,n1,n2,n3就会以常数看待，未知数就只有x,y,z 了，可以求三重积分了。

完整函数程序：

function fn(n1,n2,n3)

if n1==0

e1=1;

else if n1 0

e1=2;

end

end

if n2==0

else if n2 0

e2=2;

end

end

if n3==0

e3=1;

else if n3 0

e3=2;

end

end

f=@(x,y,z)sqrt(e1*e2*e3)*cos(n1*pi*x/12).*cos(n2*pi*y/11).*cos( n3*pi*z/9);

s=triplequad(f,-6,6,-5.5,5.5,-4.5,4.5) %求三重数值积分

将以上保存为fn.m程序文件，即m文件，然后运行：

fn(1,1,1)

s =

866.9655

[from: 211.65.33.*]

三重积分请用三重积分函数triplequad，与三个积分上下限对应，即x=triplequad(f,-6,6,-5.5,5.5,-4.5,4.5)

其中被积函数f用匿名函数来表达，即

f=@(x,y,z)sqrt(e1*e2*e3)*cos(n1*pi*x/12).*cos(n2*pi*y/11).*cos( n3*pi*z/9);

如果直接用inline或字符串的形式，则表达式中的未知数有9个，分别是e1,e2,e3,n1,n2,n3,x,y,z。而用匿名函数时，已知变量

e1,e2,e3,n1,n2,n3就会以常数看待，未知数就只有x,y,z了。

完整函数程序：

function fn(n1,n2,n3)

if n1==0

e1=1;

else if n1 0

e1=2;

end

end

if n2==0

e2=1;

else if n2 0

e2=2;

end

if n3==0

e3=1;

else if n3 0

e3=2;

end

end

f=@(x,y,z)sqrt(e1*e2*e3)*cos(n1*pi*x/12).*cos(n2*pi*y/11).*cos( n3*pi*z/9);

x=triplequad(f,-6,6,-5.5,5.5,-4.5,4.5)

fn(1,1,1)

x =

866.9655

[from: 58.192.116.*]

【篇二：matlab积分例子】

前言

本帖将通过一个典型的数值积分案例来介绍如何在 matlab 中有效地计算特殊函数（special functions）的数值积分。特殊函数的数值积分（尤其是无穷积分）比较常见的一个问题是：积分在满足收敛条件下积分结果为 nan。如果我们使用的是 integral 或 quadgk，这一问题通常会以警告的形式给出: “warning:infinite or not-a-number value encountered,” 伴随该警告的恶果是，数值积分的结果为 nan。

为了直观地说明问题，本帖将针对一个在数学、物理学、以及工程应用里较为常见的特殊函数—— modified bessel function of the first kind（数学符号记作：，matlab 中对应的函数是 besseli），介绍如何克服积分结果为 nan 的问题。本帖的方法虽然只是针对这一种特殊函数，但解决问题的思路适用于这一类数值计算问题，希望能对有兴趣的人有帮助！

最后，在分享本案例的过程中，我还会顺便指出 matlab r2014a 符号积分系统的一个计算错误，应该算是一个 bug。

特殊函数积分案例

先给出我要计算的积分：

首先，我们有必要先对的函数特性做一个基本了解（参考）。

可以看出：

是一个单调增函数

从 (2) 式可以得到，，类比，我们可以粗略的认为，近似服从指数

渐进性（严格的说，当，）。

所以，是一个近似服从指数增长的函数。要保证这样一个无穷积分

收敛，被积函数里必须要有一个衰减速度极快的因子，该因子衰减

速度远快于的增长的速率，使得被积函数整体是一个快速衰减的函数，积分才有可能收敛。所以，（1）式的被积函数表达式里有一个

快速衰减项，该项使得被积函数整体快速衰减（注：在面前，平方

项的递增影响可以忽略不记）。

下面，我们来看看直接用 matlab integral 函数计算的结果（matlab 版本低于2012a的话可以用 quadgk 代替）。integral(@(x)

x.^2.*exp(-x.^2).*besseli(0,x),0,inf)复制代码warning:infinite or not-a-number value encountered. infunfun\private\integralcalc iteratescalarvalued at 349in funfun\private\integralcalc vadapt

at132in funfun\private\integralcalc at 83in integral at 88

ans = nan

可以看到，直接计算的结果为nan。

产生nan的原因

产生nan的原因其实很简单，因为按matlab的运算规则，一旦计算过程中出现了nan，最终结果就为nan。所以，这里一定是因为计

算被积函数的函数值时得到了nan。简单的说，当 x 很大时（如1000），快速递增单独计算的结果为inf，而快速衰减项的结果为0，二者相乘，matlab 里会得到 inf*0 = nan。为了详细说明问题，我们

先来计算一下被积函数当 x = 1000 的取值。直接用matlab数值计

算

f=@(x)x^2.*exp(-x^2).*besseli(0,x)

f(1000)复制代码ans = nan

数值计算得到结果为nan。这正是因为单独的两项 besseli(0,x) 和exp(-x^2) 分别计算时得到了 inf 和 0，如下所示：exp(-1000^2)复

制代码ans =

0besseli(0,1000)复制代码ans =

inf

二者再相乘得到 nan。既然在计算函数值的过程中就出现了nan，

难怪积分结果也是nan。

我们再看看符号计算的结果：vpa(sym(1000^2*exp(-

1000^2)*besseli(0, 1000)),4)复制代码ans =

8.195e-433857

按照双精度数值计算的精度，这个结果在matlab数值计算里可以认为是0。所以，数值计算中如果能将这个结果作为0处理，也就能避免积分结果为nan了。

解决nan的方法

这里介绍的方法适用于一类具备共同特征的特殊函数积分：被积函数里有一个无上界的快速递增项和一个快速衰减项，二者共同作用的结果使得被积函数在积分区间上或者积分变量 x 取值较大时快速衰减，从而积分收敛。

最简单的解决办法：既然问题出现在大 x 下，可以适当减小积分上限，将上限 inf 换成一个有限的数，该数必须满足两个条件：1）必须足够大使得被积函数整体的取值很小（近似为0），2）又不能太大而导致数值计算里出现 inf*0 = nan。我们可以尝试将积分上限设置为100：integral(@(x) x.^2.*exp(-x.^2).*besseli(0,x),0,100)复制代码ans = 0.637956643215776

比较通用的解决办法：有时候直接将被积函数里快速衰减因子和快速增长因子的乘积重新定义一下可能会更方便。针对之前的案例，我们可以考虑重新定exp(-x.^2).*besseli(0,x)，使得新的定义里 x 很大时直接返回结果0。为此，我们可以先大致估算一下

besseli(0,x)在 x 700 时得到inf，故我们可以定义：function y = expbesseli0(x)

y = zeros(size(x));

id = x

y(id) = exp(-x(id).^2).*besseli(0,x(id));复制代码然后，重新使用下列命令计算积分integral(@(x) x.^2.*expbesseli0(x),0,inf)复制代码ans = 0.637956643215780

注意，上述代码里，积分上限依然是inf，因为我们已经重定义了expbesseli0。

matlab/mupad符号积分的错误

为了验证上述解决办法的合法性，我特意使用了 mupad 符号积分来验证。出乎意料的是，mupad 符号积分的结果完全错误！

syms x

int(x^2*exp(-x^2)*besseli(0,x),x,0,inf)复制代码ans =-

(3*exp(1/4))/8

符号积分的结果为负数，这显然是不可能的。被积函数在积分区间上是非负的，积分结果不可能是负的，这足以说明 mupad 符号积分计算错误。但为了进一步验证，我们可以看看 maple 的计算结果：