计算声学第七章 数值微积分

1、第七章第七章 数值微积分数值微积分球面波在水平界面上的反射声场求解球面波在水平界面上的反射声场求解 球面波在平面界面上的反射和折射问题处理起来比较复杂，可以采用波动声学的方法来进行处理。首先对点声源作空间傅立叶变换，即把点声源辐射的球面波分解成按照不同方向传播的平面波，表达式为平面波的双积分形式；然后应用平面波在分界面上的反射理论，计算得到各平面波分量的反射波；最后把所有反射波分量通过积分进行叠加，从而得到球面波在平面界面上的反射声场。 通过坐标变换和化简，得到反射声场的表达式为 jjrdeJVjkpz200sin第七章第七章 数值微积分数值微积分 求函数 在区间 上的定积分是微积分学中的基本

2、问题之一，也是实际问题中经常遇到的计算问题。 在微积分中，牛顿-莱布尼兹公式是计算定积分的一种有效工具。对于上面的定积分，如果 在区间 上连续，并且 的原函数为 ，则可计算定积分)(xf,ba badxxffI)()(xf,ba)(xf)(xF )()()(aFbFdxxffIba第七章第七章 数值微积分数值微积分l实际计算积分时可能存在的困难：实际计算积分时可能存在的困难：(1) 本身形式复杂，求原函数更为困难，如(2) 虽然形式简单，但是其原函数不能用初等函数表示成有限的形式，如(3) 的原函数能用初等函数表示，但较为复杂，如 )(xf)(xf)(xfcbxaxxf2)(,ln1)(xxf

3、,11)(4xxf,2xe,sin2xxxsin 12arctan12arctan2211212ln24122xxxxxxxF第七章第七章 数值微积分数值微积分(4) 本身没有具体的解析表达式，仅有表格或图形给出的实验测量数据点。 以上几种情况下就无法应用牛顿-莱布尼兹公式，所以，需要研究计算定积分的近似近似方法数值积分法数值积分法；在微分学中，函数 的导数是通过极限定义的，如果函数以一组观测数据的形式给出，或者表达式过于复杂时，就需要用数值方法求解函数的导数或微分导数或微分。 )(xf)(xf第七章第七章 数值微积分数值微积分l数值积分的基本思想数值积分的基本思想 找到一个符合精度要求的简单

4、函数 代替原来的函数 ，得到由于简单函数 积分容易计算，所以积分值容易得到。 采用数值方法计算积分值，首先需要把积分区间细分，在每一个小区间内用简单函数来代替原来的复杂函数进行积分计算。本章学习利用代数插值多项式代数插值多项式来代替被积函数进行积分计算。 xP xf babadxxPdxxf xP xf第七章第七章 数值微积分数值微积分 对于积分其几何意义为由 所围成曲边梯形的面积。计算曲边梯形面积比较困难，因为有一条曲边 badxxffI)( xfyybxax, 0, xfy 第七章第七章 数值微积分数值微积分 如果用直线、抛物线等代替曲边梯形的曲边，则面积计算变得容易。用容易计算面积的图形

5、来代替曲边梯形，就可以求出曲边梯形面积的近似值，从而得到积分的近似值。 按照这种思想，可以构造一些求积分值的近似公式，如梯形公式：中矩形公式： bfafabdxxfba2 2bafabdxxfba第七章第七章 数值微积分数值微积分 当函数 已知时，讨论定积分的计算问题为了避开求原函数的困难，通过被积函数 的值来求解积分值。由积分中值定理积分中值定理，对连续函数 ，有但是 的值一般是不知道的， 难以准确计算。称 为 在区间 的平均高度平均高度。寻找寻找 的一种的一种近似算法近似算法，便可得到一种数值积分公式。 badxxffI)()(xf)(f)(f,ba)(f )()()()(bafabdxx

6、ffIba)(xf)(xf)(xf第七章第七章 数值微积分数值微积分例如左矩形公式左矩形公式：右矩形公式右矩形公式：中矩形公式中矩形公式： )()()(afabdxxffIba )()()(bfabdxxffIba 2)()(bafabdxxffIba第七章第七章 数值微积分数值微积分 对于一般情况，由定积分定义， 在 内n+1个节点 处高度为 ，通过加权 后再求和，从而得到积分的近似值，这类公式一般形式为其中 称为求积节点求积节点； 为求积系数求积系数（或节点 的权权），它仅与节点值及区间 有关，而与被积函数 的形式无关。 记称 为积分公式的余项余项或是截断误差截断误差。)(xf,baix)

