数值积分与微分方法

数值积分与微分

摘要

本文首先列举了一些常用的数值求积方法，一是插值型求积公式，以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 。

1.2 Newton Cotes -公式 1.

2.1 Newton Cotes -公式

将积分区间等分，并取分点为求积公式，这样构造出来的插值型求积公式就是Newton Cotes -公式。

()0()()n

b

n k k a

k f x dx b a C f x ==-∑?

（）

其中，

()0

()n

n k k b a b a C =-=-∑

且Cotes 系数满足重要的关系式：

()

1 (k=0,1,2,,n)n

n k

k C

==∑

1n =时，求积公式为梯形公式（两点公式）：

()[()()]2

b a

b a f x dx f a f b -≈+? 梯形公式具有1阶代数精度，余项为：

3

()[]() [,]

12T b a R f f a b ηη-''=-∈

n =2时，求积公式为Simpson 公式(三点公式)：

()[()4()()]62

b a

b a a b f x dx f a f f b -+≈++? Simpson 公式具有3阶代数精度，余项为：

5(4)

1[]()() [,]902

S b a R f f a b ηη-=-∈

n =4时，求积公式为Cotes 公式（五点公式）

： 01234()[7()32()12()32()7()]90

b a

b a f x dx f x f x f x f x f x -≈++++? 其中，

4

k b a

x a k -=+

Cotes 公式具有5次代数精度，余项为：

7(6)

8[]()() [,]9454

C b a R f f a b ηη-=-∈

1.2.2 复合Newton Cotes -公式

当积分区间过大时，直接使用Newton Cotes -公式所得的积分的近似值很难得到保证，因此在实际应用中为了既能够提高结果的精度，又使得算法简便且容易在计算机上实现，往往采用复合求积的方法。

所谓复合求积，就是先将积分区间分成几个小区间，并从每个小区间上用低阶Newton Cotes -公式计算积分的近似值，然后对这些近似值求和，从而得到所求积分的近似值，由此得到一些具有更大实用价值的数值求积公式，统称为复合求积公式。

将[,]a b 区间n 等分，记分点为(0,1,,)k x a kh k n =+= ，其中，b a

h n

-=称为步长，

然后在每个小区间内利用梯形公式，即可导出复合梯形公式：

1

11

1()()[()()]2k k n n

b x k k a x k k h f x dx f x dx f x f x --===≈+∑∑?? 若将所得积分近似值记为n T ，并注意到0,n x a x b ==，则复合梯形公式为：

1

1

()=[()2()()]2n b

n k a

k h

f x dx T f a f x f b -=++∑?

其余项为：

22

10()()()= () [,]1212n n k T k f b a h b a h R f a b n ηηη-=''--''=--∈∑

类似可得复合Simpson 公式：

111012

()[()4()2()()]6n n b n k a k k k h

f x dx S f a f x f x f b --+===+++∑∑? 其中，1/21

2

k k x x h +=+.其余项为：

4(4)

()= ()() [,]1802

n S b a h R f a b ηη--∈

1.2.3 Newton Cotes -公式在MATLAB 中的实现

1）复合梯形数值积分：

调用形式：Z=trapz(X,Y)

其中，X,Y 分别代表数目相同的向量或者数值，Y 与X 的关系可以是函数形态或者离散形态；Z 代表返回的积分值。 2）自适应Simpson 公式

基本调用格式：q=quad （fun ，a ，b ，tol ，trace ，p1，p2） 其中：fun 代表被积函数；a ，b 为积分的上下限； q 为积分结果；tol 为默认误差限，默认了1.e-6；

trace 表示取0表示不用图形显示积分过程，非0表示用图形显示积分过程； p1，p2为直接传递给函数fun 的参数 3）自适应Lobatto 法数值积分：quadl （）

Quadl 是高阶的自适应Newton Cotes -数值积分法函数，比quad 函数更有效，精度更高，使用方法与quad 完全相同。

1.3 Gauss Ledendre -求积公式

1、精度较高Gauss Ledendre -公式

（1）Ledendre 多项式。以Gauss 点k x 为零点的n 次多项式：

12()()()()n n p x x x x x x x =---

上式称为Ledendre 多项式

（2）Gauss Ledendre -求积公式。以Ledendre 多项式的n 个实根为节点的插值求积公式为Gauss Ledendre -求积公式。

考虑在[1,1]-上Gauss 求积公式的构造 1）一个节点

1

1

()2(0)f x dx f -≈?

2）两个节点

二次Ledendre 正交多项式

221

()(31)2

p x x =

- 所以两点的Gauss Ledendre -求积公式为：

1

1

()f x dx f f -≈?

对于一般区间的积分，可以用22

a b b a

