数值计算06-数值积分与数值微分

1数值计算第六章数值积分与数值微分§6.1引言一、数值积分的必要性讨论如下形式的一元函数积分在微积分里，按Newton-Leibniz公式求定积分要求被积函数☞有解析表达式；☞

的原函数

为初等函数．

实际问题1.

的原函数

不能用初等函数表示考虑一个实际问题:建筑上用的一种铝制波纹瓦是用一种机器将一块平整的铝板压制而成的.假若要求波纹瓦长4英尺,每个波纹的高度(从中心线)为1英寸,且每个波纹以近似英寸为一个周期.求制做一块波纹瓦所需铝板的长度L.从到英寸间的弧长L.这个问题就是要求由函数给定的曲线,

由微积分学我们知道,所求的弧长可表示为:上述积分称为第二类椭圆积分。What’stheOriginalfunction?!It’ssocomplexthatwecannotgetit.类似的，下列函数也不存在由初等函数表示的原函数:2.

有些被积函数其原函数虽然可以用初等函数表示成有限形式,但表达式相当复杂,计算极不方便.例如函数:并不复杂,但它的原函数却十分复杂:3.

没有解析表达式，只有数表形式:1423454.5688.5原来通过原函数来计算积分有它的局限性。那……怎么办呢？呵呵…这就需要积分的数值方法来帮忙啦。二、数值积分的基本思想1、定积分的几何意义2、数值积分的理论依据依据积分中值定理,对于连续函数,在内存在一点,使得称

为在区间上的平均高度.3、求积公式的构造

若简单选取区间端点或中点的函数值作为平均高度，则可得一点求积公式如下：左矩形公式：中矩形公式：右矩形公式：左矩形公式：中矩形公式：右矩形公式：

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

若取三点，并令

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

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

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

确定求积系数

和求积节点

(iii)

求积公式的误差估计和收敛性分析.(ii)确定衡量求积公式好坏的标准；22求积公式的构造方法数值积分的基本问题

针对某类函数，选择合适的求积结点和求积系数，使得求积公式(1)具有尽可能小的截断误差或尽可能高的代数精度。

寻找被积函数f(x)在区间[a,b]上的一个容易求积的近似函数(如某个多项式)p(x)，然后用p(x)在区间[a,b]上的积分值去近似f(x)在区间[a,b]上的积分值。最常用的近似函数p(x)是f(x)在区间[a,b]上插值多项式函数。§6.2数值积分一、定义在积分区间上，取个节点作

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

无关。称为插值型求积公式例如Newton-Cotes公式一、Cotes系数取节点为等距分布：由此构造的插值型求积公式称为Newton-Cotes公式,此时求积系数：令Cotes系数二、Newton-Cotes公式1、定义：记则求积公式变为称上式为n阶闭型Newton-Cotes求积公式。注意:由式确定的Cotes系数只与和有关,

与

和积分区间无关，且满足:2、截断误差Newton-Cotes公式的误差为:与x有关三、几种常用的低阶求积公式梯形公式：Simpson公式：

Simpson3/8公式：Cotes公式例3例4、用n=2和n=3的Newton-Cotes公式解：求的近似值。1.n=2时2.n=3时（的真实值为0.7668010）39§6.3求积分的MATLAB命令40格式：S=trapz(x,y)例1：用梯形法求区间内，函数

的定积分值>>x1=[0:pi/30:pi]';>>y=[sin(x1)cos(x1)sin(x1/2)];>>S1=trapz(x1,y)S1=1.99820.00001.9995trapz()函数例2：求解描述被积函数：第一种，一般函数方法functiony=c3ffun(x)y=2/sqrt(pi)*exp(-x.^2);第二种：inline函数f=inline(‘2/sqrt(pi)*exp(-x.^2)’,‘x’);函数定义被积函数：>>y=quad('c3ffun',0,1.5)y=0.9661用inline函数定义被积函数：>>f=inline('2/sqrt(pi)*exp(-x.^2)','x');>>y=quad(f,0,1.5)y=0.9661矩形区域上的二重积分的数值计算格式：矩形区域的双重积分：

y=dblquad(Fun,xm,xM,ym,yM)

限定精度的双重积分：

y=dblquad(Fun,xm,xM,ym,yM,

)例：求解>>f=inline('exp(-x.^2/2).*sin(x.^2+y)','x','y');>>y=dblquad(f,-2,2,-1,1)y=1.57449318974494§6.4

数值微分数值微分是用离散方法近似地计算函数在某点的导数值。在实验和工程应用中，常常得到的离散形式的数据，这时，数值微分显得尤为重要。若在x=a可导，根据导数的定义，可以用差商近似代替微商（导数），有

(1)(2)（h>0且足够小）分别称为向前差商和向后差商。两式平均得(3)称为中心差商。中心差商精度较高。和差商一样，称分别为向前差分、向后差分、中心差分。高阶导数也可用差商法求得，例如二阶导数公式为

(4)2数值微分的MATLAB命令diff(x)向量x的向前差分；%diff(x)/h即为计算数值导数diff(x,n)向量x的n阶向前差分；例1：，求函数的数值导数，并画出的图像。解：f=inline('sqrt(x.^3+2*x.^2-x+12)+(x+5).^(1/6)+5*x+2');x=-3:0.01:3;df=diff(f([x,3.01]))/0.01;plot(x,df)除差商公式外，还有三点公式等其它方法计算数值导数，具体可参考相关教材。例题：教堂顶部曲面面积某个阿拉伯国家有一座著名的伊斯兰教堂，它以中央大厅的金色巨大拱形圆顶名震遐迩。因年久失修，国王下令将教堂顶部重新贴金箔装饰。据档案记载，大厅的顶部形状为半球面，其半径为30m。考虑到可能的损耗和其他技术因素，实际用量将会比教堂顶部面积多1.55%。据此，国王的财政大臣拨出了可制造5800m2有规定厚度金箔的黄金。建筑商人哈桑略通数学，他计算了一下，觉得黄金会有盈余。于是，他以较低的承包价得到了这项装饰工程。但在施工前的测量中，工程师发现教堂顶部实际上并非是一个精确的半球面而是半椭球面，其半立轴恰是30m，而半长轴和半短轴分别是30.6m和29.6m。这一来哈桑犯了愁，他担心黄金是否还有盈余？甚至可能短缺。最后的结果究竟如何呢？解：求解此例，关键在于计算教堂顶部表面积，其中椭球面中心为坐标原点建立直角坐标系，则教堂顶部半椭球面的方程可写为取引入广义极坐标变换不难得到表面积为这个积分式十分复杂，必须用数值方法来处理，由于被积表达式在用上述积分方法，处的奇性，不能直接对此式再作一个变换上式变为：现在可以用数值积分方法。设，即用以下程序计算：x=0:0.25:1;y=0:pi./2:2*pi;A=[12221;24442;24442;24442;12221];w=0;forj=1:5z=30.6*29.6*sqrt(x