数值分析方法第六章

Chapter5

数值微分与数值积分2先看一个实例:已知20世纪美国人口的统计数据为(单位:百万)年份1900191019201930194019501960197019801990人口76.092.0106.5123.2131.7150.7179.3204.0226.5251.4试计算美国20世纪的年增长率

§1数值微分3

在实际问题中，往往会遇到某函数f(x)是用表格表示的,用通常的导数定义无法求导,因此要寻求其他方法近似求导。常用的数值微分方法有:一.运用差商求数值微分运用插值函数求数值微分三.运用样条插值函数求数值微分四.运用数值积分求数值微分4运用差商求数值微分分析：微分和积分是一对互逆的数学运算。复习：导数是差商当时的极限，在精度要求不高的情形下，可以简单地取差商作为导数的近似值基本原理：用差商代替微分56向后差商近似计算：向前差商近似计算：中心差商近似计算：

差商公式7由Taylor展开误差向前差商8由Taylor展开误差向后差商9由Taylor展开误差中心差商10分析1、步长过大则截断误差显著，但如果步长太小又会导致舍入误差的增长。2、在实际计算时，在保证截断误差满足精度要求的前提下选取尽可能大的步长

3、事先给出一个合适的步长往往是很困难的。通常在在变步长的过程中实现步长的自动选择。11主要方法：插值多项式，讨论函数f(x)导数的近似值求法。

插值型求导公式

定义：若函数f(x)在节点处的函数值已知，作f(x)的n次插值多项式，并用近似代替f(x)，即由于是多项式，容易求其导数，故对应于f(x)的每一个插值多项式，建立一个数值微分公式这样建立起来的数值微分公式，称为插值型数值微分公式。12

带有余项的数值微分公式其中之间，上式第一项是的近似值。即使与f(x)的函数值处处相差甚微，但两个函数的导数在某些点上的值仍可能有很大的差异所以要对误差进行分析13125102050100025020015010050300annealingtimet,1000sSeebeckCoefficienta,mV/Ka010203040误差放大!!!14如果只计算结点处的导数值数值微分公式的余项为15

带余项的两点数值微分公式（n=1）16类似地，可以得到n=2的带余项的三点公式17类似地，可以得到n=2的带余项的二阶三点公式18利用三次样条插值函数构造数值微分公式对于给定函数表和适当地边界条件，有三次样条插值函数S(x)，并用S(x)近似代替f(x)，即由于S(x)是一个分段三次多项式，在各子区间上容易求出其导数，故建立一个数值微分公式.19注意数值微分公式20例：利用函数在节点上的函数值和边界条件S’(-1)=0.0740,S’(1)=-0.0740构造三次样条插值函数S(x)，并用它来计算f(x)和f’(x)在下列点处的近似值。2122MATLAB中的数值微分数值微分的实现在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))。23例设x由[0,2π]间均匀分布的10个点组成，求sinx的1~3阶差分。命令如下：X=linspace(0,2*pi,10);Y=sin(X);DY=diff(Y);%计算Y的一阶差分D2Y=diff(Y,2);%计算Y的二阶差分，也可用命令diff(DY)计算D3Y=diff(Y,3);%计算Y的三阶差分，也可用diff(D2Y)或diff(DY,2)24例用不同的方法求函数f(x)的数值导数，并在同一个坐标系中做出f'(x)的图像。程序如下：f=inline('sqrt(x.^3+2*x.^2-x+12)+(x+5).^(1/6)+5*x+2');g=inline('(3*x.^2+4*x-1)./sqrt(x.^3+2*x.^2-x+12)/2+1/6./(x+5).^(5/6)+5');x=-3:0.01:3;p=polyfit(x,f(x),5);%用5次多项式p拟合f(x)dp=polyder(p);%对拟合多项式p求导数dpdpx=polyval(dp,x);%求dp在假设点的函数值dx=diff(f([x,3.01]))/0.01;%直接对f(x)求数值导数gx=g(x);%求函数f的导函数g在假设点的导数plot(x,dpx,x,dx,'.',x,gx,'-');%作图25§2数值积分在高等数学中，牛顿-莱布尼兹公式

