第七章 基于MATLAB的科学计算—数值微积分1
第7章 MATLAB数值积分与数值微分
MATLAB数值积分 MATLAB数值微分
7.1 数值积分
7.1.1 数值积分基本原理
关于积分，有Newton-Leibniz公式
f ( x)dx F (b) F (a) a 但是，在很多情况下，还是要求数值积分： 1、函数由离散数据组成 2、F(x)求不出 3、F(x)非常复杂 定义数值积分如下：离散点上的函数值的线性组合
In ( f )
b
a
i 0
n
i
f ( xi )
称为积分系数，与f(x)无关，与积分区间和积分点有关
插值型 用插值函数的积分，作为数值积分
In ( f )
b b n n a
ai
b l ( x)dx f ( x ) Ln ( x)dx li (x) f ( xi )dx i i a a i 0 i 0
微积分中，关于导数的定义如下：
f ( x h) f ( x ) f ( x ) f ( x h) f ( x h) f ( x h) f ' ( x ) lim lim lim h 0 h 0 h 0 h h 2h
自然而又简单的方法就是取极限的近似值，即差商
标系中做出f'(x)的图像。程序如下：
x=-3:0.01:3; p=polyfit(x,f(x),5); dp=polyder(p); dpx=polyval(dp,x);
%用5次多项式p拟合f(x) %对拟合多项式p求导数dp %求dp在假设点的函数值
f=inline('sqrt(x.^3+2*x.^2-x+12)+(x+5).^(1/6)+5*x+2'); dx=diff(f([x,3.01]))/0.01; %直接对f(x)求数值导数
MATLAB_简介(5)MATLAB数值积分与微分(1)
>> plot(x,y,'o',x,y)
>> title('y(x) data plot') >> ylabel('y(x)'), xlabel('x') >> dy=diff(y)./diff(x); >> xd=x(1:length(x)-1); >> plot(xd,dy) >> title('Approximate derivative using diff') >> ylabel('dy/dx'), xlabel('x')
%函数计算的次数n
q= -0.4605 n= 53
再来看一个积分式
L
2
0
cos(t ) 4 sin(2t ) 1dt
2 2
Function f=len(t) f=sqrt(cos(t),^2+4*sin(2t),^2+1);
存为Len.m
>>Lenth=quad(‘Len’,0,2*pi)
以下的例子是针对数据组为离散型态，要注意 的是原数据所代表的函数分布并无明显的折角，
但是它的 一阶微分经数值微分计算后的微分函
数分布却有极大的曲折变化。 >> x=0:0.1:1;
>> y=[-.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56
9.48 9.30 11.2];
>> deno=x(3:length(f))-x(1:length(f)-2); % 注意中
央差分是 x(k+1)-x(k-1)
数值积分法matlab
数值积分法matlab数值积分法是一种通过数值近似来计算定积分的方法。
在实际问题中，很多函数的积分无法用闭合形式表达出来，这时就需要使用数值积分法来近似求解。
数值积分法的基本思想是将要积分的区间分割成若干小区间，然后在每个小区间上用一条简单的函数来逼近原函数，最后将这些小区间上的近似积分结果相加。
常用的数值积分法有矩形法、梯形法和辛普森法等。
其中，矩形法是最简单的数值积分法之一。
它将每个小区间上的函数值看作是该区间上函数的常值近似，并用矩形面积来表示该区间上的积分值。
矩形法有两种类型，即左矩形法和右矩形法。
左矩形法使用每个小区间左端点处的函数值来代表该区间上的函数值，右矩形法则使用每个小区间右端点处的函数值。
通过将所有小区间上的矩形面积相加，即可得到对整个区间上函数积分的近似值。
梯形法是数值积分法中更精确的一种方法。
它通过在每个小区间上使用梯形面积来逼近函数的积分值。
梯形法的基本思想是将每个小区间上的函数近似表示为两个端点处函数值的线性插值函数。
通过计算每个小区间上的梯形面积，并将这些面积相加，即可得到对整个区间上函数积分的近似值。
辛普森法是数值积分法中最常用的一种方法，它通过在每个小区间上使用二次多项式来逼近函数的积分值。
辛普森法的基本思想是将每个小区间上的函数近似表示为一个二次多项式，并计算该多项式对应的曲线下面积。
通过将所有小区间上的曲线下面积相加，并乘以一个系数，即可得到对整个区间上函数积分的近似值。
在使用数值积分法时，需要注意选择合适的分割数和逼近方法，以获得更精确的结果。
通常情况下，分割数越多，逼近结果越接近真实值。
但是，分割数过大也会增加计算量。
因此，需要在计算精度和计算效率之间进行权衡。
除了上述介绍的几种数值积分法外，还有其他一些方法，如高斯积分法和自适应积分法等。
这些方法在不同的情况下有着不同的适用性和计算效果。
因此，在实际问题中，需要根据具体情况选择合适的数值积分方法。
总结而言，数值积分法是一种通过数值近似来计算定积分的方法。
MATLAB数值微积分
4.1数值微积分4.1.1近似数值极限及导数Matlab 数值计算中，没有求极限指令，也没有求导指令，而是利用差分指令：用一个简单矩阵表现diff和gradient指令计算方式。
差分： Dx=diff(X)对向量： Dx=X(2:n)-X(1:n-1)对矩阵： DX=X(2:n,:)-X(1:n-1,:)长度小1.DIFF(X), for a vector X, is[X(2)-X(1) X(3)-X(2) ... X(n)-X(n-1)]. DIFF(X), for a matrix X, is the matrix of row differences, （结果缺少一行）[X(2:n,:) - X(1:n-1,:)].DIFF(X,N,DIM) is the Nth difference function along dimension DIM.If N >= size(X,DIM), DIFF returns an empty array （N阶差分）梯度：FX=gradient(F)Fx(1)=Fx(2)-Fx(1);F=[1,2,3;4,5,6;7,8,9]Dx=diff(F) （按行）Dx_2=diff(F,1,2) （按列）[FX,FY]=gradient(F)Fx(1)=Fx(2)-Fx(1),Fx(end)=F(end)-F(end-1)FX与F维数相同。
[FX_2,FY_2]=gradient(F,0.5)%采样间隔0.5 即：Fx(1)=(Fx(2)-Fx(1))/2F =1 2 34 5 67 8 9Dx =3 3 33 3 3Dx_2 =1 11 11 1FX =1 1 11 1 1 1 1 1 （列方向） FY =3 3 3 3 3 3 3 3 3 （行方向） FX_2 =2 2 2 2 2 2 2 2 2 FY_2 =6 6 6 6 6 6 6 6 6【例4.1-1】设xx xx f sin 2cos 1)(1-=，xxx fsin )(2=，试用机器零阈值eps 替代理论0计算极限)(lim )0(11x f L x →=，)(lim )0(22x fL x →=。
第7章 MATLAB数值积分与数值微分PPT教学课件
MATLAB数值积分 MATLAB数值微分
2020/12/10
1
7.1 数值积分
7.1.1 数值积分基本原理 求解定积分的数值方法多种多样，如简单方法：
梯形法、辛普生(Simpson)法、牛顿－柯特斯 (Newton-Cotes)法等方法。它们的基本思想都是 将整个积分区间[a,b]分成n个子区间[xi,xi+1]， i=1,2,…,n，其中x1=a，xn+1=b。求定积分问题就分 解为求和问题。
其中，fname是被积函数名。a和b分别是定积分的下 限和上限。tol用来控制积分精度，默认时取tol=10-6。 trace控制是否展现积分过程，若取非0则展现积分过 程，取0则不展现，默认时取trace=0。返回参数I即 定积分值，n为被积函数的调用次数。
2020/12/10
3
例7-1 求定积分。
2020/12/10
10
例7-4 用trapz函数计算定积分。 命令如下：
X=1:0.01:2.5; Y=exp(-X); %生成函数关系数据向量 trapz(X,Y) ans =
0.2858
2020/12/10
Байду номын сангаас11
7.1.3 多重定积分的数值求解
MATLAB提供的dblquad函数用于求二重积分的数值 解，triplequad函数用于求三重积分的数值解。函 数的调用格式为
2020/12/10
2
7.1.2 数值积分的实现方法
1．变步长辛普生法
基于变步长辛普生法，MATLAB给出了quad函数和 quadl函数来求定积分。函数的调用格式为：
[I,n]=quad(@fname,a,b,tol,trace)
第七章 基于MATLAB的科学计算—数值微积分1
科学计算—理论、方法及其基于MATLAB 的程序实现与分析数 值 微 积 分§1 数值微分对于给定的函数()x f y =，如果 １、()x f y =的函数关系式比较复杂时；２、()x f y =未知，而仅仅知道()x f y =在1+n 个相异点k x ，n k ,,1,0 =处 的函数值k y ；则希望能用相对简单的计算方法，求得()x f y =导数的（近似）值。
基于上述考虑，选择的方法之一是利用函数()x f y =的插值多项式的导数作为函数()x f y =导数的近似值，例如Lagrange 插值多项式，由于()()()()∑==≈nk k k n x f x l x L x f 0 （１）因而有()()x L x f n'≈' （２） 这里需要说明一点的是，尽管()x f 和()x L n 的函数值可能相差不多，但是它们的导数有可能相差很大，如下面的例子例1： 考虑函数()22511x x f +=在区间[-1，1]的插值问题，取区间[-1，1]的11个点k x k ⨯+-=2.01，10,,1,0 =k ，作函数()22511x x f +=的10次插值多项式：()18552.163597.1234338.3819095.4949417.220246810+-+-+-=x x x x x x L n函数()x f 和插值多项式()x L n 的导数分别为()()2225150xxdx x df +-=()x x x x x dxx dLn 7.334.4936.22883.39594.22093579-+-+-= 对函数()x f 和插值多项式()x L n 及其导数分别比较，结果如图所示：Derivative_Runge下面我们基于Taylor 公式()()()()()()()()()()⎥⎦⎤⎢⎣⎡≅++++''+'+=++++1112!1!!2n n n n n h O h n f h n x f h x f h x f x f h x f ξ (3) 讨论导数的近似计算问题 1 一阶导数的近似计算令k k x x h -=+1,可得一阶向前有限差商公式(First Forward Finite Divided Difference ):()()()()()()()()()()()()()()⎪⎩⎪⎨⎧=-≈'⇒+''--='⇒+''+'+=+++h O error h x f x f x f h O h x f h x f x f x f h O h x f h x f x f x f k k k k k k k k k k k 121321!2!2 (4)类似地,当1--=k k x x h 时,可得一阶向后有限差商公式(First BackwardFinite Divided Difference ):()()()()()()()()()()()()()()⎪⎩⎪⎨⎧=-≈'⇒+''+-='⇒+''+'-=---h O error h x f x f x f h O h x f h x f x f x f h O h x f h x f x f x f k k k k k k k k k k k 121321!2!2 (5)由近似计算公式(4)和(5),可得可得一阶中心有限差商公式(FirstCentered Finite Divided Difference ):()()()()()()()()()()⎪⎩⎪⎨⎧=-≈'⇒+-='⇒+-+-+2112112254h O error hx f x f x f h O hx f x f x f k k k k k k (6)2 二阶导数的近似计算在Taylor 公式中，用h 2代替h 可得()()()()()()3222!22h O h x f h x f x f x f k k k k +''+'+=+ (7) 式(7)与式(4)结合可得二阶向前有限差商公式(Second Forward Finite Divided Difference ):()()()()()32122h O h x f x f x f x f k k k k ++''+-=-++ (8)()()()()()⎪⎩⎪⎨⎧=+-≈''⇒++h O error h x f x f x f x f k k k k 2122 ()()()()()()⎪⎩⎪⎨⎧=---≈''⇔+++h O error h h x f x f h x f x f x f k k k k k 112 (9) 类似地利用()()()()()()3222!22h O h x f h x f x f x f k k k k +''+'-=- ()()()()()321!2h O h x f h x f x f x f k k k k +''+'-=-可得二阶向后有限差商公式(Second Backward Finite Divided Difference ):()()()()()32122h O h x f x f x f x f k k k k ++''+-=---()()()()()⎪⎩⎪⎨⎧=+-≈''⇒-+h O error h x f x f x f x f k k k k 2212 (10) 进一步底，由公式()()()()()()()()()()()()43214321!3!2!3!2h O h x f h x f h x f x f x f h O h x f h x f h x f x f x f k k k k k k k k k k +'''-''+'-=+'''+''+'+=-+ (11)可得二阶中心有限差商公式(Second Centered Finite DividedDifference ):()()()()()⎪⎩⎪⎨⎧=+-≈''-+22112hO error h x f x f x f x f k k k k (12)§2 数值积分§2.1 Newton-Cotes 数值积分公式基于数值微分的同样原因，当函数()x f y =在],[b a 上的定积分不易计算时，利用函数()x f y =的Lagrange 插值多项式()()()∑==nk k k n x f x l x L 0，有()()()()k nk b a k ban b ax f dx x l dx x L dx x f I ∑⎰⎰⎰=⎥⎦⎤⎢⎣⎡=≈=0 (13) 由于Lagrange 插值多项式()()()∑==nk k k n x f x l x L 0的插值基函数()x l k 只依赖插值节点k x ，n k ,,1,0 =，所以当k x ，n k ,,1,0 =取定后，()x l k 就是完全确定的多项式函数，令()⎰=bak k dx x l A ，n k ,,1,0 = (14)则由式(13)得到Newton-Cotes 求积公式：()()()k nk k ban ba x f A dx x L dx x f I ∑⎰⎰==≈=0(15)特别地，当取插值节点为kh a x k +=，n k ,,1,0 =，nab h -=(16) 时有１） 两点公式（Trapezoidal Rule ）：⇒=1n 2210a b h A A -=== ()()()()()[]b f a f ab x f A dx x L dx x f I k n k k b a ba+-==≈=∑⎰⎰=201 (17) 2） 三点公式（Simpson's 1/3 Rule ）：34,322120hA h A A a b h n ===⇒-=⇒=， ()()()()()()⎥⎦⎤⎢⎣⎡+⎪⎭⎫ ⎝⎛++=⎥⎦⎤⎢⎣⎡+⎪⎭⎫ ⎝⎛++-=≈=⎰⎰b f b a f a f h b f b a f a f a b dx x L dx x f I b ab a2432462 (18) 3） 四点公式（Simpson's 3/8 Rule ）：89,83332130hA A h A A a b h n ====⇒-=⇒= ()()()()()()[]()()()()[]b f h a f h a f a f hb f h a f h a f a f ab dxx L dx x f I b ab a+++++=+++++-=≈=⎰⎰2338323383 (19) 4） 误差估计：利用Lagrange 插值的误差公式：()()()()()()()∏=+-+=-=nk k n n n x x n f x L x f x R 01!1ξ (20) 容易得到Trapezoidal Rule 的误差估计：()()()()()()()()()()()32321121222a b M a b f dx ab x b a x f dx b x a x f dxx L dx x f b a ba bab a-≤-''=++-''≈--''=-⎰⎰⎰⎰ξξξ (21) Simpson's 1/3 Rule 以及Simpson's 3/8 Rule 的误差估计分别为：()()()()5445229012901⎪⎭⎫⎝⎛-≤⎪⎭⎫ ⎝⎛--≈-⎰⎰a b M f a b dx x L dx x f bab a ξ (22)()()()()5445338033803⎪⎭⎫⎝⎛-≤⎪⎭⎫ ⎝⎛--≈-⎰⎰a b M f a b dx x L dx x f bab aξ (23) 其中()(){}x f M k b a k ],[max =.注１ 数值积分公式(13)表明，数值积分方法最大的优点是将复杂的函数积分转化为相对简单的加、减、乘、除运算；注２ 从几何的角度看，上述三种求积公式（17）、（18）和（19）（在积分区间],[b a 上）分别用直线()()b f ab ax a f b a b x y --+--=（24） 抛物线()()()()()()()()()b f h h a x a x h a f h b x a x a f h h a x h a x y 2222---++-------=（25）三次立方曲线()()()()()()()()()()()()()()()()h a f h h a x h a x a x h a f h h a x h a x a x h a f h h a x h a x a x a f h h a x h a x h a x y 332332323333+-----++------+-----+-------= （26）代替曲线()x f y =:例1 计算积分⎰=10dx e I x ，并估计误差. 解 1） 用梯形公式计算，8591409.1)(20110≈+-=≈e e T I .由于x e x f =)(，所以xe xf ='')(，于是，梯形公式的误差23.012)(max 12110≈=''≤-=≤≤ex f T I R x T .2） 用辛普森公式计算，7188612.1)4(60112/10≈++-=≈e e e S I . 由于x ）（e xf =)(4，于是，辛普森公式的误差00094.02880)(max 28801410≈=≤-=≤≤ex f S I R ）（x S .【注】 （1） 当0)()2(=x f 时，梯形求积公式准确成立；当)()4(x f =0时，辛普森公式准确成立.即梯形公式对一次多项式准确成立，而辛普森公式对三次多项式准确成立.（2） 一般地，由余项公式(7.9)知，当)(x f 是n 次多项式时，积分余项为零，从而牛顿-柯特斯求积公式准确成立.（3） 数值求积方法是一种近似方法，因此，要求求积公式对尽可能多的被积函数)(x f 能准确计算积分值.作为衡量公式逼近好坏的标准之一，下面给出代数精度(the Degree of Precision of the Quadrature Formula)的概念.§2.1.2 代数精度【定义1】 如果某求积公式对于次数小于等于m 的多项式能准确求出积分值，而对某个1+m 次多项式就不能准确求出积分，则称该求积公式具有m 次代数精度.由定义1，具有m 次代数精度的求积公式（7.3）对mx x x f ,,,1)( =时精确成立，而对1)(+=m x x f 不能精确成立.为了构造形如公式（7.3）的求积公式，通过解方程组⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧+-==-==-==∑⎰∑⎰∑⎰=++==n k b a m m mm k k n k b a k k nk b a k m a b dx x x A a b xdx x A a b dx A 0110220,1,2,1可得求积节点k x 与求积系数k A ，由此求得具有m 次代数精度的求积公式.不难验证，梯形公式和矩形公式均具有一次代数精度，而辛普森公式具有三次代数精度.对牛顿-柯特斯求积公式，有下面结论.【定理1】 牛顿-柯特斯求积公式（7.4）至少具有n 次代数精度，而当n 为偶数时，至少具有1+n 次代数精度.证明 当)(x f 为任何次数不高于n 的多项式时，0)()1(≡+x f n ，所以0)(=⎰ban dx x R ,显然结论成立.当n 为偶数时，只须对1)(+=n x x f 时的结论验证.因为)!1()()1(+=+n x f n ，由截断误差公式⎰⎰⎰∏∏=+=-=-=bab an nj n nj jn dtj t hdx xx dx x R 002)()()(，若令2n u t +=，则有⎰⎰∏-=+-+=ban n n j n n du j nu hdx x R 2202)2()(.注意到∏=-+nj j n u 0)2(是u 的奇函数（n 为偶数），因此⎰=b a n dx x R 0)(，结论成立.例2 对于求积公式20011()(1)(2)'(2)f x dx A f B f B f -=-++⎰，试确定100,,B B A ，使此求积公式有最高的代数精度，是几阶的？§2.2 复化求积公式由前面分析知道，在整个积分区间],[b a 上用直线或用抛物线代替原曲线()x f y =做积分，尽管计算简单，但精度较差，结果不好，因此，一个自然的想法是：在尽可能保持公式(17)、和(18)计算简单、易于计算机实现的优点的前提下，寻找提高误差精度的途径和方法。
matlab数值微积分
人工智能与机器学习
结合人工智能和机器学习技术，自动选择合适的 数值微积分算法，实现自适应计算。
ABCD
并行计算
利用多核处理器和分布式计算资源，实现数值微 积分的并行计算，加速大规模问题的求解。
云平台集成
将Matlab数值微积分与云平台集成，实现数据 共享、远程计算和动态资源调度。
未来展望
01
更广泛的应用领域
工程领域
Matlab广泛应用于数学、物理、化学、生 物等领域的数值计算和数据分析。
Matlab在机械、电子、控制、航空航天等 工程领域有广泛应用，支持各种工程设计 和仿真。
金融领域
图像处理和计算机视觉
Matlab在金融领域主要用于数据分析、统 计建模、风险评估等方面。
Matlab提供了图像处理工具箱和计算机视 觉工具箱，广泛应用于图像处理和计算机 视觉领域。
quad: 这是一个用于数值积分 的函数，可以计算一维函数的 定积分。与integral函数不同 的是，它使用自适应Simpson 方法进行积分。例如，对于函 数f(x)，可以使用以下代码计 算定积分
数值积分函数
```matlab
result = quad(@(x) f(x), a, b);
数值积分函数
互操作性和兼容性。
THANKS
感谢观看
将连续的问题离散化，用 有限个点来近似表示连续 的函数。
逼近
通过选取适当的离散点， 使用数学方法逼近真实的 函数值。
迭代
通过不断迭代逼近真实值， 提高计算的精度。
数值微积分的计算方法
差分法
01
通过差分代替导数，将微分问题转化为差分问题，进而求解。
辛普森法则
02
利用区间中点的函数值和区间端点的函数值来近似计算定积分。
第七章 MATLAB微积分数值计算
相同维数的向量U，且 1 U1 2 ( V4 4V3 5V2 2V1 ) h 1 U i 2 (Vi 1 2Vi Vi 1 ), (1 i n) h 1 U n 2 (2Vn 5Vn 1 4Vn 2 Vn 3 ) h 默认的步长为1。 U=4*del2(V,h1,h2)，对矩阵V，横向(x方向)以步长h1，纵向(y方
diff调用格式为：
Dy=diff(Y)：计算向量Y的向前差分，并把结果赋值给向量Dy Dy(i)=Y(i+1)-Y(i)，i=1,2,…,n-1。注意向量Dy元素个数比Y少一个
Dy=diff(Y,n)：计算向量Y的n阶向前差分。注意向量Dy元素个数比Y少n个.例如：
diff(Y,2)=diff(diff(Y))=DX(i+1)-DX(i)= Y(i+2)-2Y(i+1)+Y(i) ， i=1,2 ……n-2。
计算积分，可以采取逐步缩小步长h的办法。即先任
取步长h进行计算，然后取较小步长 h’ 进行计算，如果两
次计算结果相差较大，则取更小步长进行计算，如此下去， 直到相邻两次计算结果相差不大为止，取最小步长算出的 结果作为积分值。这种方法称为变步长积分法。 利用两种步长计算积分时，通常取h’=h/2 。而每次
x
4
精度为O(∆X4)的高阶中心差分算法
yi yi yi yi 2 8 yi 1 8 yi 1 yi 2 12 x yi 2 16 yi 1 30 yi 16 yi 1 yi 2 12 x
2
yi 3 8 yi 2 13 yi 1 13 yi 1 8 yi 2 yi 3 8 x
7.1 数值微分 7.2 数值积分 7.3 常微分方程的数值解法
第7章 MATLAB数值微分与积分_习题答案
1第7章 MATLAB 数值微分与积分习题7一、选择题1．diff([10,15])的值是（ ）。
AA ．5B ．10C ．15D ．252．数值积分方法是基于（ ）的事实。
DA ．求原函数很困难B ．原函数无法用初等函数表示C ．无法知道被积函数的精确表达式D ．A ，B ，C 三个选项3．求数值积分时，被积函数的定义可以采取（ ）。
DA ．函数文件B ．内联函数C ．匿名函数D ．A ，B ，C 三个选项4．以下选项不能用来求数值积分的函数是（ ）。
BA ．quadgkB ．quad2C ．integralD ．integral25．以下选项不是离散傅里叶变换的函数是（ ）。
CA ．fftB ．fft2C ．fft1D ．fftn二、填空题1．在MATLAB 中，没有直接提供求 的函数，只有计算 的函数diff 。
数值导数，向前差分2．基于变步长辛普森法，MATLAB 给出了 函数和 函数来求定积分。
quad ，quadl3．MA TLAB 提供了基于全局自适应积分算法的 函数来求定积分，该函数的积分限 （可以或不可以）为无穷大。
integral ，可以4．MATLAB 提供的 、 、 函数用于求二重积分的数值解， 、 函数用于求三重积分的数值解。
integral2，quad2d ，dblquad ，integral3，triplequad5．MA TLAB 提供了离散傅里叶变换函数fft ，对应的逆变换函数是 。
ifft三、应用题1．求函数在指定点的数值导数。
（1）2346x x x x f 22ππππ,,,,cos sin)(=+= （2）321x 1x x f 2,,,)(=+=2(1):(2):直接用导数函数求:f=inline('x./sqrt(x.^2+1)');f(1)用拟合函数求:f=inline('sqrt(x.^2+1)');x=0:0.001:5;p=polyfit(x,f(x),5);dp=polyder(p);dpx=polyval(dp,1)2．求定积分。
Matlab基础及其应用 第7章 数值微积分与常微分
6
6
6
1 −19
6 −19
6
7.2 数值积分
MATLAB基础与应用教程
7.2.1 数值积分的原理
基本思想是将整个积分区间[a，b]分成n个子区间[xi，x
i + 1]，i = 1，2，…，n，其中x1 = a，xn + 1 = b。这样求
定积分问题就分解为下面的求和问b 题：
S a f (x)dx
MATLAB基础与应用教程
7.2 数值积分
MATLAB基础与应用教程
7.2.2 定积分的数值求解实现
1．自适应积分算法 基于全局自适应积分算法的integral函数来求定积分。函数 的调用格式为 q = integral(@fun,xmin,xmax) q = integral(@fun,xmin,xmax,Name,Value) fun是被积函数，xmin和xmax分别是定积分的下限和上限。 选项用于设置积分器的属性
n
xi1
i 1 xi
f
(x)dx
矩形法是用矩形面积近似曲边梯形的面积，如图7.2（a）所 示；梯形法是用斜边梯形面积近似曲边梯形的面积，如图7. 2（b）所示；而辛普生法是用抛物线近似曲边。
7.2 数值积分
7.2.1 数值积分的原理
1．梯形法 2．自适应辛普生法 3．高斯-克朗罗德（Gauss-KrLeabharlann nrod）法7.2 数值积分
【例7.5】求定积分
1ex ln xdx 0
>> format long; >> I=quadgk(@(x)exp(x).*log(x),0,1) I=
−1.317902162414081
MATLAB基础与应用教程
MATLAB 数值微积分与微分方程式求解
Matlab求解常微分方程式
Ex:
dy1 dt
=
y1
+
y2e−t
dy2 dt
=
− y1 y2
+ cos(t)
1. edit fun.m
function dydt=fun(t,y)
dydt(1) = y(1)+y(2)*exp(-t);
dydt(2) = -y(1)*y(2)+cos(t);
Stiff ODE 指的是其內部某些狀態響應快速，而某些則 相對具較緩慢動態
勁度
勁度系統 (stiff system)表示其具有快速變化以及緩慢變化 的部分。 勁度系統例子： dy = −1000y + 3000 − 2000e−t
dt
假使y(0)=0，其解 y = 3 − 0.998e−1000t − 2.002e−t
用前述方法求解用前述方法求解用前述方法求解用前述方法求解0010ysin??yettt21221sincosttyyyyyteeyty?????12122210100sincosttyy?yyyteeyty???????11matlabode指令matlab用于求解起始值常微分方程式問題的指令stiffode指的是其內部某些狀態響應快速而某些則相對具較緩慢動態ode23tbode23tode23sode15sstiffodeode113ode23ode45nonstiffode指令問題形式勁度勁度系統勁度系統勁度系統勁度系統stiffsystem表示其具有快速變化以及緩慢變化的部分
k4 = f (ti + h, yi + k3h)
方程式系統
許多實際的工程及科學問題需要求解的是聯立常微分方 程式系統，而不只是單一方程式。
第7章MATLAB微积分运算
例 7-6
编制如下程序
clear;clc; [X, Y, Z] = sphere(20); surfnorm(X, Y, Z) axis square box on
运行结果如图7-5所示
图 7-5
7
7.1.2 符号微分
例 7-7
编制如下命令文件
clear;clc; F='sin(3*x)' Fx=diff(F,'x') Fxx=diff(F, 'x', 2)
运行结果如图7-1所示 图 7-1
3
例 7-3
编制如下程序
clear;clc; n = 360 x = linspace(0, 2*pi, n); y = exp(-0.2*x).*(cos(x)+sin(x)); dy = diff(y); dx = diff(x); dyx = dy./dx; plot(x,y,'r-') hold on plot(x(1:n-1),dyx,'b-.') legend('y','dy/dx')
运算结果为
F= sin(3*x) Fx =
3*cos(3*x)
Fxx =
-9*sin(3*x)
8
例 7-8
编制如下程序
clear;clc; Z = 'exp(a*x+b*y)*(sin(x*y)+sin(x)+cos(y))' Zx = diff(Z, 'x') Zy = diff(Z, 'y') 运行结果为
cos(a*x)+b
F=
(sin(pi*a)+b*a*pi)/a
Matlab数值积分与数值微分
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+π)dx3π的值.先建⽴函数⽂件fesin.mfunction f=fesin(x)f=exp(-0.5*x).*sin(x+(pi/6));再调⽤quad函数[I,n]=quad(@fesin,0,3*pi,1e-10)I=0.9008n=365--------------------------------------------------------------------- 例如：分别⽤quad函数和quadl函数求积分∫e0.5x sin(x+π6)dx3π的近似值，⽐较函数调⽤的次数.先建⽴函数⽂件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⾼斯-克朗罗德法Matlab提供了⾃适应⾼斯-克朗罗德法的quadgk函数来求震荡函数的定积分，函数的调⽤格式为：[I,err]=quadgk(@fname,a,b)Err返回近似误差范围，其他参数的意义与quad函数相同，积分上下限可以是-Inf或Inf，也可以是复数，若为复数则在复平⾯上求积分.--------------------------------------------------------------------- 例如：求积分∫xsinx1+cos2xdx π的数值.先编写被积函数的m⽂件fsx.mfunction f=fsx(x)f=x.*sin(x)./(1+cos(x).^2);再调⽤quadgk函数I=quadgk(@fsx,0,pi)I=2.4674--------------------------------------------------------------------- 例如：求积分∫xsinxdx +∞∞的值.先编写被积函数的m⽂件fsx.mfunction f=fsx(x)f=x.*sin(x)./(1+cos(x).^2); 再调⽤quadgk函数I=quadgk(@fsx,-Inf,Inf)I=-9.0671e+017---------------------------------------------------------------------1.3梯形积分法对于⼀些不知道函数关系的函数问题，只有实验测得的⼀组组样本点和样本值，由表格定义的函数关系求定积分问题⽤梯形积分法，其函数是trapz函数，调⽤格式为：I=Traps(X,Y)X,Y为等长的两组向量，对应着函数关系Y=f(X) X=(x1,x2,…,x n)(x1分区间是[x1,x n]--------------------------------------------------------------------- 例如：已知某次物理实验测得如下表所⽰的两组样本点.现已知变量x和变量y满⾜⼀定的函数关系，但此关系未知，设y=f(x)，求积分13.39∫f(x)dx1.38的数值.X=[1.38,1.56,2.21,3.97,5.51,7.79,9.19,11.12,13.39];Y=[3.35,3.96,5.12,8.98,11.46,17.63,24.41,29.83,32.21]; I=trapz(X,Y) I=217.1033---------------------------------------------------------------------例如：⽤梯形积分法求积分：∫e ?x dx 2.51的数值.x=1:0.01:2.5; y=exp(-x); I=trapz(x,y) I= 0.2858---------------------------------------------------------------------2. 多重数值积分的实现重积分的积分函数⼀般是⼆元函数f(x,y)或三元函数f(x,y,z)；形如：∫∫f (x,y )dxdy ba dc∫∫∫f(x,y,z)dxdydz b a d cf eMatlab 中有dblquad 函数和triplequad 函数来对上述两个积分实现.调⽤格式为： I=dblquad(@fun,a,b,c,d,tol)I=triplequad(@fun,a,b,c,d,e,f,tol)Fun 为被积函数，[a,b]为x 的积分区间；[c,d]为y 的积分区间；[e,f]为z 的积分区间.Dblquad 函数和triplequad 函数不允许返回调⽤的次数，如果需要知道函数调⽤的次数，则在定义被积函数的m ⽂件中增加⼀个计数变量，统计出被积函数被调⽤的次数.---------------------------------------------------------------------例如:计算⼆重积分I =∫∫√dxdy π2π2π2π2的值.先编写函数⽂件fxy.mfunction f=fxy(x,y) global k; k=k+1;f=sqrt(x.^2+y.^2);再调⽤函数dblquadglobalk; k=0;I=dblquad(@fxy,-pi/2,pi/2,-pi/2,pi/2,1.0e-10) I= 11.8629 k k= 37656---------------------------------------------------------------------例如：求三重积分∫∫∫4xze ?z2y?x 2dxdydz ππ1的值.编写函数⽂件fxyz1.mfunction f=fxyz1(x,y,z)global j;j=j+1;f=4*x.*z.*exp(-z.*z.*y-x.*x);调⽤triplequad函数editglobalj;j=0;I=triplequad(@fxyz1,0,pi,0,pi,0,1,1.0e-10)I=1.7328jj=1340978---------------------------------------------------------------------Matlab数值微分1.数值微分与差商导数的三种极限定义f′(x)=limn→0f(x+h)?f(x)hf′(x)=limn→0f(x)?f(x?h)f′(x)=limn→0f(x+h2)?f(x?h2)h上述公式中假设h>0，引进记号：f(x)=f(x+h)f(x)f(x)= f(x)f(xh)δf(x)= f(x+h)?f(x?h)称上述?f(x)、?f(x)、δf(x)为函数在x点处以h(h>0)为步长的向前差分、向后差分、中⼼差分，当步长h⾜够⼩时，有：f′(x)≈?f(x) hf′(x)≈f(x) f′(x)≈δf(x)f(x) h 、?f(x)h、δf(x)h也分别被称为函数在x点处以h(h>0)为步长的向前差商、向后差商、中⼼差商.当h⾜够⼩时，函数f(x)在x点处的导数接近于在该点的任意⼀种差商，微分接近于在该点的任意⼀种差分.2.函数导数的求法2.1⽤多项式或样条函数g(x)对函数f(x)进⾏逼近(插值或拟合)，然后⽤逼近函数g(x)在点x处的导数作为f(x)在该点处的导数.2.2⽤f(x)在点x处的差商作为其导数.3.数值微分的实现⽅法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按⾏计算差分.--------------------------------------------------------------------- 例如：⽣成6阶范德蒙德矩阵，然后分别按⾏、按列计算⼆阶向前差分A=vander(1:6)A=111111321684212438127931102425664164131256251252551777612962163661D2A1=diff(A,2,1)D2A1=180501220057011018200132019424200255030230200D2A2=diff(A,2,2)D2A2=000084211083612457614436920004008016540090015025--------------------------------------------------------------------- 例如：设f(x)=√x3+2x2?x+12+√(x+5)6+5x+2求函数f(x)的数值导数，并在同⼀坐标系中作出f’(x)的图像.已知函数f(x)的导函数如下：f′(x)=3x2+4x?12√x3+2x2?x+12+16√()56+5编辑函数⽂件fun7.m和fun8.m functionf=fun7(x)f=sqrt(x.^3+2*x.^2-x+12)+(x+5).^(1/6)+5*x+2;functionf=fun8(x)f=(3*x.^2+4*x-1)/2./sqrt(x.^3+2*x.^2-x+12)+1/6./(x+5).^(5/6)+5 ;x=-3:0.01:3;p=polyfit(x,fun7(x),5);⽤5次多项式拟合曲线dp=polyder(p);对拟合多项式进⾏求导dpx=polyval(dp,x);对dp在假设点的求函数值dx=diff(fun7([x,3.01]))/0.01;直接对dx求数值导数gx=fun8(x);求函数f的函数在假设点的导数plot(x,dpx,x,dx,'.',x,gx,'-')可以发现，最后得到的三条曲线基本重合.--------------------------------------------------------------------- 练习：A.⽤⾼斯-克朗罗德法求积分∫dx1+x2 +∞∞的值并讨论计算⽅法的精确度.（该积分值为π）function f=fun9(x)f=1./(1+x.^2);formatlong[I,err]=quadgk(@fun9,-Inf,Inf)I=err=B.设函数f(x)=sin x⽤不同的办法求该函数的数值导数，并在同⼀坐标系中作出f′(x)的图像.已知f′(x)=x cos x+cos x cos2x?sin x+2sin x sin2x()2function f=fun10(x)f=sin(x)./(x+cos(2*x));function f=fun11(x)f=(x.*cos(x)+cos(x).*cos(2*x)-sin(x)-2*sin(x).*sin(2*x))/(x+cos(2 *x)).^2; x=-3:0.01:3;p=polyfit(x,fun10(x),5);dp=polyder(p);dpx=polyval(dp,x);dx=diff(fun10([x,3.01]))/0.01;gx=fun11(x);plot(x,dpx,'r:',x,dx,'.g',x,gx,'-k')。
matlab中的微分方程的数值积分
MATLAB是一种流行的数学软件，用于解决各种数学问题，包括微分方程的数值积分。
微分方程是许多科学和工程问题的数学描述方式，通过数值积分可以得到微分方程的数值解。
本文将介绍在MATLAB中如何进行微分方程的数值积分，以及一些相关的技巧和注意事项。
一、MATLAB中微分方程的数值积分的基本方法1. 常微分方程的数值积分在MATLAB中，常微分方程的数值积分可以使用ode45函数来实现。
ode45是一种常用的数值积分函数，它使用4阶和5阶Runge-Kutta 方法来求解常微分方程。
用户只需要将微分方程表示为函数的形式，并且提供初值条件，ode45就可以自动进行数值积分，并得到微分方程的数值解。
2. 偏微分方程的数值积分对于偏微分方程的数值积分，在MATLAB中可以使用pdepe函数来实现。
pdepe可以求解具有定解条件的一维和二维偏微分方程，用户只需要提供偏微分方程的形式和边界条件，pdepe就可以进行数值积分，并得到偏微分方程的数值解。
二、在MATLAB中进行微分方程数值积分的注意事项1. 数值积分的精度和稳定性在进行微分方程的数值积分时，需要注意数值积分的精度和稳定性。
如果数值积分的精度不够，可能会导致数值解的误差过大；如果数值积分的稳定性差，可能会导致数值解发散。
在选择数值积分方法时，需要根据具体的微分方程来选择合适的数值积分方法，以保证数值解的精度和稳定性。
2. 初值条件的选择初值条件对微分方程的数值解有很大的影响，因此在进行微分方程的数值积分时，需要选择合适的初值条件。
通常可以通过对微分方程进行分析，或者通过试验求解来确定合适的初值条件。
3. 数值积分的时间步长在进行微分方程的数值积分时，需要选择合适的时间步长，以保证数值积分的稳定性和效率。
选择时间步长时，可以通过试验求解来确定合适的时间步长，以得到最优的数值解。
三、MATLAB中微分方程数值积分的实例以下通过一个简单的例子来演示在MATLAB中如何进行微分方程的数值积分。
MATLAB中的微积分运算（数值符号）
MATLAB中的微积分运算（数值符号）显然这个函数是单词differential（微分）的简写，⽤于计算微分。
实际上准确来说计算的是差商。
如果输⼊⼀个长度为n的⼀维向量，则该函数将会返回长度为n-1的向量，向量的值是原向量相邻元素的差，于是可以计算⼀阶导数的有限差分近似。
(1)符号微分1.常⽤的微分函数函数:diff(f) 求表达式f对默认⾃变量的⼀次微分值diff(f,x) 求表达式f对⾃变量x的⼀次积分值diff(f,n) 求表达式f对默认⾃变量的n次微分值diff(f,t,n)求表达式f对⾃变量t的n次微分值>> x=1:10x =1 2 3 4 5 6 7 8 9 10>> diff(x)ans =1 1 1 1 1 1 1 1 1例1:求矩阵中各元素的导数求矩阵[1/(1+a) (b+x)/cos(x)1/(x*y) exp(x^2)]对x的微分,可以输⼊以下命令A = sym('[1/(1+a),(b+x)/cos(x);1,exp(x^2)]');B = diff(A,'x')可得到如下结果:例2:求偏导数求的偏导数。
syms x y;f = x*exp(y)/y^2;fdx = diff(f,x)fdy = diff(f,y)可得到如下结果：例3:求复合函数的导数求的导数sym('x');y = 'x*f(x^2)'y1 = diff(y,'x')得到结果如下：例4:求参数⽅程的导数对参数⽅程求导syms a b tf1 = a*cos(t);f2 = b*sin(t);A = diff(f2)/diff(f1) %此处代⼊了参数⽅程的求导公式B = diff(f1)*diff(f2,2)-diff(f1,2)*diff(f2)/diff(f1)^3 %求⼆阶导数可得到如下结果：例5:求隐函数的导数求的⼀阶导数syms x yp = 'x*y(x)-exp(x+y(x))'%隐函数可进⾏整体表⽰%注意y（x）这种写法，它代表了y是关于x的函数p1 = diff(p,x)可得到如下结果：2.符号积分1符号函数的不定积分函数:int功能:求取函数的不定积分语法:int(f)int(f,x)说明:第⼀个是求函数f对默认⾃变量的积分值;第⼆个是求⾃变量f对对⾃变量t的不定积分值。
第7章 MATLAB数值微分与积分_习题答案
1第7章 MATLAB 数值微分与积分习题7一、选择题1．diff([10,15])的值是（ ）。
AA ．5B ．10C ．15D ．252．数值积分方法是基于（ ）的事实。
DA ．求原函数很困难B ．原函数无法用初等函数表示C ．无法知道被积函数的精确表达式D ．A ，B ，C 三个选项3．求数值积分时，被积函数的定义可以采取（ ）。
DA ．函数文件B ．内联函数C ．匿名函数D ．A ，B ，C 三个选项4．以下选项不能用来求数值积分的函数是（ ）。
BA ．quadgkB ．quad2C ．integralD ．integral25．以下选项不是离散傅里叶变换的函数是（ ）。
CA ．fftB ．fft2C ．fft1D ．fftn二、填空题1．在MATLAB 中，没有直接提供求 的函数，只有计算 的函数diff 。
数值导数，向前差分2．基于变步长辛普森法，MATLAB 给出了 函数和 函数来求定积分。
quad ，quadl3．MA TLAB 提供了基于全局自适应积分算法的 函数来求定积分，该函数的积分限 （可以或不可以）为无穷大。
integral ，可以4．MATLAB 提供的 、 、 函数用于求二重积分的数值解， 、 函数用于求三重积分的数值解。
integral2，quad2d ，dblquad ，integral3，triplequad5．MA TLAB 提供了离散傅里叶变换函数fft ，对应的逆变换函数是 。
ifft三、应用题1．求函数在指定点的数值导数。
（1）2346x x x x f 22ππππ,,,,cos sin)(=+= （2）321x 1x x f 2,,,)(=+=2(1):(2):直接用导数函数求:f=inline('x./sqrt(x.^2+1)');f(1)用拟合函数求:f=inline('sqrt(x.^2+1)');x=0:0.001:5;p=polyfit(x,f(x),5);dp=polyder(p);dpx=polyval(dp,1)2．求定积分。
matlab数值微分
h h 余项主部 f ( x + ) − f ( x − ) h2 2 2 2 − ⇒ f ′( x) = − f ′′′( x) + o(h2 ) — 一阶中心公式 24 h h 余项分别为： 余项分别为： − h f ′′( x) + o(h) 或 − h f ′′(ξ )； ξ − f ′′(ξ0 ) , ξ0 ∈[ x0 , x1 ] 2 2 2 2 h h2 h2 (3) 2 ξ f ′′′( x ) + o( h ) 或 − f ′′′(ξ )； − − f (ξ 1 ) , ξ 1 ∈[ x0 , x1 ] 24 24 24 2 2 h h h f ′′( x) + o(h) 或 f ′′(ξξ)。 + f ′′(ξ1 ) , ξ1 ∈ [ x0 , x1 ] 2 2 2
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
科学计算—理论、方法及其基于MATLAB 的实现与分析数值微积分§1 数值微分对于给定的函数()x f y =，如果 １、()x f y =的函数关系式比较复杂时； ２、()x f y =未知，而仅仅知道()x f y =在1+n 个相异点kx ，n k ,,1,0Λ=处的函数值k y ；则希望能用相对简单的计算方法，求得()x f y =导数的（近似）值。
基于上述考虑，选择的方法之一是利用函数()x f y =的插值多项式的导数作为函数()x f y =导数的近似值，例如Lagrange 插值多项式，由于()()()()∑==≈nk k k n x f x l x L x f 0（１）因而有()()x L x f n'≈' （２）这里需要说明一点的是，尽管()x f 和()x L n 的函数值可能相差不多，但是它们的导数有可能相差很大，如下面的例子 例1： 考虑函数()22511x x f +=在区间[-1，1]的插值问题，取区间[-1，1]的11个点kx k ⨯+-=2.01，10,,1,0Λ=k ，作函数()22511x x f +=的10次插值多项式：()18552.163597.1234338.3819095.4949417.220246810+-+-+-=x x x x x x L n函数()x f 和插值多项式()x L n 的导数分别为 ()()2225150x xdxx df +-=()x x x x x dxx dLn 7.334.4936.22883.39594.22093579-+-+-= 对函数()x f 和插值多项式()x L n 及其导数分别比较，结果如图所示：Derivative_Runge下面我们基于Taylor 公式()()()()()()()()()()⎥⎦⎤⎢⎣⎡≅++++''+'+=++++1112!1!!2n n n n n h O h n f h n x f h x f h x f x f h x f ξΛ (3)讨论导数的近似计算问题 1 一阶导数的近似计算令k k x x h -=+1,可得一阶向前有限差商公式(First Forward Finite Divided Difference):()()()()()()()()()()()()()()⎪⎩⎪⎨⎧=-≈'⇒+''--='⇒+''+'+=+++h O error h x f x f x f h O h x f h x f x f x f h O h x f h x f x f x f k k k k k k k k k k k 121321!2!2 (4)类似地,当1--=k k x x h 时,可得一阶向后有限差商公式(First Backward Finite Divided Difference):()()()()()()()()()()()()()()⎪⎩⎪⎨⎧=-≈'⇒+''+-='⇒+''+'-=---h O error h x f x f x f h O h x f h x f x f x f h O h x f h x f x f x f k k k k k k k k k k k 121321!2!2 (5)由近似计算公式(4)和(5),可得可得一阶中心有限差商公式(First Centered Finite Divided Difference):()()()()()()()()()()⎪⎩⎪⎨⎧=-≈'⇒+-='⇒+-+-+2112112254h O error hx f x f x f h O hx f x f x f k k k k k k (6)2 二阶导数的近似计算在Taylor 公式中，用h 2代替h 可得()()()()()()3222!22h O h x f h x f x f x f k k k k +''+'+=+ (7)式(7)与式(4)结合可得二阶向前有限差商公式(Second Forward Finite Divided Difference):()()()()()32122h O h x f x f x f x f k k k k ++''+-=-++ (8)()()()()()⎪⎩⎪⎨⎧=+-≈''⇒++h O error h x f x f x f x f k k k k 2122 ()()()()()()⎪⎩⎪⎨⎧=---≈''⇔+++h O error h h x f x f h x f x f x f k k k k k 112 (9) 类似地利用()()()()()()3222!22h O h x f h x f x f x f k k k k +''+'-=-()()()()()321!2h O h x f h x f x f x f k k k k +''+'-=-可得二阶向后有限差商公式(Second Backward Finite Divided Difference):()()()()()32122h O h x f x f x f x f k k k k ++''+-=---()()()()()⎪⎩⎪⎨⎧=+-≈''⇒-+h O error h x f x f x f x f k k k k 2212 (10) 进一步底，由公式()()()()()()()()()()()()43214321!3!2!3!2h O h x f h x f h x f x f x f h O h x f h x f h x f x f x f k k k k k k k k k k +'''-''+'-=+'''+''+'+=-+ (11)可得二阶中心有限差商公式(Second Centered Finite Divided Difference):()()()()()⎪⎩⎪⎨⎧=+-≈''-+22112hO error h x f x f x f x f k k k k (12)§2 数值积分§2.1 Newton-Cotes 数值积分公式基于数值微分的同样原因，当函数()x f y =在],[b a 上的定积分不易计算时，利用函数()x f y =的Lagrange 插值多项式()()()∑==nk k k n x f x l x L 0，有()()()()k nk ba k ban b ax f dx x l dx x L dx x f I ∑⎰⎰⎰=⎥⎦⎤⎢⎣⎡=≈=0 (13) 由于Lagrange 插值多项式()()()∑==nk k k n x f x l x L 0的插值基函数()x l k 只依赖插值节点k x ，n k ,,1,0Λ=，所以当k x ，n k ,,1,0Λ=取定后，()x l k 就是完全确定的多项式函数，令()⎰=bak k dx x l A ， n k ,,1,0Λ= (14)则由式(13)得到Newton-Cotes 求积公式：()()()k nk k ban b ax f A dx x L dx x f I ∑⎰⎰==≈=0(15)特别地，当取插值节点为kh a x k +=，n k ,,1,0Λ=，nab h -=(16) 时有１） 两点公式（Trapezoidal Rule ）：⇒=1n 2210ab h A A -=== ()()()()()[]b f a f ab x f A dx x L dx x f I k nk k baba+-==≈=∑⎰⎰=21 (17)2） 三点公式（Simpson's 1/3 Rule ）：34,322120hA h A A a b h n ===⇒-=⇒=，()()()()()()⎥⎦⎤⎢⎣⎡+⎪⎭⎫ ⎝⎛++=⎥⎦⎤⎢⎣⎡+⎪⎭⎫ ⎝⎛++-=≈=⎰⎰b f b a f a f h b f b a f a f a b dxx L dx x f I baba2432462 (18)3） 四点公式（Simpson's 3/8 Rule ）：89,83332130hA A h A A a b h n ====⇒-=⇒=()()()()()()[]()()()()[]b f h a f h a f a f hb f h a f h a f a f ab dxx L dx x f I bab a+++++=+++++-=≈=⎰⎰2338323383 (19)4） 误差估计：利用Lagrange 插值的误差公式：()()()()()()()∏=+-+=-=nk k n n n x x n f x L x f x R 01!1ξ (20)容易得到Trapezoidal Rule 的误差估计：()()()()()()()()()()()32321121222a b M a b f dx ab x b a x f dx b x a x f dxx L dx x f b a ba bab a-≤-''=++-''≈--''=-⎰⎰⎰⎰ξξξ (21) Simpson's 1/3 Rule 以及Simpson's 3/8 Rule 的误差估计分别为：()()()()5445229012901⎪⎭⎫⎝⎛-≤⎪⎭⎫ ⎝⎛--≈-⎰⎰a b M f a b dx x L dx x f bab a ξ (22) ()()()()5445338033803⎪⎭⎫⎝⎛-≤⎪⎭⎫ ⎝⎛--≈-⎰⎰a b M f a b dx x L dx x f babaξ (23)其中()(){}x f M k b a k],[max =.注１ 数值积分公式(13)表明，数值积分方法最大的优点是将复杂的函数积分转化为相对简单的加、减、乘、除运算； 注２ 从几何的角度看，上述三种求积公式（17）、（18）和（19）（在积分区间],[b a 上）分别用直线()()b f ab ax a f b a b x y --+--=（24）抛物线()()()()()()()()()b f h h a x a x h a f h b x a x a f h h a x h a x y 2222---++-------=（25）三次立方曲线()()()()()()()()()()()()()()()()h a f h h a x h a x a x h a f h h a x h a x a x h a f h h a x h a x a x a f h h a x h a x h a x y 332332323333+-----++------+-----+-------= （26）代替曲线()x f y =:例1 计算积分⎰=10dx e I x ，并估计误差.解 1） 用梯形公式计算，8591409.1)(20110≈+-=≈e e T I .由于x e x f =)(，所以x e x f ='')(，于是，梯形公式的误差23.012)(max 12110≈=''≤-=≤≤ex f T I R x T .2） 用辛普森公式计算，7188612.1)4(60112/10≈++-=≈e e e S I .由于x ）（e x f=)(4，于是，辛普森公式的误差00094.02880)(max 28801410≈=≤-=≤≤ex f S I R ）（x S .【注】 （1） 当0)()2(=x f 时，梯形求积公式准确成立；当)()4(x f=0时，辛普森公式准确成立.即梯形公式对一次多项式准确成立，而辛普森公式对三次多项式准确成立.（2） 一般地，由余项公式(7.9)知，当)(x f 是n 次多项式时，积分余项为零，从而牛顿-柯特斯求积公式准确成立. （3） 数值求积方法是一种近似方法，因此，要求求积公式对尽可能多的被积函数)(x f 能准确计算积分值.作为衡量公式逼近好坏的标准之一，下面给出代数精度(the Degree of Precision of the Quadrature Formula)的概念.§2.1.2 代数精度【定义1】 如果某求积公式对于次数小于等于m 的多项式能准确求出积分值，而对某个1+m 次多项式就不能准确求出积分，则称该求积公式具有m 次代数精度.由定义1，具有m 次代数精度的求积公式（7.3）对m x x x f ,,,1)(Λ=时精确成立，而对1)(+=m x x f 不能精确成立.为了构造形如公式（7.3）的求积公式，通过解方程组⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧+-==-==-==∑⎰∑⎰∑⎰=++==n k b a m m mm k k n k b a k k nk b a k m a b dx x x A a b xdx x A a b dx A 0110220,1,2,1Λ可得求积节点k x 与求积系数k A ，由此求得具有m 次代数精度的求积公式.不难验证，梯形公式和矩形公式均具有一次代数精度，而辛普森公式具有三次代数精度.对牛顿-柯特斯求积公式，有下面结论.【定理1】 牛顿-柯特斯求积公式（7.4）至少具有n 次代数精度，而当n 为偶数时，至少具有1+n 次代数精度.证明 当)(x f 为任何次数不高于n 的多项式时，0)()1(≡+x f n ，所以0)(=⎰ban dx x R ,显然结论成立. 当n为偶数时，只须对1)(+=n x x f 时的结论验证.因为)!1()()1(+=+n x fn ，由截断误差公式⎰⎰⎰∏∏=+=-=-=baban nj n nj jn dtj t hdx xx dx x R 002)()()(，若令2n u t +=，则有⎰⎰∏-=+-+=ban n n j n n du j nu hdx x R 2202)2()(.注意到∏=-+nj j n u 0)2(是u的奇函数（n为偶数），因此⎰=ban dx x R 0)(，结论成立.例2 对于求积公式20011()(1)(2)'(2)f x dx A f B f B f -=-++⎰，试确定100,,B B A ，使此求积公式有最高的代数精度，是几阶的？§2.2 复化求积公式由前面分析知道，在整个积分区间],[b a 上用直线或用抛物线代替原曲线()x f y =做积分，尽管计算简单，但精度较差，结果不好，因此，一个自然的想法是：在尽可能保持公式(17)、和(18)计算简单、易于计算机实现的优点的前提下，寻找提高误差精度的途径和方法。