微积分的数值计算方法.docx

微积分的数值计算方法 第七章 微积分的数值计算方法 7.1 微积分计算存在的问题/数值积分的基本概念 1. 微分计算问题 求函数的导数(微分)，原则上没有问题。当然，这是指所求函数为连续形式且导数存在的情形。但如果函数一表格形式给出，要求函数在某点的导数值；或者是希望某点的导数值只用其附近离散点上的函数值近似地表示，这就是新问题了，它称为微分的数值计算，或称为数值微分。 2定积分计算问题 计算函数f在a,b上的定积分 I=òf(x)dx ab当被积函数f的原函数能用有限形式F(x)给出时，可用积分基本公式来计算： I=òaf(x)dx=F(b)-F(a) 然而，问题在于： f的原函数或者很难找到，或者根本不存在；f可能给出一个函数表；仅仅知道f是某个无穷级数的和或某个微分方程的解等等。这就迫使人们不得不寻求定积分的近似计算，也称数值积分。 3数值积分的基本形式 数值积分的基本做法是构造形式如下的近似公式 bòf(x)dx»åAf(x) (7.1.1) akkk=0bn或记成 nòf(x)dx=åAf(x)+Rf (7.1.2) akknk=0bnI=åAkf(xk) 和 Rnf 分别成为a,b上的f的数值求积公式及其k=0*余项(截断误差)，xk和Ak(k=0,1,L,n)分别称为求积节点和求积系数(求积系数与被积函数无关)。 这种求积公式的特点是把求积过(极限过程)程转化为乘法与加法的代数运算。构造这种求积公式需要做的工作是：确定节点xk及系数Ak(k=0,1,L,n)，估计余项Rnf以及讨论I的算法设计及其数值稳定性。 4插值型求积公式 如何构造求积公式呢？基本的技术是用被积函数f的Lagrange插值多项式*Ln(x)近似代替f，也即对a,b上指定的n+1个节点1 a£x0<x1<¼¼<xn£b及相应的函数值f(x0),f(x1),L,f(xn)，作 nf(x)=Ln(x)+Rn(x)=ålk(x)f(xk)+k=01(n+1)!f(n+1)(x)w(n+1)(x) 代入(7.1.2)式等号左边有 òf(x)dx=òL(x)dx+òR(x)dx aananbbbn =åòalk(x)dxf(xk)+k=0b1(n+1)!òafb(n+1)(x(x)w(n+1)(x)dx 或写成形如(7.1.2)式的一般形式： 其中 Ak=bòf(x)dx=åAf(x)+Rf (7.1.4) akknk=0kbnòl(x)dx (k=0,1,L,n) (7.1.5) a Rnf=1(n+1)!òafb(n+1)(x(x)w(n+1)(x)dx (7.1.6) 称(7.1.4)为插值型求积公式。下面将要介绍的几种实用求积公式无不都是插值型公式的某种具体形式。 由上述定义，用余项公式可以衡量数值求积公式的精确程度。不过，由余项公式(7.1.6)可见，如果f为次数£n的多项式，则Rnf中有fbn(n+1)=0,故Rnf=0，从而公式(7.1.4)成为等式 òaf(x)dx=åAkf(xk)，这就是说，k=0当被积函数f为次数£n的多项式时，其相应的插值型求积公式不是近似公式，而是准确公式。据此，人们引入了另一个衡量数值求积公式近似程度好坏的“代数精度” 概念。 5. 代数精度 定义7.1.1 如果数值求积公式òaf(x)dx»mbnåAf(x)，当f(x)是0m次kkk=0bn多项式(可表示为f(x)=1,x,L,x)时，均有òaf(x)dx=åAf(x)，也kkk=0 2 n即有Rnf=0，便称求积公式åAkf(xk)至少具有m次代数精度；如果当k=0m+1)时，只能是f(x)是m+1次多项式(可表示为f(x)=x也即有Rf¹0 ，便称求积公式åAf(x)具òf(x)dx»åAf(x)，akkbnnnkkk=0k=0有 m次代数精度。 定理 7.1.1 有n+1个节点的插值型求积公式(7.1.4)，至少具有n次代数精度。 例7.1.1 确定节点x1,x2,x3和系数A，使得下列形式的求积公式 òf(x)dx=Af(x)+f(x)+f(x) -11231具有3次代数精度。 解 利用代数精度定义 取 f(x)=1,x,x,x令 23ìò-11dx=A(1+1+1)ï1ïò-1xdx=A(x1+x2+x3) í12222ïò-1xdx=A(x1+x2+x3)ï13333xdx=A(x+x+x)123îò-1计算得 A=123,x1=-1223,x2=0,x3=12即得至少具有3次代数精度的求积公式为 òf(x)dx»-11f(-412)+f(0)+f(12) 验证上式两端对f(x)=x，有 左边=2/5 ，右边=1/3 两者不相等，故可知所得公式正好为3次代数精度。 7.2 Newton-Cotes型求积公式 构造具体形式的插值型求积公式(6.1.4)的第一个想法是考虑求积节点为等距节点的情况(这时公式可能比较简单，而且节点也就定下来了)，这就是Newton-Cotes(牛顿-柯特斯公式)，而Newton-Cotes型公式当节点只有两点和三点时，就是熟知的梯形公式和Simpson(辛普生)公式。 1. Newton-Cotes 求积公式 3 考虑积分区间a,b上的节点为等距节点 xk=a+kh(h=b-anb,k=0,1,L,n) (7.2.1 ) 由(7.1.5)式可得插值型求积系数 Ak=òbnalk(x)dx=òÕaj=0j¹kx-xjxk-xjdx (通过变换x=a+th) =(b-a)ònk!(n-k)!n0(-1)n-kn0t(t-1)L(t-k+1)(t-k-1)L(t-n)dt记C(n)k=ònk!(n-k)!(-1)n-kt(t-1)L(t-k+1)(t-k-1)L(t-n)dt (7.2.2) *并称它为Cotes系数，从而得等距节点的插值型的n阶Newton-Cotes求积公式INC，òbnaf(x)dx=(b-a)åCkk=0(n)f(xk)记为INC*。(7.2.3) (n)显然，对不同的n以及k=0,1,L,n,按公式(7.2.2)算出Ck阶的Newton-Cotes公式。 ，按公式可构造出不同下面是Cotes系数表 n Ck k=0,1,L,n 1 1/2 1/2 2 1/6 4/6 1/6 3 1/8 3/8 3/8 1/8 4 7/90 16/45 2/15 16/45 7/90 8 989/28350 5888/28350 -928/28350 10496/28350 -4540/28350 989/28350 2. 梯形公式/Simpson公式 当n=1，即两个节点，取积分区间a,b两个端点a,b为节点时，由上表(7.2.1)可构造求积公式 称为梯形公式。 (n)òbaf(x)dx»b-a2f(a)+f(b)记为T (7.2.4) 4 当n=2，即3个等距节点，取积分区间a,b两端点a,b及区间中点点时，由表7.2.1可构造求积公式 a +b 2为节òbaf(x)dx»(b-a6)f(a)+4f(a+b2)f(b)记为S (7.2.5) 称为Simpson(辛普生)公式或抛物线求积公式。 类似地，当n=3可构造出所谓Simpson3/8公式；当n=4可构造出Milne(米尔尼)公式(有时也称Cotes公式)等等。 3. 误差分析 关于Simpson-Cotes型求积公式的截断误差或称余项 bRf=òaf(x)dx-INC *的分析，计算数学家为我们证明了一个一般性定理。 定理 7.2.1 Newton-Cotes公式(7.2.3)的余项Rnf可表示为： (1) 对n为奇数的情形，设fÎC(n+1)a,b，则 (h) hÎ(a,b) Rnf=rnhn+2f(n+1)rn=1(n+1)!òm(m-1)L(m-n)dm 0n(2) n为偶数的情形，设fÎC(n+2)a,b，则 (h) hÎ(a,b) Rnf=rnhn+3fn(n+2)rn=1(n+2)!ò0m(m-1)L(m-n)dm 2由定理可知，Newton-Cotes公式作为插值型公式的特例，当n为奇数时，保持至少n次代数精度不变，只在n为偶数时，代数精度才略加1次。 梯形公式具有1次代数精度，余项为 RTf=-(b-a)123f(h)hÎ(a,b) (7.2.6) ''Simpson公式具有3次代数精度，余项为 Rsf=-190(b-a2)f5(4)(h)hÎ(a,b) (7.2.7) 当n为较大数时，Newton-Cotes公式的数值稳定性还可能有问题，正因为如此，对n>2的公式我们不感兴趣。 5 例 7.2.1 用不同的方法计算并比较下列积分： 用传统的方法 ò0edx=e|0=e-1=1.71828(e=2.71828L) 1xx1用梯形公式 ò0edx»1x12 e+e=1.859140101|RTf|=|-用Simpson公式 112(1-0)e£1013h112e=0.2265235ò0edx»111x165e+4e2+e=1.7188612 h |RSf|=|- e£e=0.00094385902288014. 数值的稳定性 若取f=1，可推出(n³1) òa1dx=(b-a)åCk´1Þb-a=(b-a)åCk k=0k=0nbn(n)n(n)ÞåCkk=0(n)=1 假设用n阶的Newton-Cotes公式做实际计算，而且f(xk)可能使用近似值f(xk)，这反映到计算上就有误差 nnn(b-a)åCkf(xk)-(b-a)åCkk=0k=0(n)(n)f(xk)=(b-a)åCkf(xk)-f(xk)k=0(n)记 e=max|f(xk)-f(xk)| 则有 0£k£nn|åCkf(xk)-åCkk=0k=0(n)n(n)nf(xk)|£å|Ckk=0(n)|f(xk)-f(xk)£eå|Ckk=0n(n)| 由此可见 (1) 若 Ck>0 (k=0,1,L,n),则 nnnå|Ck=0(n)kn|=åCkk=0(n)=1，于是有 |åCkf(xk)-åCkk=0k=0(n)n(n)f(xk)|£e 即计算公式是稳定的。 6 (2) 若存在Ck<0，则有 nnå|Ck|>åCkk=0k=0(n)n(n)=1 这样的话，初始数据误差有可能引起计算结果误差的增大，即计算的数值稳定性得不到保证。 为了方便以后应用，把上面的论证归结为更一般的结果： 定理 7.2.2 插值型求积nåAkf(xk)中的求积系数Ak之和为k=0nåAk=b-a(7.2.11) k=0n定理 7.2.3 若插值型求积åAkf(xk)中的求积系数Ak>0(k=0,1,L,n) k=0则该求积公式公式是稳定的。 5. 复化梯形公式/复化Simpson公式 现在，在每个子区间上利用梯形公式，有 -1òbnxn-1k+1af(x)dx=åf(x)dx»åh(xk)+f(xk+1) k=0òxkk=02fn-1 =h2f(x0)+2åf(xk)+f(xn)记为Tn (7.2.12) k=1T2n称为复化梯形公式。其余RTn由(7.2.6)式当fÎCa,b，有 n-1Rh3Tn=-12åf''(hk) hkÎxk,xk+1 k=02n =-(b-a)hf''(hk)12åk=0n=-(b-a)h2 12f''(h) hÎa,b (7.2.13) 如果在每个区间上利用Simpson公式，有 n-1òbn-1f(x)dx=xk+1haåk=0òxf(x)dx»(x(xk+xk+1kåk)+4fk=06f2)f(xk+1) -1n-1 =hn6f(x0)+4åf(x)+2åf(xk)+f(xn)记为Sn k=0k+12k=1 7 其中 xk+12=xk+h2，其余RS由(7.2.7)式当fÎCa,b，类似地可得 n4RSn(b-a)h4=-f1802(b-a)h28804(4)(h) hÎa,b (7.2.15a) 或 =-f(4)(h) hÎa,b (7.2.15b) 6. 复化公式的收敛性 复化求积方法带来新的收敛性问题，即随着子区间越分越细，也即理论上n®¥时，复化公式的值是否越来越接近积分的准确值？ 事实上，只需 12n-1fÎCa,b,由Tn公式(7.2.12)显然有 nTn=åf(xk)h+k=0åk=1f(xk)h¾n¾¾®®¥12òf(x)dx+abòbaf(x)dx=òbaf(x)dx即复化梯形公式Tn的收敛性是有保证的。类似地，复化Simpson公式的收敛性也是有保证的。 定义 7.2.1 如果一种复化求积公式In，当n*®¥(h®0)时有渐进关系式 *limn®¥I-Inhp*=C(C¹0) 则称求积公式In是 对复化梯形公式 limI-Inh2*p阶收敛的。 h®0=lim(-h®0å121ba1n-1f''(hk)h) k=0=-12òf''(x)dx=-112f'(b)-f'(a) 可知复化梯形公式有2阶收敛性，且在h很小时有误差估计式 I-Tn»-112hf'(a)-f'(b) 2对于复化Simpson公式有4阶收敛性，且在h很小时有误差估计式 I-Sn»-h4f'''(a)-f'''(b) 18021 8 例 7.2.2 求 f(x)=sinxx在0,1上的积分 I=ò1sinxx0dx，要求 -3(1) 利用复化梯形公式计算上述积分，要求截断误差不超过(1/2)´10 (2) 利用(1)中复化梯形公式所用的节点，改用复化 Simpson公式做积分积分计算，并估计截断误差。 解 (1) 要使用复化梯形公式计算上述积分的截断误差不超过(1/2)´10-3，即要求 |Rb-a22TNf|=|-12hf''(h)|£112h|f''(h)|£12´10-3对此，先估计|f''(h)|。注意到 f(x)=sinx1x=ò0cos(xt)dt 因而有 f(k)(x)=ò1dk1kkp0dxkcos(xt)dt=ò0tcos(xt+2)dt 于是可得 |f(k)(x)|=ò1kkp|dt£ò1k10t|cos(xt+2)0tdt£k+1利用此结果得，可知要求 |R11TNf£|1212hf''(h)|£12h22+1£12´10-3只需h£0.1342 或 n³b-a0.1342=10.1342=7.4516 可知只需取8个等分节点的复化梯形公式计算上述积分： I»T8=0.9456911 (2)利用同样的点做复化Simpson公式计算，这时可取h=14，则 I»S4=0.9460832 |R14(4)1SNf£|2880hf(h)|£142880(4)14+1£0.271´10-6 9