F(x)是被积函数f(x)的原函数完美地解决了定积分问题。26

在工程技术和科学研究中，常常遇到如下情况：1.

f(x)的结果复杂，求原函数困难。2.

f(x)的原函数不存在，或不能用初等函数表示。3.

f(x)无函数式，只有函数表。27R1R2F颗粒尺寸，R分布密度，f例：计算R1

R

R2的颗粒数量(分数)原始方法：曲线下面面积所占方格数目数值积分雏形

数格子28当f

(x)

是闭区间上的连续函数时，如何计算下列定积分？积分中值定理定积分的计算29上述公式的几何意义底为而高为的矩形面积恰好等于所求曲边梯形的面积。问题：如何确定？I(f)Oxy30左矩形右矩形矩形中矩形定积分常用的近似公式左矩形公式：中矩形公式：右矩形公式：

若取两点，并令，则可得梯形公式（两点求积公式）则可得Simpson公式(三点求积公式)

若取三点，并令37

线性近似：梯形公式ab梯形近似xyy=f(x)简单，但是

ab梯形近似xyy=f(x)38

抛物线近似：辛普森(Simpson)公式xyy=f(x)ab抛物线近似通常，抛物线近似优于梯形近似但是，有特例

abxyy=f(x)39利用梯形公式，Simposn公式公式计算，并与精确值进行比较。解：梯形公式Simposn公式精确值

一般地，取区间内个点处的高度通过加权平均的方法近似地得出平均高度这类求积方法称为机械求积:§3插值型求积公式

或写成:数值积分公式求积系数

求积节点记称为数值求积公式称为求积公式余项(误差).（1）（2）构造或确定一个求积公式，要解决的问题包括:(i)

确定求积系数

和求积节点

(iii)

求积公式的误差估计和收敛性分析.(ii)确定衡量求积公式好坏的标准；定义：在积分区间上，取个节点作

的次代数插值多项式（拉格朗日插值公式）:则有其中，为插值余项。于是有：取Ak由节点决定，与

无关。称为插值型求积公式

称求积公式具有m次代数精度,如果它满足如下两个条件:定义：(i)对所有次数≤m次的多项式,有(ii)存在m+1次多项式,使得求积公式的代数精度上述定义中的条件(i),(ii)等价于:注：梯形公式与中矩形公式都只具有1次代数精度。一般的，若要使求积公式（1）具有m次代数精度，则只要使求积公式对f(x)=1，x，x2，…,xm都准确成立，即代数精度推论求积系数满足:

形如的求积公式至少有n

次代数精度

该公式为插值型（即：）定理截断误差51几个常用的求积公式的代数精度1.梯形公式的代数精度52532.辛普森公式的代数精度5455可以证明对等式不成立辛普森公式具有三次代数精度56n+1个结点的插值型数值积分公式至少有n阶代数精度。因此对于次数不大于n的多项式其余项57

求插值型求积公式

并确定其代数精度。分析：实际上该题目是求A0，A1，并确定其代数精度。从而求积公式为有2个结点,因而代数精度大于等于1，从m=2开始验证代数精度解§4Newton-Cotes公式一、Cotes系数取节点为等距分布：由此构造的插值型求积公式称为Newton-Cotes公式,此时求积系数：令Cotes系数x–xj=(t-j)hdx=hdt二、Newton-Cotes公式1、定义：记则求积公式变为称上式为n阶闭型Newton-Cotes求积公式。注意:由式确定的Cotes系数只与和有关,

与

和积分区间无关，且满足:61

科特斯系数表62且满足

注：Cotes系数仅取决于n

和i，可通过查表得到。与被积函数f(x)

及积分区间[a,b]均无关。

牛顿-科特斯公式(2)(1)63梯形公式

当n=1,x0=a,x1=b时，有

64辛普森公式

当n=2,x0=a,x1=(a+b)/2,x2=b时，有65柯特斯公式其中66

定理：若在[a,b]上连续，梯形公式的截断误差为若在[a,b]上连续，辛普生公式的截断误差为若在[a,b]上连续，柯特斯公式的截断误差为67实际计算中，当积分区间较大时，直接使用这些积分公式精度难以保证。