x t +-=+将[,]a b 区间转化为[1,1]-，即

1111()()222b a a b b a f x dx f t dt ---+-=+??

然后用相应的Gauss Ledendre -求积公式计算。 （3）一般形式的Gauss Ledendre -求积公式为：

1

()()()

n

b

j j a

j w x f x dx A f x =≈∑?

其中()w x 是一个权重函数，j A 为系数，j x 为横坐标上的节点。

因为()1w x =，所以，一个n 点的求积公式具有如下形式：

1

1

1

()()

n

j j j f x dx A f x -=≈∑?

其中，()j f x 是函数()f x 在节点j x 处的值，节点j x 是Ledendre 正交多项式的根。

222

[()](1)j n

j j A p x x =

'-

MATLAB 没有提供Gauss Ledendre -的有关计算函数，此处给出一部分的编程代码：

function q=gaussL （f ，a ，b ，x ，A ） N=length （x ）； T=zeros （1，N ）；

T=（a+b ）/2+((b-a)/2)*x;

q=((b-a)/2)*sum(A.*feval(‘f’.T));

其中，f 为被积函数；x 和A 的值可有上表查到。

二、数值微分

数值微分的建立常用的三种思路：

1、直接从微分的定义出发，通过近似的处理（泰勒展开），得到数值微分的近似公式；

2、利用插值的基本思想，采用插值近似公式，对插值公式的近似求导得到原数值微分的近似公式

3、根据已知数据，利用最小二乘拟合的方法，得到近似的函数，然后对此近似函数求微分就可以得到数值微分的近似公式。

2.1 差分法近似微分

1、计算公式

在微积分中，一阶微分的计算可以在相邻点x h +和x 间函数取得极限求得。

000()()

()()()()22()lim lim lim h h h h h f x f x f x h f x f x f x h f x h h h

→→→+--+---'=== 所以给出下列差分近似式子： 一阶向前差分：

()()

()f x h f x f x h

+-'=

一阶向后差分：

()()

()f x f x h f x h

--'=

精度较高的一阶中心差分：

()()

22'()h h f x f x f x h

+--=

2、在MATLAB 中的实现 调用形式：Y=diff(X,n)

其中：X 表示求导变量，可以是向量或者矩阵。如是矩阵形式则按照各列做差分；n 表示n 阶差分，即差分n 次；

用diff 函数进行离散数据的近似求导与向前差分近似，但误差较大。可以将数据利用插值或者拟合得到多项式，然后对近似多项式进行微分。

2.2 插值型近似微分

1、方法概述

Lagrange 插值公式，使得 ()()()n n f x L x R x =+

其中，

()()()n

n k k k L x l x f x ==∑，(1)1[()]

()()(+1)n n n f x R x w x n ξ++=！

利用插值公式近似替代原函数，再对插值公式求导，可得插值型求导公式为：

()()()()m m n f x L x =

余项为：

(1)()

1[()]

()[()](+1)k n k n n k d f x R x w x dx n ξ++=！

特别的，n=1时，可得一阶微分两点公式为：

()()

()f x h f x f x h +-'=

()()

()f x f x h f x h --'=

n=2时，

'0200121

'()()[3()4()()]2f x L x f x f x f x h

≈=-+-

'121021

'()()[()()]2f x L x f x f x h ≈=-+

'2220121

'()()[()4()3()]2f x L x f x f x f x h

≈=-+

下面给出一个常用的五点公式：

021121

'()[()8()8()()]

12f x f x f x f x f x h --≈-+-

2、三次样条插值函数求微分的MATLAB 函数

由于三次样条插值的导数近似被插值函数导数的效果很好，此处给出三次样条插值函数的MATLAB 调用步骤：

Step1：对离散数据用csapi 函数（或者spline 函数），得到其三次样条插值函数

调用形式pp=csapi （x ，y ）

其中，x ，y 分别为离散数据对的自变量和因变量；pp 为得到的三次样条插值函数 Step2：用fnder 函数求三次样条插值函数的导数

调用形式 fprime=fnder （f ，dorder ）

其中，f 为三次样条插值函数，dorder 为三次样条插值函数的求导阶数； fprime 为得到的三次样条插值函数的导数值 Step3：用fnval 函数求导函数在未知点处的导数值

调用形式 v=fnval （fprime ，x ）

其中，fprime 为三次样条插值函数导函数；x 为未知点处自变量值；v 为未知点处的导数值。

三、数值积分与微分的比较

1、数值解法

微积分是高等数学的重要内容，在实际工程中有许多重要的应用。微积分的数值解法，是不同于高等数学中的解析方法，适合求解没有或者很那求出微分或者积分解析表达式的实际问题的计算。

