matlab积分例子[共9页]

1、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 上可积时，函数返回一个定积

2、分结果。当 a,b 中有一个是 inf 时，函数返回一个广义积分。当 a,b 中有一个符号表达式时，函数返回一个符号函数。例：求函数 x2+y2+z2 的三重积分。内积分上下限都是函数，对 z 积分下限是 sqrt(x*y) ，积分上限是 x2*y ；对 y 积分下限是 sqrt(x) ，积分上限是 x2 ；对 x 的积分下限 1，上限是 2，求解如下：syms x y z % 定义符号变量f2=int(int(int(x2+y2+z2,z,sqrt(x*y),x2*y),y,sqrt(x),x2),x,1,2)%注意定积分的书写格式 f2 = 1610027357/6563700- 6072

3、064/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.数值积分的

4、实现方法基于变步长辛普生法， 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 即定积

5、分值， n 为被积函数的调用次数。例：求函数 exp(-x*x) 的定积分，积分下限为 0，积分上限为 1。fun=inline(exp(-x.*x),x); % 用内联函数定义被积函数 fname isim=quad(fun,0,1) % 辛普生法isim =0.746824180726425il=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.746

6、8或：format long gx=0:0.001:1; %x 向量，也可以是不等间距y=exp(-x.2); %y 向量，也可以不是由已知函数生成的向量s=trapz(x,y); % 求向量积分 s =0.746824071499185int 的积分可以是定积分，也可以是不定积分（即有没有积分上下限都可以积）可以得到解析的解，比如你对 x2 积分，得到的结果是1/3*x3 ，这是通过解析的方法来解的。如果 int(x2,x,1,2) 得到的结果是 7/3quad 是数值积分，它只能是定积分（就是有积分上下限的积分），它是通过 simpson 数值积分来求得的（并不是通过解析的方法得到 解析解

7、，再将上下限代入，而是用小梯形的面积求和得到的）。如果 f=inline(x.2);quad(f,1,2) 得到的结果是 2.333333 ，这个数并不是 7/3int 是符号解，无任何误差，唯一问题是计算速度； quad 是数值解，有计算精度限制，优点是总是能有一定的速度，即总能在一定时间内给出一个一定精度的解。from: 58.192.116.*对于 y=exp(-(x.2+x+1)/(1+x) ，被积函数之原函数无 封闭解析表达式 ,符号计算无法解题，这是符号计算有限性，结果如下： syms xy=exp(-(x.2+x+1)/(1+x) s=int(y,x,0,inf)y =exp(-

8、x2-x-1)/(1+x)warning: explicit integral could not be found.in at 58s =int(exp(-x2-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 % 所求定积分值或进行编程，积分上限人工输入，程序

9、如下：%表达式保存为函数文件 function y=fxy(x)y=exp(-(x.2+x+1)./(1+x); % save fxy.m% main - 主程序clear,clch=.001;p=0;a=0;r=input( 请输入积分上限， r=) while a rp=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=.30676016860-2: int=.45996331590-5: int=.5583068217

10、0-10: int=.56409289750-100: int=.56413460550-1000: int=.5641346055from: 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,

11、y,z了，可以求三重积分了。 完整函数程序：function fn(n1,n2,n3)if n1=0 e1=1;else if n1 0e1=2;end endif n2=0e2=1;else if n2 0e2=2; end endif n3=0e3=1;else if n3 0e3=2; end endf=(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 文件，然

12、后运行： fn(1,1,1)s =866.9655from: 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

13、,n2,n3 就会以常数看待，未知数就只有 x,y,z 了。完整函数程序：function fn(n1,n2,n3)if n1=0 e1=1;else if n1 0e1=2; end endif n2=0e2=1;else if n2 0e2=2;endendif n3=0e3=1;else if n3 0e3=2; end endf=(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.965

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

15、、物理学、以及工程应用里较为常见的特殊函数 modified bessel function of thefirst kind （数学符号记作：， matlab 中对应的函数是 besseli ），介绍如何克服积分结果为 nan 的问题。本帖的方法虽然只是针对这 一种特殊函数，但解决问题的思路适用于这一类数值计算问题，希 望能对有兴趣的人有帮助！最后，在分享本案例的过程中，我还会顺便指出 matlab r2014a 符号积分系统的一个计算错误，应该算是一个 bug 。特殊函数积分案例先给出我要计算的积分：首先，我们有必要先对的函数特性做一个基本了解（参考 ）。可以看出：是一个单调增函数从 (2

16、) 式可以得到，类比，我们可以粗略的认为， 近似服从指数渐进性（严格的说，当 ，）。所以， 是一个近似服从指数增长的函数。要保证这样一个无穷积分收敛，被积函数里必须要有一个衰减速度极快的因子，该因子衰减速度远快于的增长的速率，使得被积函数整体是一个快速衰减的函数，积分才有可能收敛。所以，（ 1）式的被积函数表达式里有一个快速衰减项，该项使得被积函数整体快速衰减（注：在面前，平方项的递增影响可以忽略不记）。下面，我们来看看直接用 matlab integral 函数计算的结果（ matlab版本低于 2012a 的话可以用 quadgk 代替）。 integral(x)x.2.*exp(-x.2

17、).*besseli(0,x),0,inf) 复制代码 warning:infinite ornot-a-number value encountered. infunfunprivateintegralcalciteratescalarvalued at 349in funfunprivateintegralcalc vadaptat132in funfunprivateintegralcalc at 83in integral at 88ans = nan 可以看到，直接计算的结果为 nan 。产生 nan 的原因产生 nan 的原因其实很简单，因为按 matlab 的运算规则，一旦计算过

18、程中出现了 nan ，最终结果就为 nan 。所以，这里一定是因为计算被积函数的函数值时得到了 nan 。简单的说，当 x 很大时（如1000 ），快速递增单独计算的结果为 inf ，而快速衰减项的结果为 0，二者相乘， matlab 里会得到 inf*0 = nan 。为了详细说明问题，我们先来计算一下被积函数当 x = 1000 的取值。直接用 matlab 数值计算f=(x)x2.*exp(-x2).*besseli(0,x) f(1000) 复制代码 ans = nan数值计算得到结果为 nan 。这正是因为单独的两项 besseli(0,x) 和exp(-x2) 分别计算时得到了 i

19、nf 和 0，如下所示： exp(-10002) 复制代码 ans =0besseli(0,1000) 复制代码 ans = inf二者再相乘得到 nan 。既然在计算函数值的过程中就出现了 nan ，难怪积分结果也是 nan 。我们再看看符号计算的结果： vpa(sym(10002*exp(-10002)*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) 复制代