为了提高计算精度，采用复合求积的方法。

将区间[a,b]适当分割成若干个子区间，对每个子区间使用求积公式，构成所谓的复化求积公式，这是提高积分精度的一个常用的方法。68xyy=f(x)ab将[a,b]

分成若干小区间[xk-1,xk],(k=1,2,,n)x0x1x2xnxn-1h前提：节点函数值f(xk)

已知复合梯形公式69将积分区间[a,b]划分为n等分,记子区间

的中点为，在每个小区间上应用辛普森公式，则有记

复合辛普森公式70复合柯特斯公式其中71误差分析A、复合梯形公式的误差B、复合Simpson公式的误差72C、复合Cotes公式的误差73

用复化梯形公式计算定积分才能使误差不超过解:取,则,又区间长度b-a=1，对复化梯形公式有余项即,n≥212.85，取n=213，即将区间[0,1]分为213等份时，用复化梯形公式计算误差不超过

。

问区间[0,1]应分多少等份74复合求积方法对提高精度是一种行之有效的方法，但是在使用求积公式之前必须先给出步长。步长取得太大精度难以保证，步长太小则会导致计算量的增加，而事先给出一个合适的步长往往是困难的。为了解决这个问题，较好的方法就是采用变步长的计算方法即在步长逐次减半的过程中，反复利用复合求积公式进行计算，直到所求得的积分值满足精度要求为止。7576

设将积分区间[a,b]n等分，即分成n个子区间，一共有n+1个节点，即x=a+kh,k=0,1,…，n，步长。对于某个子区间,利用梯形公式计算积分近似值有

对整个区间[a,b]有77将子区间再二等份,取其中点作新节点,此时区间数增加了一倍为2n,对某个子区间,利用复化梯形公式计算其积分近似值。对整个区间[a,b]有比较和有78

当把积分区间分成n等份，用复化梯形公式计算积分I的近似值时，截断误差为

若把区间再分半为2n等份，计算出定积分的近似值，则截断误差为

当在区间[a,b]上变化不大时,有所以79

可见,当步长二分后误差将减至,将上式移项整理，可得验后误差估计式

上式说明，只要二等份前后两个积分值和相当接近，就可以保证计算结果的误差很小，使接近于积分值I。80变步长的梯形求积算法实现（1）变步长的梯形求积法的计算步骤①变步长梯形求积法。它是以梯形求积公式为基础，逐步减少步长，按如下递推公式求二分后的梯形值其中Tn和T2n分别代表二等分前后的积分值

②如果,(ε为给定的误差限)则T2n作为积分的近似值,否则继续进行二等分,即转①再计算，直到满足所要求的精度为止，最终取二分后的积分值T2n

作为所求的结果81例用变步长梯形求积法计算定积分解:先对整个区间

0,1

用梯形公式,对于

所以有然后将区间二等份,由于,故有

进一步二分求积区间,并计算新分点上的函数值

82有

这样不断二分下去。积分的准确值为0.946083。

83龙贝格求积公式

变步长梯形求积法算法简单，但精度较差，收敛速度较慢，但可以利用梯形法算法简单的优点，形成一个新算法，这就是龙贝格求积公式。龙贝格公式又称逐次分半加速法。根据积分区间分成n等份和2n等份时的误差估计式可得

积分值的误差大致等于,如果用对进行修正时，与之和比更接近积分真值84考察与n等份辛普森公式之间的关系。将复化梯形公式

梯形变步长公式代入表达式得

故85

(6.9)

这就是说，用梯形法二分前后两个积分值和作线性组合，结果却得到复化辛普森公式计算得到的积分值。

86再考察辛普森法。其截断误差与成正比，因此，如果将步长折半，则误差减至，即有

由此可得

可以验证,上式右端的值其实等于Cn，就是说，用辛普森公式二等份前后的两个积分值Sn和S2n

作线性组合后，可得到柯特斯公式求得的积分值Cn，即有87用同样的方法，根据柯特斯公式的误差公式，可进一步导出龙贝格公式