7、, 1 , 0( )(nixfiiA niiibaxfAdxxffI0)()(ixiAix,ba)(xfniiibanxfAdxxffR0)()( fRn第七章第七章 数值微积分数值微积分 数值求积方法的特点是直接利用积分区间 上一些离散节点处的被积函数值进行线性组合来近似计算定积分的值，从而将定积分的计算归结为函数值的计算函数值的计算，这就避开了牛顿-莱布尼兹公式需要求解原函数的困难，并为计算机计算积分提供了可行性。 数值求积公式的节点可以包含积分区间的端点，也可以不包含。求积节点包含积分区间端点时称为闭型求积公式闭型求积公式，如梯形公式；求积节点不包含积分区间的端点时，称为开型开型求积公式

8、求积公式，如中矩形公式。ba,1 牛顿-柯斯特（Newton-Cotes）公式 用拉格朗日插值多项式 作为 的近似函数，设 上的节点为则有其中则计算定积分时，有)(xLn)(xf,babxxxxan210niiinxfxlxL0)()()(nijjjijixxxxxl0)( niibaibaniiibanbaxfdxxldxxfxldxxLdxxffI00)()()()()()(1 牛顿-柯斯特（Newton-Cotes）公式 记则有公式上式称为插值型求积公式插值型求积公式。其中 只与插值节点 有关，与被积函数 无关。截断误差截断误差：式中 ， 。dxxlAbaii)( niiibaxfAdx

9、xffI0)()(iAix)(xfbanxnbanndxxfndxxLxffR)()()!1(1)()(1)1(niinxxx01)()(),(bax1 牛顿-柯斯特（Newton-Cotes）公式l牛顿牛顿- -柯特斯公式柯特斯公式:将 区间n等分，记节点 ，则对求积系数 的公式做变量替换 ，得到记则有,ba), 1 , 0(niihaxinabhiAthax nnijjinnnijjbanijjjijidtjtininabdtjijthdxxxxxA00000)()!( !) 1()()( nnijjinninidtjtininC00)(), 1 , 0(,)()!( !) 1(1)()(

10、niiCabA1 牛顿-柯斯特（Newton-Cotes）公式 于是上式称为牛顿牛顿- -柯特斯求积公式柯特斯求积公式。其中 是不依赖于 和区间 的常数，称为柯特斯系数柯特斯系数。柯特斯系数的性质：柯特斯系数的性质： (1)对称性： ； (2)权性： 。 niinibaxfCabdxxffI0)()()()()(niC)(xf,ba)()(ninniCC10)(niniC1 牛顿-柯斯特（Newton-Cotes）公式 根据区间的等分数n不同，计算柯特斯系数，得到不同的积分公式：1）n=1：梯形公式梯形公式，实际就是两点两点线性插值积分2）n=2：辛普森辛普森（Simpson）公式或抛物线抛物

11、线公式，实际就是三点三点抛物插值积分2110)1(1)1(0tdtCC niinibabfafabxfCabdxxffI0)()()(2)()()(61) 1(4120)2(2)2(0dtttCC64)2(4120)2(1dtttC niinibabfbafafabxfCabdxxffI0)()(24)(6)()()(1 牛顿-柯斯特（Newton-Cotes）公式3）n=4：柯特斯公式柯特斯公式，实际就是五点五点插值积分404)4(4)4(0907)4)(3)(2)(1(! 4! 04) 1(dtttttCC403)4(3)4(19032)4)(3)(2(! 3! 14) 1(dtttttC

12、C402)4(29012)4)(3)(1(! 2! 24) 1(dtttttC )(7)(32)(12)(32)(790 )()()(432100)(xfxfxfxfxfabxfCabdxxffIniiniba1 牛顿-柯斯特（Newton-Cotes）公式定理定理：若 在 上连续，则梯形公式的截断误差梯形公式的截断误差若 在 上连续，则辛普森公式的截断误差辛普森公式的截断误差若 在 上连续，则柯特斯公式的截断误差柯特斯公式的截断误差)(xf ,ba,),(12)(31bafabfR )()4(xf,ba,)(2901)(2880)()4(5)4(52bafabfabfR,)(44958)(1

13、013760)()6(7)6(74bafabfabfR)()6(xf,ba1 牛顿-柯斯特（Newton-Cotes）公式例例1 1 分别用梯形公式、辛普森公式和柯特斯公式计算积分 ，并与精确值进行比较。解解：（1）梯形公式（2）辛普森公式15 . 0dxxI4267767. 0) 15 . 0(25 . 0)()(2bfafabI43093403. 0) 175. 045 . 0(65 . 0)(24)(6bfbafafabI1 牛顿-柯斯特（Newton-Cotes）公式（3）柯特斯公式（4）精确值43096407. 0 )7875. 03275. 012625. 0325 . 07(90

14、5 . 0)(7)(32)(12)(32)(790bfefcfdfafabI43096441. 03215 . 03xI1 牛顿-柯斯特（Newton-Cotes）公式例例：用辛普森公式计算积分 的近似值，并求截断误差。解解： 又5 . 111sindxx0.36076 0.61837)0.71736441470.0833(0.8 5 . 11sin25. 11sin41sin65 . 01sin5 . 11dxxxxxxxxxf1sin3611cos1224)(6875)4(733611224max)(max68755 . 11)4(5 . 11xxxxxfxx1 牛顿-柯斯特（Newton

15、-Cotes）公式截断误差：000792. 07328805 . 0)(28805 . 02)4(22ffR1、分别用梯形公式和辛普森公式计算积分并估计误差（计算取5位小数）。 2、证明柯特斯系数具有性质 。1 牛顿-柯斯特（Newton-Cotes）公式10dxeIx 10niniC2 龙贝格（Romberg）求积公式l牛顿牛顿- -柯特斯公式问题柯特斯公式问题： 当积分区间比较大时，在整个区间上用牛顿-柯特斯公式求积分，精度难以保证。如果增加节点，使用高阶的牛顿-柯特斯公式，一方面系数计算复杂，舍入误差增大；另一方面，高次等距插值多项式 可能引起严重的振荡现象，使得求积公式数值不稳定。l复

16、化求积基本思想复化求积基本思想： 将区间 等分为n个子区间，在每个子区间 上用基本求积公式进行积分，然后再累加得出新的求积公式（即复化求积公式）。这样既可以提高结果的精度，算法实现也简便。)(xLn,ba), 2 , 1(,1nixxii2 龙贝格（Romberg）求积公式将区间 分为n等分，记分点 ，步长 ，子区间为 。l复化梯形公式复化梯形公式：在每个子区间 上应用梯形公式，积分值记为,ba), 1 , 0( niihaxinabh), 2 , 1( ,1nixxii,1iixxiI), 2 , 1()()(2)(111nixfxfxxdxxfIiiiixxiiinniiniiiniini

17、xxbaTbfxfafhxfxfhIdxxfdxxfIii )()(2)(2)()(2 )()(11111112 龙贝格（Romberg）求积公式l复化辛普森公式复化辛普森公式：l复化柯特斯公式复化柯特斯公式：)()(24)(6111021bfxfxfafhSniiniin221hxxii)(7)(14321232)(79011104310211041bfxfxfxfxfafhCniiniiniiniin43 ,2 ,4432141hxxhxxhxxiiiiii2 龙贝格（Romberg）求积公式l复化求积公式误差复化求积公式误差复化梯形公式复化梯形公式：复化辛普森公式复化辛普森公式：复化柯特

18、斯公式复化柯特斯公式：)(12)(1223fhabfnhfRT )(2880)(2880)4(4)4(5fhabfnhfRS)(1013760)(1013760)6(6)6(7fhabfnhfRC2 龙贝格（Romberg）求积公式例例1 计算积分 ，若要求误差不超过 ，分别用复化梯形公式和复化辛普森公式计算，问至少需各取多少个节点 ？解：解：（1）复化梯形公式的截断误差 ，依题意即为已知 ， ，所以有得 ，取 10dxeIx41021)(122fhabfRT 421021)(12 fhabfRT1, 0baxexf)(eexfxfx )()()4(42210211212)( nenffRT3

19、 .67n68n2 龙贝格（Romberg）求积公式（2）复化辛普森公式截断误差依题意有解得 ，取)(2880)(2880)4(4)4(5fhabfnhfRS444)4(102128802880)(nenffRS1 . 2n3n5、分别用复化梯形公式和复化辛普森公式计算积分（取11个节点，6位小数计算）。6、用积分 计算 ，要使计算误差不超过，问用复化梯形公式时至少取多少个节点？ 2 龙贝格（Romberg）求积公式1002dxeIx2ln2182dxx2ln510212 龙贝格（Romberg）求积公式l变步长求积公式变步长求积公式复化求积公式特点：复化求积公式特点： 事先给出合适步长； 误

20、差的估计需要已知被积函数，且有时较难估计上界。变步长求积公式特点：变步长求积公式特点： 步长逐次分半，反复利用复化求积公式进行计算，并可查看相邻两次计算结果的误差是否达到要求，直到满足计算精度要求。2 龙贝格（Romberg）求积公式 设求积区间 分为n等分，复化梯形公式为将每个子区间分半，即取区间的中点，则求积区间 分为2n等份，复化梯形公式,ba)()(2)(211bfxfafhTniin,ba10211021111021112)(221)(2)()(2)(221)()(2)(2)(4niinniiniiniiniinxfhTxfhbfxfafhbfxfxfafhT2 龙贝格（Romber

21、g）求积公式 即特点特点： 当计算 时， 为已知数据，只需要计算新增的分点 的函数值 ； 计算过程常用 是否满足作为控制计算精度的条件。10212)(221niinnxfhTTnT221ix)(21ixfnTnnTT22 龙贝格（Romberg）求积公式l龙贝格求积公式龙贝格求积公式基本思想基本思想：用复化求积公式的截断误差对递推化求积公式进行修正，希望提高该公式的收敛速度。 由复化梯形公式截断误差可得表示二分一次，误差变为原来的1/4倍。由上式得到41)()(41212 ffTITInnnnnSTTI313422 龙贝格（Romberg）求积公式同样可以得到从而由得到 ，即为龙贝格公式。16

22、12nnSISInnnCSSI151151626412nnCICInnnRCCI631636422 龙贝格（Romberg）求积公式因此在步长二分的过程中，分别应用公式进行修正，逐步将粗糙的梯形公式积分值 加工成精度较高的龙贝格公式积分值 ，这种方法称为龙贝格算法龙贝格算法。nnnTTS31342nnnSSC15115162nnnCCR63163642nTnR2 龙贝格（Romberg）求积公式龙贝格计算流程图：2 龙贝格（Romberg）求积公式例例2 用龙贝格算法计算积分 ，要求误差不超过 。解：解：计算过程中，1）先计算当m=2时， ，10214dxxI51021214)(, 1, 0

23、xxfbanT3)1 ()0(21)()(21ffbfafabT11h1 . 32 . 3215 . 1)5 . 0(221112fhTT2 龙贝格（Romberg）求积公式当m=3时，当m=4时， 5 . 02h25. 03h13117647. 3)75. 0()25. 0(221224ffhTT13898849. 3 )875. 0()625. 0()375. 0()125. 0(221348ffffhTT2 龙贝格（Romberg）求积公式2）应用公式计算 ，最后得到nnCS ,nR133333. 33311 . 3343134121TTS141569. 31 . 331131177.

24、3343134242TTS142118. 3 133333. 3151141569. 315161511516121SSC141595. 3141569. 3151141593. 315161511516242SSC2 龙贝格（Romberg）求积公式计算结果如下表所示kkn2kT212kS22kC32kR013123.13.133333243.1311773.1415693.142118383.1389893.1415933.1415953.1415874163.1409423.1415933.1415933.1415935323.1414303.1415933.1415933.141593

25、141587. 3142118. 3631141595. 363646316364121CCR7、用龙贝格公式计算积分要求误差不超过 。（0.71327）2 龙贝格（Romberg）求积公式102dxeIx5104 数值微分 实际问题中变量之间的关系通常是由离散的数据 给出，或者函数的表达式过于复杂,如果求解变量的变化率，只能用数值方法求它的近似导数数值微分。即：给定函数表 ，求函数 在节点 处的导数值。 )10( )(,(,n, ,ixfxii)10( )(,(,n, ,ixfxii)(xfix4 数值微分 在微积分中，导数表示函数在某点上的瞬时变化率瞬时变化率，它是平均变化率的极限；在几何

26、上可解释为曲线的斜率斜率，物理上可解释为物体的变化速率变化速率。 有如下三种导数 的定义形式最简单的计算数值微积分的方法就是用函数的差商近似函数的导数，即取极限的近似值。 xf hhxfhxfhhxfxfhxfhxfxfhhh2limlimlim0004 数值微分 用向前平均变化率(差商)近似导数有式中 的位置在 的前面，因此称为向前差商向前差商。 根据泰勒展开式向前差商的截断误差截断误差为hxfhxfxfkkk hxxfhxf hxfhxfkkkkk , ,! 22hxkkx hOfhhxfhxfxfxRkkk 24 数值微分 同样得到向后差商向后差商及其截断误差截断误差hhxfxfxfkk

27、k hOhhxfxfxfxRkkk4 数值微分 中心差商为：截断误差截断误差为hhxfhxfxfkkk2 132! 3! 2fhxfhxf hxfhxfkkkk 232! 3! 2fhxfhxf hxfhxfkkkk ,612 222212hOfhffhhhxfhxfxfxRkkk hxhx00,4 数值微分 差商的几何意义：几何意义：由微积分中的极限定义表达式表示 在 处切线的斜率，即图中直线P的斜率。 而差商表示过点 和 的弦线Q的斜率。 因此数值微分是用近似值内接弦的斜率代替准确值切线的斜率。hxfhxfxfh0000lim xf0 xx hxfhxf0000,xfxhxfhx00,4

28、数值微分4 数值微分 类似可以得到 上述差商型求导公式简洁，可利用节点值快速计算节点的导数，但计算精度较低，在实际计算中，步长值h应取一个合适合适的值。)()()(2)()(2211hOhxfxfxfxfkkkk 4 数值微分例例：用中心差商公式计算 在 处的一阶导数的近似值，取 分别计算，并比较结果的精度。解：准确值为对于 ，其中心差商公式为xxf)(2x00005. 0 ,0005. 0 ,005. 0 ,05. 0 , 5 . 0h3535534. 0221)2( fxxf)(hhhf222)2(4 数值微分计算结果如下表 由上例看出，步长过大，则截断误差较大；但步长太小又会导致舍入误差

29、的增长，在实际计算时，希望在保证截断误差满足精度要求的前提下选取尽可能大的步长 。hh2h2hhh222hhhf222)2(0.51.5811391.2247450.356394-0.0028410.051.4317821.3964240.353581-0.0000280.0051.4159801.4124450.353554-0.00000030.00051.414391.414040.3500000.000051.414231.414200.3000004 数值微分l插值型求导公式插值型求导公式 构造数值微分公式的普遍方法是用一个易于计算其微分的函数近似代替原问题中的函数 ，插值多项式是最

30、简单的函数。构造n次拉格朗日插值多项式 ，并取 的值作为 的近似值，建立的数值公式 称为插值插值型求导公式型求导公式。截断误差)(xf)(xLn)(xLn)(xf )()(xLxfn)()!1()()()(1)1(xnfxLxfnnn4 数值微分其中 且依赖于 （ 是 的函数）。上式求导等式右侧第二项无法求得，为此，限定求节点处的导数值，于是),(0nxxx)()!1()()()!1()()()()1(11)1(nnnnnfdxdnxxnfxLxf), 1 , 0(nixi), 1 , 0( ),()!1()()()(1)1(nixnfxLxfinninix4 数值微分当节点等距分布时，常用的

31、数值微分公式有：（1）一阶两点公公式10101011)(iiyxxxxyxxxxxL)(1)(11iiyyhxL)(1)(2)(1)(111iiiiiyyhfhyyhxf )(1)(2)(1)(121iiiiiyyhfhyyhxf ) 1, 1 , 0(),(11nixxii), 1(),(12nixxii4 数值微分（2）一阶三点公式（3）二阶三点公式 ,3)43(21)(13221fhyyyhxfiiii)2, 1 , 0(ni ,6)(21)(23211fhyyhxfiii) 1, 2 , 1(ni ,3)34(21)(33212fhyyyhxfiiii), 3 , 2(ni ,12)2(1)(42112fhyyyhxfiiii ) 1, 2 , 1(ni4 数值微分例例：根据下列 的数值表用中点公式计算 的近似值，并估计误差，同时把结果和精确值比较。