2、数值积分与微分与插值和拟合的关系 数值微分与数值积分依赖插值和拟合，二者之间密不可分。比如在进行数值微分时，针对离散的数据点，常常利用插值和拟合来减少数据误差。数值积分的基本思路也来自于插值法。比如当所积函数的形式比较复杂或是通过表格形式给出，则可以通过构造插值多项式来代替原函数，简化问题。

插值型求积公式是以构造插值函数代替原函数进行积分：

(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 ξ+==++

+?

?∑??（）

插值型微分公式是利用插值函数代替原函数进行求导：

()()()()m m n f x L x =

3、精度比较

插值积分的余项为：

(1)

+11[][()](1)!b n n n a R f f x w x dx n ξ+=+?（）

插值微分的精度为：

(1)1[()]

()()(+1)n n n f x R x w x n ξ++=！

可见二者的精度形式完全一样。因为都是根据插值基函数演变而来，都是基于泰勒展开式展开的。

四、参考文献

[1]杜廷松，覃太贵，《数值分析及实验》科学出版社 2012.10

[2]数值微分与数值积分https://www.doczj.com/doc/5910834251.html,/p-620586060.html 2015.12.7

数值分析第四章数值积分与数值微分习题复习资料

第四章 数值积分与数值微分 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(-++=，则可得如下递推公式

Matlab数值积分与数值微分

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高斯-克朗罗德法

数值微分与数值积分

专题六数值微积分与方程求解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 --==?+?+???==?+?+?? ???

数值微分与数值积分练习题

第五章 数值微分与数值积分 一．分别用向前差商，向后差商和中心差商公式计算()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

实验4_数值积分与数值微分

数值分析实验报告四 数值积分与数值微分实验（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 #include /*被积函数的定义*/ float f(float x) {

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 #include /*被积函数的定义*/ float f(float x) { float y; y=sqrt(x); return y; } void main()

数值积分与微分方法

数值积分与微分 摘要 本文首先列举了一些常用的数值求积方法，一是插值型求积公式，以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章 数值积分与数值微分

第七章数值积分与数值微分 积分问题最早来自于几何形体的面积、体积计算，也是经典力学中的重要问题（例如计算物体的重心位置）. 在现实应用中，很多积分的结果并不能写成解析表达式，因此需要通过数值方法来计算. 数值微分是利用一些离散点上的函数值近似计算某一点处的函数导数，它针对表达式未知的函数. 本章介绍一元函数积分（一重积分）和微分的各种数值算法，它们也是数值求解积分方程、微分方程的基础. 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）

数值积分与微分方程

2.3 数值积分 2.3.1 一元函数的数值积分 函数1 quad 、quadl 、quad8 功能 数值定积分，自适应Simpleson 积分法。 格式 q = quad(fun,a,b) %近似地从a 到b 计算函数fun 的数值积分，误差为10-6。 若给fun 输入向量x ，应返回向量y ，即fun 是一单值函数。 q = quad(fun,a,b,tol) %用指定的绝对误差tol 代替缺省误差。tol 越大，函数计 算的次数越少，速度越快，但结果精度变小。 q = quad(fun,a,b,tol,trace,p1,p2,…) %将可选参数p1,p2,…等传递给函数 fun(x,p1,p2,…)，再作数值积分。若tol=[]或trace=[]，则用缺省值进行计算。 [q,n] = quad(fun,a,b,…) %同时返回函数计算的次数n … = quadl(fun,a,b,…) %用高精度进行计算，效率可能比quad 更好。 … = quad8(fun,a,b,…) %该命令是将废弃的命令，用quadl 代替。 例2-40 >>fun = inline(‘3*x.^2./(x.^3-2*x.^2+3)’); equivalent to: function y=funn(x) y=3*x.^2./(x.^3-2*x.^2+3); >>Q1 = quad(fun,0,2) >>Q2 = quadl(fun,0,2) 计算结果为： Q1 = 3.7224 Q2 = 3.7224 补充：复化simpson 积分法程序 程序名称 Simpson.m 调用格式 I=Simpson('f_name',a,b,n) 程序功能 用复化Simpson 公式求定积分值 输入变量 f_name 为用户自己编写给定函数()y f x 的M 函数而命名的程序文件名 a 为积分下限 b 为积分上限 n 为积分区间[,]a b 划分成小区间的等份数 输出变量 I 为定积分值 程序 function I=simpson(f_name,a,b,n) h=(b-a)/n; x=a+(0:n)*h; f=feval(f_name,x);

数值积分与数值微分 编程计算

解：卫星轨道的示意图如右上图所示，,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 是由测量或数值计算给出的数据表，是列表函数，也无法用牛顿—莱布尼兹公式计算. 对于上述这些情况，都要求建立定积分的近似计算方法——数值积分法.