在变步长的过程中运用上述公式，就能将粗糙的梯形值Tn逐步加工成精度较高的辛普森值Sn、柯特斯值Cn和龙贝格值Rn或者说，将收敛缓慢的梯形值序列Tn加工成收敛迅速的龙贝格值序列Rn，这种加速方法称为龙贝格算法88龙贝格求积法计算步骤用梯形公式计算积分近似值按变步长梯形公式计算积分近似值将区间逐次分半,令区间长度

计算③按加速公式求加速值梯形加速公式：辛普森加速公式：

龙贝格求积公式：89④精度控制；直到相邻两次积分值

（其中ε为允许的误差限）则终止计算并取Rn作为积分的近似值，否则将区间再对分，重复②，③，④的计算，直到满足精度要求为止。90例

用龙贝格算法计算定积分要求相邻两次龙贝格值的偏差不超过解:由题意

9192由于，于是有

前面介绍的n+1个节点的Newton-Cotes求积公式，其特征是节点是等距的。这种特点使得求积公式便于构造，复化求积公式易于形成。但同时也限制了公式的精度。n是偶数时，代数精度为n+1，n是奇数时，代数精度为n

。

我们知道n+1个节点的插值型求积公式的代数精确度不低于n

。设想：能不能在区间[a,b]上适当选择n+1个节点x0x1,x2,……,xn,使插值求积公式的代数精度高于n？

答案是肯定的，适当选择节点，可使公式的精度最高达到2n+1，这就是本节所要介绍的高斯求积公式。5高斯(Gauss)求积公式考虑更一般形式的数值积分问题定义：若求积公式对一切不高于m次的多项式p(x)都等号成立，即R(p)=0;而对于某个m+1次多项式等号不成立，则称此求积公式的代数精度为m.一、构造高斯型求积公式的基本原理和方法

数值分析定理1：设节点x0,x1…,xn∈[a,b]，则求积公式

的代数精度最高为2n+1次。

分别取

f(x)=1,x，x2，...xr

代入公式，并让其成为等式，得：

A0+A1+……+An=∫ab1dx.=b-ax0A0+x1A1+……+xnAn=∫abxdx.=（b2-a2)/2......x0

rA0+x1

rA1+……+xn

rAn=∫abxrdxr

=(br+1-ar+1)

(r+1)数值分析

事实上,取

2n+2次多项式g(x)=(x-x0)2(x-x1)2….(x-xn)2

代入求积公式,这里x0,x1…,xn是节点，有左右,故等式不成立,求积公式的代数精度最高为2n+1次。证毕.

上式共有r+1个等式，2n+2个待定系数(变元),要想如上方程组有唯一解，应有方程的个数等于变元的个数,即r+1=2n+2,这样导出求积公式的代数精度至少是2n+1,下面证明代数精度只能是2n+1.

定义:使求积公式达到最高代数精度2n+1的求积公式称为Guass求积公式。Guass求积公式的节点xk称为Guass点,系数Ak称为Guass系数.因为Guass求积公式也是插值型求积公式,故有结论:

n+1个节点的插值型求积公式的代数精度d

满足：

n

d2n+1。数值分析例：选择系数与节点，使求积公式（1）成为Gauss公式。解：n=1,由定义，若求积公式具有3次代数精度，则其是Gauss公式。为此，分别取

f(x)=1,x，x2，x3

代入公式，并让其成为等式，得c1+

c2=2c1x1+

c2x2=0c1x12+

c2x22=2/3c1x13+

c2x23=0求解得：所求Gauss公式为：(1)用待定系数法构造高斯求积公式

设Pn(x)，n=0,1,2,…,为正交多项式序列，Pn(x)具有如下性质：1）对每一个n,Pn(x)是n次多项式。n=0,1,…2）(正交性)3）对任意一个次数≤n-1的多项式P(x)，有4）Pn(x)在(a,b)内有n个互异零点。（2）利用正交多项式构造高斯求积公式数值分析定理2

设x0,x1,…,xn是n+1次正交多项式Pn+1(x)的n+1

