定积分的近似计算

1、数学软件与数学实验专题 定积分的近似计算数学软件与数学实验数学软件与数学实验2q 问题背景和实验目的专题专题 定积分的近似计算定积分的近似计算u 定积分计算的基本公式是牛顿莱布尼兹公式。但当被积函数的原函数不知道时，如何计算？这时就需要利用近似计算。特别是在许多实际应用中，被积函数甚至没有解析表达式，而是一条实验记录曲线，或一组离散的采样值，此时只能用近似方法计算定积分。u 本实验主要研究定积分的几种近似计算算法：矩形法、梯形法和抛物线法；蒙特卡洛随机投点法和蒙特卡洛样本均值法。同时介绍 Matlab 计算定积分的相关函数。数学软件与数学实验数学软件与数学实验3q 矩形法u 定积分的定义：专题

2、专题 定积分的近似计算定积分的近似计算badxxf)(,1iiixx0 x1x2x1nxnx 1x2x1ix ix nx,)(iixf ,1iiixxxiixxmax0limxnixni 1数学软件与数学实验数学软件与数学实验4矩形矩形法法1( )(), nbiiaif x dxfx n 充分大，x 充分小u 定积分的近似：l 通常我们取nxxx21nabh左点法右点法中点法l 点 可以任意选取，常见的取法有： 左端点 ，右端点 和中点 。,1iiixx1ixix2/ )(1iixx数学软件与数学实验数学软件与数学实验5()/ixhban, 1 2ixaihi, ,n步长节点u 右点法：11(

3、 )()()nnbiiaiif x dxfxhf xixu 中点法：1111( )()22(nnbiiiaiixxf x dxfxhfiixx111( )()()nnbiiaiif x dxfxhf x-1ixu 左点法：左点左点法、右点法和中点法法、右点法和中点法数学软件与数学实验数学软件与数学实验6解：矩形法举例矩形法举例h =1/100=0.01, xi = i*h, a=0, b=1, n=100 u 例：用不同的矩形法计算下面的定积分 ( 取 n=100 )， 并比较这三种方法的相对误差。1021xdxl 左点法：niiniixhxfhxdx1121110211 )(1 0.7878

4、9399673078l 右点法：1201()1niidxhf xx 0.78289399673078 0.78540024673078niiixxfhxdx11102)2(1l 中点法：(i = 0,1,2,.,100)数学软件与数学实验数学软件与数学实验711020arctan1dxxx40.78789399673078/4/40.003178l 理论值：l 左点法相对误差：0.78289399673078/4/40.0031880.78540024673078/4/4 -62.65310 u 误差分析矩形法举例矩形法举例l 右点法相对误差：l 中点法相对误差：不同的方法有不同的计算精度有没

5、有更好的近似计算定积分的方法有没有更好的近似计算定积分的方法 ?数学软件与数学实验数学软件与数学实验8)(xfab1ixixxyobadxxfS)(1S2SiSnSniibaSdxxfS1)(定积分几何意义定积分几何意义数学软件与数学实验数学软件与数学实验9iSiiiixyyS21nixfyii, 2 , 1 ),(u 曲边小梯形的面积可以由直边小梯形的面积来近似u 整个曲边梯形的面积：badxxfS)(iniiiniixyyS1112梯形法梯形法数学软件与数学实验数学软件与数学实验10u 如果我们 n 等分区间 a,b，即令：badxxfS)(niiiiniiiniiyyhxyyS11111

6、22nxxx21nabh则则22)(110nnbayyyyhdxxf梯形公式梯形公式梯形法梯形法梯形公式与中点公式有什么区别梯形公式与中点公式有什么区别 ?数学软件与数学实验数学软件与数学实验11解：101120122nndxyyhyyx 0.78539399673078u 例：用梯形法计算下面定积分 ( 取 n=100 )， 并计算相对误差梯形法举例梯形法举例1021xdxa=0, b=1, n=100, f (x) = 1/( 1+x2 ) h =1/100=0.01, xi = i*h, yi = f (xi) l 相对误差：0.78539399673078/4/4-65.30510 数

7、学软件与数学实验数学软件与数学实验12u 2n 等分区间 a,b ，得该直线用抛物线代替，计算精度是否会更好？11, , 0,1,22ibahxihinnu 计算每个节点上的函数值：nixfyii2 , 1 , 0 ),(抛物线法抛物线法u 在区间 x0, x2 上，用过以下三点),( ),( ),(222111000yxPyxPyxP的抛物线来近似原函数 f (x) 。数学软件与数学实验数学软件与数学实验13u 设过以上三点的抛物线方程为：则在则在区间 x0, x2 上，有上，有y = x2 + x + = p1(x) 2020)()(1xxxxdxxpdxxf20) (2xxdxxx20

8、2 3 23xxxxx20012 (4)6xxyyy012 (4)6bayyyn抛物线法抛物线法数学软件与数学实验数学软件与数学实验14u 同理可得：)4(6)( )4(6)(2122243222242nnnxxxxyyynabdxxfyyynabdxxfnnu 相加即得：2221222121( )( ) (4) 6iinbxaxiniiiif x dxf x dxbayyyn抛物线法抛物线法数学软件与数学实验数学软件与数学实验15u 整理后可得：)(2 )(46 )(2242123120nnnbayyyyyyyynabdxxf或辛普森 (Simpson) 公式抛物线法公式抛物线法抛物线法数学

9、软件与数学实验数学软件与数学实验16)(461123120102nnyyyyynabxdx 0.78539816339745)(22242nyyyu 例：用抛物线法计算下面定积分 ( 取 n=100 )， 并计算相对误差1021xdx解：a=0, b=1, n=100, yi = f (xi) = 1/( 1+xi2 ) 0.78539816339745/4/4-162.82710 l 相对误差：抛物线法抛物线法数学软件与数学实验数学软件与数学实验17u 梯形法：trapztrapz(x,y)x 为分割点（节点）组成的向量，y 为被积函数在节点上的函数值组成的向量。22)(110nnbayyy

10、ynabdxxf,x10nxxx01y (),(),()nf xf xf xq Matlab 近似计算定积分的相关函数Matlab 计算定积分函数介绍计算定积分函数介绍数学软件与数学实验数学软件与数学实验18前面的做法u 例：用梯形法计算下面定积分 ( 取 n=100) 1021xdx解：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函数1012120(/2/2)1nndxbayyyyyxntrapz(x,1./(1+x.2)trapz 举例举例数学软件与数学实验数学软

11、件与数学实验19quad(f,a,b,tol)f = f(x) 为被积函数，a,b 为积分区间，tol 为计算精度将自变量看成是将自变量看成是向量向量badxxf)(u 抛物线法：quadl 不用自己分割积分区间l 可以指定计算精度，若不指定，缺省精度是 10-6l 精度越高，函数运行的时间越长l 此处的函数 f 是数值形式，应该使用数组运算，即 点运算：.*，./ ，. ，. 注：抛物线法抛物线法数学软件与数学实验数学软件与数学实验201021xdx解：解： quad(1./(1+x.2),0,1) quad(1./(1+x.2),0,1,10e-10) quad(1./(1+x.2),0,

12、1,10e-16)函数表达式一定要用函数表达式一定要用 单引号单引号 括起来！括起来！涉及的运算一定要用涉及的运算一定要用 数组运算数组运算！u 例：用 quad 计算定积分：quad 举例举例数学软件与数学实验数学软件与数学实验21q 抛物线法计算二重积分： dblquaddblquad(f,a,b,c,d,tol)u tol 为计算精度，若不指定，则缺省精度为 10-6 badcdxdyyxf),(u f(x,y) 可以由 inline 定义，或通过一个函数句柄传递u a,b 是第一积分变量的积分区间，c,d 是第二积分变量 的积分区间按字母顺序，大写字母排在小写字母的前面二重积分的计算二

13、重积分的计算数学软件与数学实验数学软件与数学实验2221201(43)Ixyydxdy f=inline(4*x*y+3*y2); I=dblquad(f, -1,1,0,2)u f(x,y) 中关于第一自变量的运算是数组运算， 即把 x 看成是向量，y 看成是标量。 也可以全部采用数组运算例2：计算二重积分 20112)34(dxdyxxy dblquad(inline(4*x*y+3*x2),-1,1,0,2) dblquad(inline(4*x*y+3*x.2),-1,1,0,2)X例1：计算二重积分dblquad 举例举例数学软件与数学实验数学软件与数学实验23例：计算二重积分 20

14、112)34(dxdyxxy 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)q 被积函数 f (x,y) 的另一种定义方法：匿名函数dblquad 举例举例数学软件与数学实验数学软件与数学实验24int(f,a,b) 计算 f 关于默认自变量 的定积分，积分区间为a,b。int(f) 计算 f 关于默认自变量 的不定积分。int(f,v,a,b) 计算函数 f 关于自变量 v 的定积分，积分区间为 a, bint(f,v) 计算函数 f 关

15、于自变量 v 的不定积分badvvf)( )f v dvfindsym(f,1)q 符号积分： intint 符号积分符号积分数学软件与数学实验数学软件与数学实验25 syms x y; f=y*sin(x); int(f,x) int(f,y) int(f) int(a+b)ans=-y*cos(x)ans=1/2*y2*sin(x)ans=-y*cos(x)ans=a*b+1/2*b2u 例：指出下面各条语句的输出结果int 举例举例数学软件与数学实验数学软件与数学实验26u 例：用 int 函数计算定积分：解： syms x; f=1/(1+x2); int(f,x,0,1) f=sym

16、(1/(1+x2); int(f,x,0,1) int(1/(1+x2),x,0,1)或 int(1/(1+x2),0,1)1021xdx或或int 举例举例数学软件与数学实验数学软件与数学实验27double(a) 将将 a 转化为双精度型，若转化为双精度型，若 a 是字符，则取对应的是字符，则取对应的 ASCII 码码 a=3; double(a) double(a)例：ans = 3ans = 97其它相关函数其它相关函数数学软件与数学实验数学软件与数学实验28221dxex x=1:0.001:2; y=exp(x.(-2); trapz(x,y)l 梯形法：l 抛物线法： quad(

17、exp(x.(-2),1,2,10e-10)l 符号积分法： syms x int(exp(x(-2),x,1,2)例 1：用 Matlab 函数近似计算积分数值实验数值实验数学软件与数学实验数学软件与数学实验2920112)(dyyxdxl 抛物线法： dblquad(inline(x+y2),0,2,-1,1)l 符号积分法： f=int(x+y2,y,-1,1); int(f,0,2)数值实验数值实验例 2：用 Matlab 函数近似计算二重积分数学软件与数学实验数学软件与数学实验30q 用蒙特卡罗用蒙特卡罗 ( Monte Carlo ) 法近似计算法近似计算n=100000; %总模

18、拟次数总模拟次数m=0; %记录有利实验次数的变量初始化记录有利实验次数的变量初始化for i=1:n %开始模拟实验开始模拟实验 x=rand(1); %产生随机点产生随机点 y=rand(1); %产生随机点产生随机点 if ( 1/(x2 + 1) = y ) %判断试验点是否落在积分线下方判断试验点是否落在积分线下方 m=m+1; %有利实验次数加有利实验次数加1 endendfprintf(面积的近似值为：面积的近似值为：%fn, m/n)蒙特卡罗法近似计算定积分蒙特卡罗法近似计算定积分-投点法投点法1021xdx数学软件与数学实验数学软件与数学实验31蒙特卡罗法近似计算定积分蒙特卡

19、罗法近似计算定积分-样本平均值样本平均值法法设设g(x)是是( (a, b) )上的一个密度函数，改写上的一个密度函数，改写 dxxfba XgXfEdxxgxgxfbaniiixgxfnXgXfE1)()(1)()(n基本原理基本原理：对积分：对积分 其中，其中，X是服从是服从g(x)的随机变量可见，积分可以表示的随机变量可见，积分可以表示为为X的函数的期望。由矩法，若有的函数的期望。由矩法，若有n个来自个来自g( (x) )的观测值的观测值x1 1, , ,xn，则可给出，则可给出 的一个矩估计。的一个矩估计。数学软件与数学实验数学软件与数学实验32特别地，若特别地，若a，b有限，可取有限，可取 g(x) 为为 a,ba,b 上均匀分布此时上均匀分布此时, ,设设x1, xn是来自是来自U(a,b)的随机数，则的随机数，则 的一个估计为的一个估计为 niiniiixfnabxgxfn111 用上式来估计，和计算，随机数