第7章 数值积分与数值微分
第七章数值积分与数值微分
积分问题最早来自于几何形体的面积、体积计算，也是经典力学中的重要问题（例如计算物体的重心位置）. 在现实应用中，很多积分的结果并不能写成解析表达式，因此需要通过数值方法来计算. 数值微分是利用一些离散点上的函数值近似计算某一点处的函数导数，它针对表达式未知的函数. 本章介绍一元函数积分（一重积分）和微分的各种数值算法，它们也是数值求解积分方程、微分方程的基础.
7.1数值积分概论
7.1.1基本思想
考虑如下定积分的计算：
I(f)≡∫f(x)dx
b
a
,(7.1) 其中函数f: ?→?，首先应想到的是微积分中学习过的牛顿-莱布尼兹（Newton-Leibniz）公式：
∫f(x)dx
b
a
=F(b)?F(a) ,
其中F′(x)=f(x),即F(x)为f(x)的原函数. 但是，诸如e x2,sinx
x
,sinx2等表达式很简单的函数却找不到用初等函数表示的原函数，因此必须研究数值方法来近似计算积分. 另一方面，某些函数的原函数虽然可以解析表示，但其推导、计算非常复杂，此时也需要使用数值积分方法.
一般考虑连续的、或在区间[a,b]上可积①的函数f(x)，则根据积分的定义有：
lim n→∞,?→0∑(x i+1?x i)f(ξi)
n
i=0
=I(f) , (7.2)
其中a=x0 上述讨论表明，近似计算积分I(f)的数值积分方法（numerical quadrature）一般具有如下形式： I n(f)≡∑A k f(x k) n k=0 , (7.3) 其中a≤x0 推导求积公式的一种方法是用多项式函数p(x)来近似f(x)，则可期望有以下的近似关系：①连续函数在闭区间内一定有界、可积，而可积函数则可能不连续. ∫f (x )dx b a ≈∫p (x )dx b a , 其中多项式函数的积分很容易通过牛顿-莱布尼兹公式求出. 假设使用拉格朗日插值法构造p (x )，区间[a, b]内的插值节点为x 0,x 1,…,x n ，则 p (x )=L n (x )=∑f (x k )l k (x )n k=0 , l k (x)为拉格朗日插值基函数. 由此得到求积公式为： I n (f )=∫∑f (x k )l k (x )n k=0dx b a =∑f (x k )∫l k (x )dx b a n k=0 . (7.4) 由于l k (x )为拉格朗日插值基函数，一旦插值节点确定，积分∫l k (x )dx b a 可方便地计算出来. 这 种用多项式插值近似被积函数得到的求积公式(7.4)被称为插值型求积公式（interpolatory quadrature ），易知它也是一种机械求积公式，其积分节点就是插值节点，而积分系数 A k =∫l k (x )dx b a ,(k =0,1,?,n ). (7.5) 下面的例子推导n=0, n=1两种情况下的插值型求积公式，它们分别称为中矩形公式（midpoint rule ）和梯形公式（trapezoid rule ）. 例7.1 (中矩形公式与梯形公式)：根据n=0, n=1两种情况对应的拉格朗日插值推导相应的求积公式，假设插值节点分别为区间[a, b]的中点和两个端点. [解] 当n=0时，按题意设x 0=(a +b)/2, 由于0次拉格朗日插值多项式为常数，则 L 0(x )=f (x 0) 因此， I 0(f )=∫f (x 0)dx b a =(b ?a )f ( a + b 2) . (7.6) 当n=1时，按题意设x 0=a , x 1=b , 利用线性拉格朗日插值基函数和公式(7.5)，求出： A 0=∫l 0(x )dx b a =∫x ?b dx b a =b ?a , A 1=∫l 1(x )dx b a =∫ x ?a b ?a dx b a =b ?a 2 . 因此， I 1(f )=∑A k f (x k )1 i=0= b ?a 2[f (a )+f(b)] . (7.7) 中矩形公式(7.6)和梯形公式(7.7)具有很直观的几何意义，即分别用矩形面积和梯形面积来近似函数曲线和横轴围成区域的面积（如图7-1）. 7.1.2求积公式的积分余项与代数精度 定义7.1：对于计算积分I(f)的求积公式I n (f)，称I(f)?I n (f)为该公式的积分余项，常记为R [f ]. 积分余项反映了求积公式的截断误差，是衡量求积公式准确度的重要依据. 假设I n (f)为某个插值函数p(x)的积分，则 R [f ]=∫[f (x )?p (x )]dx b a , 即积分余项等于插值余项的积分. 对于插值型求积公式(7.4)，有 R [f ]=∫[f (x )?L n (x )]dx b a =∫f (n+1)(ξ)(n +1)!ωn+1(x )dx b a , (7.8) 其中ξ依赖于x . 下面介绍代数精度的概念，它是衡量求积公式准确度的另一个重要标准. 定义7.2：如果某求积公式对于次数不超过m 的多项式均准确成立，但对于m+1次多项式可能不准确，则称该求积公式具有m 次代数精度（degree of exactness ）. 上述定义表明，一个求积公式具有较高次的代数精度，就意味着它能准确计算次数较高的多项式的积分②. 应注意，在某些情况下代数精度并不是越高越好. 要判断一个机械求积公式的代数精度，最直接的方法是考察当被积函数分别为1,x,…,x m 时求积公式的准确性. 下面给出一个定理，其证明留给感兴趣的读者思考. 定理7.1：机械求积公式I n (f )=∑A k f (x k )n k=0至少有m 次代数精度的充要条件是当f (x )分别为1,x,…,x m 时， I (f )=I n (f ) . 我们讨论的所有求积公式都至少具有0次代数精度，因此根据定理7.1它们应对f (x )=1的积分准确，则推出： ∑A k n k=0=∫1dx b a = b ?a , (7.9) 这说明积分系数之和等于区间长度. 根据插值型求积公式的含义，容易得出如下定理，其证明留给读者思考. 定理7.2：机械求积公式I n (f )=∑A k f (x k )n k=0是插值型求积公式(7.4)的充要条件是它至少有n 次代数精度. 考察例7.1中的中矩阵公式和梯形公式，可得出它们都具有1次代数精度. 一般地，在部分积分节点、积分系数已知的情况下，利用定理7.1可建立方程求解剩余的积分系数或节点，使其达到一定的代数精度. 而且定理7.1中求积公式的形式还可以更一般，只要是函数值或其导数的线性组合即可，下面的例子说明了这种情况. 例7.2(求积公式的代数精度)：用形如H 2(f )=A 0f (0)+A 1f (1)+B 0f ′(0)的求积公式近似积分I (f )=∫f(x)dx 10，试确定系数A 0, A 1, B 0，使公式具有尽可能高的代数精确度. [解] 根据题意可令f (x )=1，x ，x 2分别代入求积公式使H 2(f )=I(f)精确成立. ② 当然，这没有什么实际意义，因为很容易得到多项式函数的原函数. (a) 中矩形公式 (b) 梯形公式 图7-1 中矩形公式(a)和梯形公式(b)的示意图. 当f (x )=1时，得 A 0+A 1=∫1?dx 1 =1 当f (x )=x 时，得 A 1+ B 0=∫xdx 1 =12 当f (x )=x 2时，得 A 1=∫x 2dx =1310 联立上述三个方程，解得A 1=13,A 0=23,B 0=16 . 当f (x )=x 3时，容易验证上述求积公式不准确，因此H 2(f )最多具有2次代数精度. 7.1.3求积公式的收敛性与稳定性 实际使用的求积公式都是机械求积公式(7.3)，下面针对它给出求积公式的收敛性和稳定性的概念. 定义7.3：对于n 的值可以任意的一系列机械求积公式I n (f)=Σk=0n A k f (x k )，a ≤x 0 ? lim n→∞,?→0I n (f)=∫f (x )dx b a , 其中?=max 1≤k≤n (x k ?x k?1)，则称这一系列求积公式具有收敛性. 收敛性说明求积公式在积分节点逐渐增多、且节点间距逐渐变小时，其结果收敛到准确的积分值. 这个概念不同于公式(7.2)，后者反映的是被积函数具有可积性. 在实际应用中，求积公式具有收敛性非常重要，后面还将针对具体的公式加以讨论. 在讨论求积公式的稳定性之前，先分析数值积分问题的敏感性和条件数. 假设f(x)为准 确的被积函数， f ?(x )为实际计算时受扰动影响的被积函数，扰动的大小为δ=‖f (x )?f ?(x )‖∞ =max a≤x≤b |f (x )?f ?(x )|，则扰动对积分计算的影响为： |∫f (x )dx b a ?∫f ?(x )dx b a |≤∫|f (x )?f ?(x )|dx b a ≤(b ?a )δ, (7.10) 这说明，积分计算结果的误差最多为扰动的(b ?a)倍，积分区间的长度(b ?a)是绝对条件数的上限. 一般来说，数值积分问题是不太敏感的. 这一点不难理解，因为积分运算本身就是一个平均的过程，它不易受被积函数的某些微小变化的影响. 求积公式的稳定性反映计算过程中的扰动是否被放大、以及放大的程度. 具体来说，在计算机械求积公式时，需考虑积分节点的函数值出现误差时，它对结果产生的影响. 假设节 点函数值由f (x k )变为f ?(x k )，则数值积分的结果由I n (f )变为I n (f ?)，两者之差满足： |I n (f )?I n (f ?)|=|∑A k [f (x k )?f ?(x k )]n k=0 |≤∑|A k |n k=0?|f (x k )?f ?(x k )|≤(∑|A k |n k=0 )ε, (7.11) 其中ε=max 0≤k≤n |f (x k )?f ?(x k )|≤δ. 根据(7.9)式，若同时有A k >0,k =0,…,n ，则不等式 (7.11)变为： |I n (f )?I n (f ?)|≤(b ?a )ε≤(b ?a )δ . (7.12) (7.12)式表明求积公式的结果受扰动影响的程度与积分问题敏感性的结果(7.10)式一致，这是控制数值计算误差能达到的最佳情况. 将它作为一个标准，可定义求积公式的稳定性. 定义7.4：若对k =0,1,…,n ，均有A k >0, 则机械求积公式I n (f)=Σk=0n A k f (x k )是稳定的. 利用定义7.4很容易直接判断求积公式的稳定性. 在实际情况中，不稳定的求积公式的 积分系数绝对值之和Σk=0 n |A k |可能远大于b ?a ，从而导致函数值的扰动在计算结果上被放大很多. 本节介绍了求积公式的基本形式，以及积分余项、代数精度、收敛性和稳定性的概念，其中收敛性是针对一系列公式（积分节点数目逐渐增多）而言的. 后面介绍具体公式时，将考察单个公式的积分余项、代数精度和稳定性，并讨论积分节点数目逐渐增多时的收敛性. 此外，还应注意具体公式中计算函数值的次数，它是度量计算量大小的标准. 7.2牛顿-柯特斯公式 在积分区间上构造等距节点的多项式插值，对应的插值型求积公式为牛顿-柯特斯公式. 7.2.1 柯特斯系数与几个低阶公式 假设将积分区间n 等分，步长?=(b ?a)/n ，插值节点为x k =a +k?,(k =0,…,n)，则可得到等距节点的拉格朗日插值多项式. 根据公式(7.4), (7.5)的推导，得到型如I n (f )=∑A k f (x k )n k=0的求积公式，其中 A k =∫l k (x )dx b a =∫(x ?x 0)?(x ?x k?1)(x ?x k+1)?(x ?x n )(x k ?x 0)?(x k ?x k?1)(x k ?x k+1 )?(x k ?x n )dx b a . 这就是n 阶牛顿-柯特斯（Newton-Cotes ）公式. 引入变量代换：x =a +t?,t ∈[0,n], 则 A k =∫∏(t ?j )n j=0,j≠k ?dt n 0=∫∏(t ?j )n j=0,j≠k b ?a dt n 0 . (7.13) 令 C k (n )=1n ∫∏(t ?j k ?j )n j=0,j≠k dt n 0 ,(k =0,?,n ), (7.14) 它仅与阶数n 有关，而与区间大小无关，则积分系数 A k =(b ?a )C k (n ). (7.15) 公式(7.14)中的C k (n )常被称为柯特斯系数，可以预先计算出不同n 值对应的柯特斯系数，制 成一个表（如表7-1所示），根据它方便地写出各阶牛顿-柯特斯公式. 表7-1 柯特斯系数表 [本章知识点]: 机械求积公式；积分余项；代数精度；插值型求积公式及其代数精度；中矩形公式、梯形公式；求积公式的稳定性；求积公式的收敛性；牛顿-柯特斯公式及其代数精度；柯特斯系数；辛普森公式、柯特斯公式；梯形公式、辛普森公式的积分余项；复合梯形公式；等距节点求积公式的准确度阶数；复合辛普森公式；步长折半情况下的复合求积公式计算；复合梯形公式的余项展开式；理查森外推法；龙贝格积分算法；自适应积分算法的原理；高斯求积公式；带权积分的高斯求积公式；高斯点的性质以及高斯求积公式的构造；高斯-勒让德公式；向前、向后、中心差分公式；二阶中心差分公式；有限差分公式的准确度阶数；插值型求导公式；数值微分的外推算法. 算法背后的历史：“数学王子”高斯 卡尔·弗里德里希·高斯（Johann Carl Friedrich Gauss，1777 年4月30日—1855年2月23日）是德国著名数学家、物理学家、 天文学家、大地测量学家. 高斯的成就遍及数学的各个领域，在数 论、非欧几何、微分几何、超几何级数、复变函数论以及椭圆函数 论等方面均有开创性贡献. 高斯被公认为是十九世纪最伟大的数 学家，与阿基米德、牛顿齐名. 高斯于1777年4月30日生于布伦瑞克的一个工匠家庭. 1792 年，15岁的高斯进入布伦瑞克学院，在那里独立发现了二项式定 理的一般形式、数论上的“二次互反律”、“质数分布定理”和“算术几 何平均”. 1795年高斯进入哥廷根大学. 1796年，19岁的高斯取得了一个数学史上极重要的结果，就是论文“正十七边形尺规作图之理论与方法”. 1801年，高斯又证明了"Fermat素数"边数的正多边形可以由尺规法作出. 高斯从1807年起担任哥廷根大学教授兼哥廷根天文台台长直至逝世. 高斯的主要贡献 ●高斯消去法 ●质数分布定理和最小二乘法 ●高斯分布（标准正态分布） ●高斯求积公式 ●无穷级数的高斯判敛法 ●发现了谷神星的运行轨迹 ●证明代数基本定理 ●微分几何最重要的创始人 高斯的故事 一天，德国哥廷根大学，一个19岁的青年吃完晚饭，开始做导师单独布置给他的数学题. 正常情况下，他总是在两个小时内完成这项作业. 像往常一样，前两道题目在两个小时内顺利地完成了. 第三道题写在一张小纸条上，是要求只用圆规和一把没有刻度的直尺做出正17边形，他没多想，埋头做起来. 然而，做着做着，他感到越来越吃力. 困难激起了他的斗志：我一定要把它做出来！天亮时，他终于做出了这道难题. 导师看了他的作业后惊呆了. 他用颤抖的声音对青年说：“这真是你自己做出来的？你知不知道，你解开了一道有两千多年历史的数学悬案？阿基米、牛顿都没有解出来，你竟然一个晚上就解出来了！……我最近正在研究这道难题，昨天不小心把写有这个题目的小纸条夹在了给你的题目里. ”多年以 后，这个青年回忆起这一幕时，总是说：“如果有人告诉我，这是一道有两千多年历史的数学难题，我不可能在一个晚上解决它. ”这个青年就是数学王子高斯. 高斯在1810年提出了消去法，他的这个方法是为了简化二次型而不是为了矩阵分解. 实际上，求解线性方程组中的高斯消去法是由20世纪40年代由Dwyer 、冯 偌伊曼（von Neumann ）等人提出的. 高斯最早提出，通过积分点的最佳分布可得到更为精确的数值积分方法，并于1814年发表了第一个这样的求积公式，高斯求积公式因此得名. 高斯名言 ● 数学是科学里的皇后，数论是数学中的女王. ● 宁可少些，但要好些. ● 数学中的一些美丽定理具有这样的特性: 它们极易从事实中归纳出来，但证明却隐 藏得极深. ● 给我最大快乐的，不是已懂得的知识，而是不断的学习；不是已有的东西，而是不 断的获取；不是已达到的高度，而是继续不断的攀登. 练习题 1. 确定下列求积公式中积分系数或积分节点的待定值，使其代数精度尽量高，并指明所构造的求积公式所具有的代数精度： (1)∫ f (x )dx ≈ A ?1f (??)+ A 0f (0)+ A 1f(?)2? ?2? (2)∫f (x )dx ≈ 1 ?1[f (?1)+2f (x 1)+3f(x 2)]3? 2. 若积分节点x k ,k =0,…,n 给定，要求机械求积公式的积分系数，使得求积公式至少具有n 次代数精度. 请根据定理7.1列出待求解的线性方程组，并判断解的存在性和唯一性. 3. 直接验证柯特斯公式： C = b ?a [7f (x 0)+32f (x 1)+12f (x 2)+32f (x 3)+7f (x 4)] 具有5次代数精度. 4. 用辛普森公式求积分∫e ?x dx 10并估计误差. 5. 证明下列等式，它们分别说明了三种矩形求积公式及其余项公式. (1).∫f (x )dx =(b ?a )f (a )+ f ′(η)2(b ?a)2b a ,η∈(a,b); (2).∫f (x )dx =(b ?a )f (b )? f ′(η)2(b ?a)2b a ,η∈(a,b) ; (3).∫f (x )dx =(b ?a )f (a +b )+ f ′′(η)(b ?a)3b a ,η∈(a,b) . 6. 对下列积分，分别用复合梯形公式和复合辛普森公式计算, 其中n 表示计算中使用n +1个区间等分点上的函数值，然后比较两种方法计算结果的准确度. (1).∫x 4+x 2dx,n =810 (2).∫√xdx, n =491 第四章 数值积分与数值微分 1.确定下列求积公式中的特定参数，使其代数精度尽量高，并指明所构造出的求积公式所具有的代数精度： 101210121 12120 (1)()()(0)(); (2)()()(0)(); (3)()[(1)2()3()]/3; (4)()[(0)()]/2[(0)()]; h h h h h f x dx A f h A f A f h f x dx A f h A f A f h f x dx f f x f x f x dx h f f h ah f f h -----≈-++≈-++≈-++''≈++-?? ?? 解： 求解求积公式的代数精度时，应根据代数精度的定义，即求积公式对于次数不超过m 的多项式均能准确地成立，但对于m+1次多项式就不准确成立，进行验证性求解。 （1）若101(1) ()()(0)()h h f x dx A f h A f A f h --≈-++? 令()1f x =，则 1012h A A A -=++ 令()f x x =，则 110A h A h -=-+ 令2 ()f x x =，则 3 221123 h h A h A -=+ 从而解得 011431313A h A h A h -?=?? ? =?? ?=?? 令3 ()f x x =，则 3()0h h h h f x dx x dx --==? ? 101()(0)()0A f h A f A f h --++= 令4 ()f x x =，则 455 1012()5 2 ()(0)()3 h h h h f x dx x dx h A f h A f A f h h ---== -++=? ? 故此时， 101()()(0)()h h f x dx A f h A f A f h --≠-++? 故 101()()(0)()h h f x dx A f h A f A f h --≈-++? 具有3次代数精度。 （2）若 21012()()(0)()h h f x dx A f h A f A f h --≈-++? 令()1f x =，则 1014h A A A -=++ 令()f x x =，则 110A h A h -=-+ 令2 ()f x x =，则 3 2211163 h h A h A -=+ 从而解得 1143 8383A h A h A h -?=-?? ? =?? ?=?? 令3 ()f x x =，则 22322()0h h h h f x dx x dx --==? ? 101()(0)()0A f h A f A f h --++= 实验三 数值积分程序设计算法 1）实验目的 通过本次实验熟悉并掌握各种数值积分算法及如何在matlab 中通过设计程序实现这些算法，从而更好地解决实际中的问题。 2）实验题目 给出积分 dx x I ? -= 3 2 2 1 1 1.用Simpson 公式和N=8的复合Simpson 公式求积分的近似值. 2.用复合梯形公式、复合抛物线公式、龙贝格公式求定积分，要求绝对误差为 7 10*2 1-= ε，将计算结果与精确解做比较，并对计算结果进行分析。 3）实验原理与理论基础 Simpson 公式 )]()2 ( 4)([6 b f b a f a f a b S +++-= 复化梯形公式 将定积分? = b a dx x f I )(的积分区间],[b a 分隔为n 等分，各节点为 n j jh a x j ,,1,0, =+= n a b h -= 复合梯形(Trapz)公式为 ])()(2)([21 1 ∑-=++-= n j j n b f x f a f n a b T 如果将],[b a 分隔为2n 等分，而n a b h /)(-=不变， 则 )]()(2)(2)([41 2 111 2b f x f x f a f n a b T n j j n j j n +++-= ∑∑-=+-= 其中 h j a h x x j j )2 1(2 12 1+ +=+ =+ ，)]()(2)(2)([41 2 11 1 2b f x f x f a f n a b T n j j n j j n +++-= ∑∑-=+ -= ∑ -=-++-+ =1 )2) 12((22 1n j n n a b j a f n a b T n=1时，a b h -=，则)]()([2 1b f a f a b T +-= )0(0T = )2 1(2 2 112h a f a b T T + -+ =)1(0T = 若12-=k n ，记)1(0-=k T T n ， ,2,1=k 1 2 --= k a b h jh a x j +=1 2 --+=k a b j a h x x j j 2 12 1+ =+ k a b j a 2 ) 12(-++=，则可得如下递推公式 专题六数值微积分与方程求解6.1 数值微分与数值积分 ?数值微分 ?数值积分 1．数值微分 （1）数值差分与差商 微积分中，任意函数f(x)在x 0点的导数是通过极限定义的： h x f h x f x f h )()(lim )('0 0-+=→h h x f x f x f h ) ()(lim )('0 00 0--=→h h x f h x f x f h ) 2/()2/(lim )('0 --+=→ ) ()()(000 x f h x f x f -+=?) ()()(0 h x f x f x f --=?) 2/()2/()(0 h x f h x f x f --+=δ如果去掉极限定义中h 趋向于0的极限过程，得到函数在x 0点处以h （h>0）为步长的向前差分、向后差分和中心差分公式： 向前差分： 向后差分： 中心差分： 函数f(x)在点x 0的微分接近于函数在该点的差分，而f 在点x 的导数接近于函数在该点的差商。 h x f h x f x f ) ()(≈ )('0 00 -+h h x f x f x f ) ()(≈ )('0 00 --h h x f h x f x f ) 2/()2/(≈ )('0 --+向前差商： 向后差商： 中心差商： 当步长h 充分小时，得到函数在x 0点处以h （h>0）为步长的向前差商、 向后差商和中心差商公式： （2）数值微分的实现 MATLAB提供了求向前差分的函数diff，其调用格式有三种： ?dx=diff(x)：计算向量x的向前差分，dx(i)=x(i+1)-x(i)，i=1，2，…，n-1。?dx=diff(x,n)：计算向量x的n阶向前差分。例如，diff(x,2)=diff(diff(x))。?dx=diff(A,n,dim)：计算矩阵A的n阶差分，dim=1时（默认状态），按列计算差分；dim=2，按行计算差分。 注意：diff函数计算的是向量元素间的差分，故差分向量元素的个数比原向量少了一个。同样，对于矩阵来说，差分后的矩阵比原矩阵少了一行或一列。 另外，计算差分之后，可以用f(x)在某点处的差商作为其导数的近似值。 数值积分与 数值微分 习题课 一、已知012113,,424x x x ===，给出以这 3个点为求积节 点在[]0.1上的插值型求积公式 解：过这3个点的插值多项式基函数为 ()()()()()()()()()()()()()()()()1202010202121012012220211 20,0,1,2 k k x x x x l x x x x x x x x x l x x x x x x x x x l x x x x x A l x dx k --= ----= ----= --==? ()()()()()()()()()()()()111200001021102100101210120202113224111334244131441113324241142x x x x x x A dx dx x x x x x x x x x x A dx dx x x x x x x x x x x A dx x x x x ????-- ???--????=== --????-- ??? ???? ????-- ???--????===- --????-- ??? ???? ????-- ??--???==--?????102313134442dx ??= ????-- ??? ???? ? 故所求的插值型求积公式为 ()1 211 123343234f x dx f f f ??????≈- + ? ? ??????? ? 二、确定求积公式 ( )( )(1 1158059f x dx f f f -? ?≈++?? ? 的代数精度，它是Gauss 公式吗？ 证明：求积公式中系数与节点全部给定，直接检验 依次取()23451,,,,,f x x x x x x =，有 [ ](1 1111 215181519 1058059dx xdx --==?+?+???==?+?+?? ??? M a t l a b数值积分与数值微分 Matlab数值积分 1.一重数值积分的实现方法 变步长辛普森法、高斯-克朗罗德法、梯形积分法 1.1变步长辛普森法 Matlab提供了quad函数和quadl函数用于实现变步长 辛普森法求数值积分.调用格式为: [I,n]=Quad(@fname,a,b,tol,trace) [I,n]=Quadl(@fname,a,b,tol,trace) Fname是函数文件名，a,b分别为积分下限、积分上限； tol为精度控制，默认为1.0×10-6，trace控制是否展 开积分过程，若为0则不展开，非0则展开，默认不展开. 返回值I为积分数值；n为调用函数的次数. --------------------------------------------------------------------- 例如：求 ∫e0.5x sin(x+π )dx 3π 的值. 先建立函数文件 fesin.m function f=fesin(x) f=exp(-0.5*x).*sin(x+(pi/6));再调用quad函数 [I,n]=quad(@fesin,0,3*pi,1e-10) I= 0.9008 n= 365 --------------------------------------------------------------------- 例如：分别用quad函数和quadl函数求积分 ∫e0.5x sin(x+π 6 )dx 3π 的近似值，比较函数调用的次数. 先建立函数文件 fesin.m function f=fesin(x) f=exp(-0.5*x).*sin(x+(pi/6)); formatlong [I,n]=quadl(@fesin,0,3*pi,1e-10) I= n= 198 [I,n]=quad(@fesin,0,3*pi,1e-10) I= n= 365 --------------------------------------------------------------------- 可以发现quadl函数调用原函数的次数比quad少，并 且比quad函数求得的数值解更精确. 1.2高斯-克朗罗德法 第4章附录 4.2.2 复化求积分 例题4.2.5计算程序 //simp.c// # include subroutine trap(a,b,f,eps,t,n) real(8) a,b,f,t,fa,fb,h,t0,s,x fa=f(a); fb=f(b) n=1; h=b-a t0=h*(fa+fb)/2.0 5 s=0.0 do 10 k=0,n-1 x=a+(k+0.5)*h s=s+f(x) 10 continue t=(t0+h*s)/2.0 if (abs(t-t0).ge.eps) then t0=t n=n+n h=h/2.0 goto 5 end if return end %%% demo_aTrapInt.m %%% function demo_aTrapInt clc;clear; format long; [T nsub] = aTrapInt(@f01,0,1,0.000001) end function [T nsub]= aTrapInt(f,a,b,eps) tol = 1; nsub = 1; inall = 0; T = 0.5*(b-a)*(f(a)+f(b)); while tol > eps T0 = T; nsub = 2*nsub; n = nsub + 1; % total number of nodes h = (b-a)/nsub; % stepsize x = a:h:b; % divide the interval inall = inall + sum(f(x(2:2:n-1))); T = 0.5*h * (f(a)+2*inall+f(b)); tol = abs(T-T0); end end 第五章 数值微分与数值积分 一．分别用向前差商，向后差商和中心差商公式计算()f x =2x =的导数的近似值。其中，步长0.1h =。 【详解】 00()()(20.1)(2)=0.349 2410.10.1 f x h f x f f h +?+?===向前差商 00()()(2)(20.1)=0.358 0870.10.1 f x f x h f f h ????===向后差商 00()()(20.1)(20.1)= 0.353 664220.10.2f x h f x h f f h +??+??===×中心差商 二．已知数据 x 2.5 2.55 2.60 2.65 2.70 ()f x 1.58114 1.59687 2 1.62788 1.64317 求( 2.50),(2.60),(2.70)f f f ′′′的近似值。 【详解】 0.05h =，按照三点公式 3(2.50)4(2.55)(2.60)3 1.581144 1.59687 1.61245(2.50)0.316 10020.050.1 f f f f ?+??×+×?′≈==×(2.65)(2.55)1.627881.59687(2.60)0.310 10020.050.1 f f f ??′≈==× (2.60)4(2.65)3(2.70)241.6278831.64317(2.70) 4.179 90020.050.1 f f f f ?+?×+×′≈==× 三．已知如下数据 x 3 4 5 6 7 8 ()f x 2.937 6 6.963 213.600 0 23.500 8 37.318 4 55.705 6 数值分析实验报告四 数值积分与数值微分实验（2学时） 一 实验目的 1．掌握复化的梯形公式、Simpson 公式等牛顿-柯特斯公式计算积分。 2. 掌握数值微分的计算方法。 二 实验内容 1. 用复化梯形公式计算积分。 ?9 0dx x M=8 2. 用复化Simpson 公式计算积分。 ? 90dx x M=8 3. 给定下列表格值 利用四点式（n=3）求)50()50('''f f 和的值。 三 实验步骤（算法）与结果 1复化梯形公式 用C 语言编程如下： #include float y; y=sqrt(x); return y; } void main() { int i,m; float a,b,h,r; printf("输入等分数m:" ); scanf("%d",&m); printf("输入区间左端点a的值:"); scanf("%f",&a); printf("输入区间右端点b的值:"); scanf("%f",&b); float x[m+1]; h=(b-a)/m; for(i=0;i<=m;i++) x[i]=a+i*h; r=0; for(i=0;i<=m;i++) {if(i==0) r=r+h*0.5*f(x[i]); if(i>0&&i if(i==m) r=r+0.5*h*f(x[i]); } printf("输出区间[%3.1f %3.1f]的积分值:%f\n",a,b,r); } 求解结果如下： 输入等分数m:8 输入区间左端点a的值:0 输入区间右端点b的值:9 输出区间[0.0 9.0]的积分值:17.769514 2复化Simpson公式 用C语言编程如下： #include 数值积分与微分 摘要 本文首先列举了一些常用的数值求积方法，一是插值型求积公式，以N e w t o n C o t e s -公式为代表，并分析了复合型的Newton Cotes -公式；另一个是Gauss Ledendre -求积公式，并给出几个常用的Gauss Ledendre -求积公式。其次，本文对数值微分方法进行分析，主要是差分型数值微分和插值型数值微分，都给出了几种常用的微分方法。然后，本文比较了数值积分与微分的关系，发现数值积分与微分都与插值或拟合密不可分。 本文在每个方法时都分析了误差余项，并且在最后都给出了MATLAB 的调用程序。 关键词：插值型积分Gauss Ledendre -差分数值微分插值型数值微分 MATLAB 一、常用的积分方法 计算积分时，根据Newton Leibniz -公式， ()()()b a f x dx F b F a =-? 但如果碰到以下几种情况： 1）被积函数以一组数据形式表示； 2）被积函数过于特殊或者原函数无法用初等函数表示 3）原函数十分复杂难以计算 这些现象中，Newton Leibniz -公式很难发挥作用，只能建立积分的近似计算方法，数值积分是常用的近似计算的方法。 1.1 插值型积分公式 积分中的一个常用方法是利用插值多项式来构造数值求积公式，具体的步骤如下： 在积分区间上[,]a b 上取一组节点：01201,,,,()n n x x x x a x x x b ≤<<≤ 。已知()k x f 的函数值，作()x f 的n 次插值多项式，则 (1) ()10()()()()() (1)!n n x n n k k n k f f L x R x f x l x w x n ++==+=++∑ 其中，()k l x 为n 次插值基函数，则得 (1)+10()(()())1 =[()]()[()](1)!b b n n a a n b b n k k n a a k f x dx L x R x dx l x dx f x f x w x dx n ξ+==+++? ?∑??（） 公式写成一般形式： ()()[]n b k k n a k f x dx A f x R f ==+∑? 其中， 01100110 ()()()() ()()()()()b b k k k k a a k k k k k k x x x x x x x x A l x dx dx x x x x x x x x -+-+----==----?? (1)+11 [][()](1)!b n n n a R f f x w x dx n ξ+=+?（） 显然，当被积函数f 为次数小于等于n 的多项式时，其相应的插值型求积公式为准确公式，即： ()() n b k k a k f x dx A f x ==∑? 1.1.1 求积公式的代数精度 定义：求积公式对于任何次数不大于m 的代数多项式()f x 均精确成立，而对于 1()m f x x +=不精确成立，则称求积公式具有m 次代数精度。 定理：含有1n +个节点(0,1,,)k x k n = 的插值型求积公式的代数精度至少为n 。 第七章数值积分与数值微分 积分问题最早来自于几何形体的面积、体积计算，也是经典力学中的重要问题（例如计算物体的重心位置）. 在现实应用中，很多积分的结果并不能写成解析表达式，因此需要通过数值方法来计算. 数值微分是利用一些离散点上的函数值近似计算某一点处的函数导数，它针对表达式未知的函数. 本章介绍一元函数积分（一重积分）和微分的各种数值算法，它们也是数值求解积分方程、微分方程的基础. 7.1数值积分概论 7.1.1基本思想 考虑如下定积分的计算： I(f)≡∫f(x)dx b a ,(7.1) 其中函数f: ?→?，首先应想到的是微积分中学习过的牛顿-莱布尼兹（Newton-Leibniz）公式： ∫f(x)dx b a =F(b)?F(a) , 其中F′(x)=f(x),即F(x)为f(x)的原函数. 但是，诸如e x2,sinx x ,sinx2等表达式很简单的函数却找不到用初等函数表示的原函数，因此必须研究数值方法来近似计算积分. 另一方面，某些函数的原函数虽然可以解析表示，但其推导、计算非常复杂，此时也需要使用数值积分方法. 一般考虑连续的、或在区间[a,b]上可积①的函数f(x)，则根据积分的定义有： lim n→∞,?→0∑(x i+1?x i)f(ξi) n i=0 =I(f) , (7.2) 其中a=x0 第5章数值积分与数值微分方法关于定积分计算，已经有较多方法，如公式法、分步积分法等，但实际问题中，经常出现不能用通常这些积分方法计算的定积分问题。怎样把这些通常方法失效的定积分在一定精度下快速计算出来，特别是通过计算机编程计算出来就是本章研究的内容。 此外，怎样根据函数在若干个点处的函数值去求该函数的导数近似值也是本章介绍的内容。 本章涉及的方法有Newton-Cotes求积公式、Gauss求积公式、复化求积公式、Romberg求积公式和数值微分。 5.1 引例 人造地球卫星轨道可视为平面上的椭圆。我国的第一颗人造地球卫星近地点距离地球表面439km ，远地点距地球表面2384km ，地球半径为6371km ，求该卫星的轨道长度。 本问题可用椭圆参数方程 cos ,,0sin x a t a b y b t π=?≤≤>?=? (0t 2) 来描述人造地球卫星的轨道，式中a, b 分别为椭圆的长短轴，该轨道的长度L 就是如下参数方程弧长积分 但这个积分是椭圆积分，不能用解析方法计算。 5.2问题的描述与基本概念 要想用计算机来计算()b a f x dx ?，应对其做离散化处 理。注意到定积分是如下和式的极限 1 ()lim ()n b i i a i f x dx f x λξ→==?∑? 要离散化，做 1） 去掉极限号lim 2） 将i ξ取为具体的i x 值 3） 为减少离散化带来的误差，将i x ?用待定系数i A 代替 于是就得到 定义 5.1 若存在实数1212,,,;,, ,,n n x x x A A A 且任 取()[,],f x C a b ∈都有 1 ()()n b i i a i f x dx A f x =≈∑? （5.1） 解：卫星轨道的示意图如右上图所示，,a b 分别是椭圆轨道的长半轴和短半轴，地球位于椭圆的一个焦点处，焦距为c ，地球半径为r ，近地点和远地点与地球表面的距离分别是1s 和2s . 由图中可知，上述数据存在如下关系： 12122,a s s r c a s r =++=-- 由椭圆性质 b =，将12,,s s r 的数据代入以上各式可得7782.5a km =，7721.5b km =. 椭圆的参数方程为： c o s ,s i n x a t y b t == ， (02)t π≤< 根据计算参数方程弧长的公式，椭圆长度可表为如下积分： /2 22221/20 4(sin cos )L a t b t dt π=+? 由于该积分无法求得解析解，下面我们编写MATLAB 程序对其进行数值求解。 s1=439;s2=2384;r=6371; a=(s1+s2)/2+r a = 7.7825e+003 >> c=a-s1-r; >> b=sqrt(a^2-c^2) b = 7.7215e+003 y=inline('sqrt(7782.5^2*sin(t).^2+7721.5^2*cos(t).^2)'); %建立积分内联函数 >> t=0:pi/10:pi/2; y1=y(t); format long >> L1=4*trapz(t,y1) %梯形积分 L1 = 4.870744099902405e+004 >> L2=4*quad(y,0,pi/2,1e-6) %辛普森积分 L2 = 4.870744099903280e+004 求解结果显示：两种方法计算求得的积分结果相当接近，轨道长度约为：4 4.8710km ?. 第4章数值积分与数值微分 计算定积分是科学技术和工程计算中经常遇到的问题，然而绝大多数求定积分的问题都不能通过被积函数的原函数来求解，而要用数值的方法来解决. 4．1问题的提出神舟十号载人飞船于2013年6月11日17时38分在酒泉卫星发射中心成功发射，6月13日13时18分与天宫一号成功对接，在轨飞行15天，其中12天与天宫一号组成组合体在太空中飞行，6月26日8时7分，神舟十号返回舱成功返回地面. 神舟十号载人飞船发射的初始轨道为近地点约200公里、远地点约330公里的椭圆轨道，对接轨道为距地约343公里的近圆轨道，飞行速度约为每秒7.9公里.试计算神舟十号载人飞船在椭圆轨道飞行一圈的公里数，进而可以计算在轨飞行的公里数. 这里主要是计算神舟十号的椭圆轨道的周长.如图4.1所示， 图4.1神舟十号的椭圆轨道 已知地球半径=r 6371km，近地点高度为=1s 200km、远地点高度为=2s 330km，则椭圆长半轴=++=2/)2(21s s r a 6636km ，半焦距=-=2/)(12s s c 65km ，由椭圆参数方程t b y t a x sin ,cos ==，其中22c a b -=，知椭圆周长为 dt t b t a L ?+=20 2222cos sin 4π dt t c a ?-=20 222cos 4π dt t a c a ?-=2 0222cos 14π 这是一个定积分，只要求出它的值就行了. 一般的，对于定积分 ()b a I f x dx =?如果被积函数()f x 在区间[,]a b 上连续，且()f x 的原函数为()F x ，则由牛顿—莱布尼兹（Newton-Leibnitz）公式，有 ()()() b a f x dx F b F a =-?似乎问题已经解决，其实不然. （1）有很多被积函数找不到用初等函数表示的原函数F(x)，例如 等等，表面看它们并不复杂，但却无法求得用初等函 数表示的原函数F(x).这里的椭圆积分dt t a c ?-2 0222cos 1π 也是这样，（2）有的积分即使能找到用初等函数表示的原函数F(x)，但原函数非常复杂，用牛顿—莱布尼兹公式，计算也很困难，例如 24211211ln arc (21)arc (21)1422122 x x dx tg x tg x C x x x ++??=+++-+??+-+?（3）有的被积函数()f x 是由测量或数值计算给出的数据表，是列表函数，也无法用牛顿—莱布尼兹公式计算. 对于上述这些情况，都要求建立定积分的近似计算方法——数值积分法.数值分析第四章数值积分与数值微分习题复习资料
数值积分与数值微分实验报告
数值微分与数值积分
数值积分与数值微分知识题课
Matlab数值积分与数值微分
第4章 数值微分与积分(附录)
数值微分与数值积分练习题
实验4_数值积分与数值微分
数值积分与微分方法
第7章 数值积分与数值微分
数值积分与数值微分
数值积分与数值微分 编程计算
数值积分与数值微分