个零点,则插值型求积公式是Guass型求积公式。证明：只要证明求积公式的代数精确度为2n+1,即对任意一个次数≤2n+1的多项式求积公式都精确成立。设f(x)为任意一个次数≤2n+1的多项式，则有

f(x)=q(x)Pn+1(x)+r(x)，满足f(xk)=r(xk)这里，Pn+1(x)是n+1次正交多项式，q(x)、r(x)均是次数≤n的多项式。数值分析数值分析由性质(3)及(4)式，有由于n+1个节点的插值型求积公式的代数精确度不低于n，故有即对f(x)为任意一个次数≤2n+1的多项式求积公式都精确成立。证毕数值分析数值分析利用正交多项式构造高斯求积公式的基本步骤：代入积分式因此，求积系数为数值分析数值分析数值分析数值分析数值分析数值分析常用的高斯求积公式1.Gauss-Legendre求积公式

(1)其中高斯点为Legendre多项式的零点Guass点xk,Guass系数Ak都有表可以查询.数值分析数值分析数值分析数值分析数值分析数值分析数值分析数值分析一般区间的Gauss-Legendre求积公式

如果积分区间是[a,b]，用线性变换

这样就可以用Gauss-Legendre求积公式计算一般区间的积分.将积分区间从[a,b]变成[-1,1],由定积分的换元积分法有数值分析数值分析数值分析数值分析数值分析数值分析数值分析数值分析数值分析数值分析例

利用高斯求积公式计算解:

令x=1/2

(1+t),

则用高斯-Legendre求积公式计算.取n=4

积分精确值为

I=ln2=0.69314718…由此可见，高斯公式精确度是很高的.数值分析数值分析例:分别用不同方法计算如下积分,并做比较各种做法比较如下：1、用Newton-Cotes公式当n=1时，即用梯形公式，I≈0.9270354当n=2时,即用Simpson公式,I≈0.9461359当n=3时,I≈0.9461090当n=4时,I≈0.9460830当n=5时,I≈0.9460830I准=0.9460831数值分析数值分析2:用复化梯形公式

令h=1/8=0.1253：用复化辛卜生公式令h=1/8=0.125I准=0.9460831数值分析4、用Romberg公式KTn

SnCnRn00.920735510.93979330.946145920.94451350.94608690.940083030.94569060.94608330.94608310.9460831

I准=0.9460831数值分析5、用Gauss公式解：令x=(t+1)/2,

I准=0.9460831（2）用3个节点的Gauss公式（1）用2个节点的Gauss公式数值分析算法比较此例题的精确值为0.9460831..由例题的各种算法可知：对Newton-cotes公式，当n=1时只有1位有效数字，当n=2时有3位有效数字，当n=5时有7位有效数字。对复化梯形公式有2位有效数字，对复化辛卜生公式有6位有效数字。用复合梯形公式，对积分区间[0，1]二分了11次用2049个函数值，才可得到7位准确数字。Romberg公式对区间二分3次，用了9个函数值，得到同样的结果。用Gauss公式仅用了3个函数值，就得到结果。一、一维积分的高斯公式

其中f(ξi)是被积函数在积分点ξi处的数值，Hi为加数系数，n为积分点数目。对于n个积分点，只要选取适当的加数系数及积分点位置，能够使式在被积分函数为不超过(2n-1)次多项式时精确成立。由于多数函数可表示成多项式形式，这种积分适应于大多数函数。

高斯积分法与有限单元法高斯积分法例如，n=1时不论f(ξ)的次数是0还是1，只需取H１=2，ξ1＝０，上式均是精确成立的。因为

高斯积分法 当n=2时，能保证式子精确成立所允许的多项式的最高次数是3，此时，f(ξ)的通式为其精确积分为数值积分为高斯积分法为了在C0～C3取任意值(包括取零值在内)时公式(f)是精确的，显然应有所以，应取，，高斯积分法n个插值结点非等距分布结点和积分权系数可以查表高斯积分法二维积分的高斯公式 以一维高斯积分公式为基础，导出二维及三维公式。求二维重积分

的数值时，可以先对ξ、η进行积分，

或改写成这就是二维的高斯积分公式。高斯积分法三维积分的高斯公式同样，可以求得三维高斯积分公式：中的n，m，l是分别关于变量ξ，η，ζ的积分点数目。 各个维数上的积分点数目由各个自变量在被积函数中可能出现的最高次数分别决定，一般并不要求相同。但为应用方便，常常在各个方向取相同的积分数，即统一为最高值高斯积分法由前面的推导可见，当在每个方向取n个积分点时，只要多项式被积函数中自变量的次数m≤2n-1，则用高斯求积公式求得的积分值是完全精确的。反过来，对于m次多项式的被积函数，为了积分值完全精确，积分点的数目必须取。高斯积分法高斯积分方法预先定义了积分点和相应的加权系数，求出被积分的函数在指定积分点上的数值，加权后求和，就得到了该函数的积分。高斯积分方法具有最高的计算精度。采用n个积分点的高斯积分可以达到2n-1阶的精度，也就是说，如果被积分的函数是2n-1次多项式，用n个积分点的高斯积分可以得到精确的积分结果。积分阶次的选择直接影响计算的精度和计算工作量。积分阶次的选择必须保证积分的精度。（完全精确积分）很多情况下，实际选取的高斯积分点数低于精确积分的要求，往往可以取得较完全精确积分更好的精度。（减缩积分）线性单元完全精确积分

二次单元减缩积分有限元分析主要步骤

我们知道，经过单元方程的组装以后，结构静力学有限元方程如下{F}=[K]{U}其中，{F}----节点载荷向量；[K]---总体刚度矩阵；{U}---节点位移向量在引入边界条件以后，解上述方程组，就可以得到节点位移向量{U}.这是求解结构静力学方程组所得到的第一组解，它是最精确的。得到节点的位移解后，下面是求取应变解和应力解。与位移解不同，它们并不是直接在节点上获得，而是首先在积分点上获得的。

有限元分析主要步骤所谓积分点是指，在对单元建立方程时，例如刚度矩阵是需要通过积分而得到的，而积分时为了能够方便计算，大多数有限元软件采用了所谓高斯积分的方式，即在单元内分布一些高斯点这样，有限元软件会首先获得这些高斯点的应力和应变，其方法如下：在高斯积分点上，依据几何方程:{ε}={B}{U}计算出高斯积分点上的应变:ε然后基于虎克定律及几何方程推导的结果来计算高斯积分点的应力。:{σ}={D}{B}{U}有限元分析主要步骤可见，在应变和应力计算方面，高斯积分点的应变和应力是最最准确的。利用特定单元的形函数以及高斯点的应力，应变值，将这些值外推到该单元的节点上，就得到了单元上节点的应力应变值。显然，不同的单元会共用一些节点，而从不同单元内的积分点外推到这些公共节点的应变值和应力值一般不相同,将一个公共节点的多个应力进行平均，以代表该节点的应力值。有限元分析主要步骤总之，求解节点应力的步骤是：（1）根据总体方程，得到节点的位移解。（2）根据几何方程，得到单元高斯点的应变解。（3）根据物理方程，得到单元高斯点的应力解。（4）在某一个单元内，基于形函数，将高斯点的应力外推到该单元的所有节点。（5）对于某一个公共节点，将该节点关联的所有单元所推出的该节点的应力解进行平均，最终得到该节点的应力解。积分点与节点的关系我们需要对应变在单元内的面积上进行积分时，因为节点的应力、位移显然与x,y无关，我们只需要考虑对形函数积分。

采用Gauss-Legendre多项式计算积分时，我们只需要计算根据特定积分点的值（在自然坐标系下是固定的，可以查手册，这些点也叫高斯点、积分点）并加以权重就可以。这就把复杂的积分问题变成了简单的代数问题。因为形函数只有单元有关，所以积分点也只与单元形状有关。

应力一般采用多个积分点的相互插值或外延来计算节点应力。这只是为了减少误差。因为在积分点应力比节点具有更高阶的误差。136MATLAB数值积分基于变步长辛普生法，MATLAB给出了quad函数来求定积分