matlab数值微分与数值积分

从本节开始，我们将进入matlab的数值微积分与方程求解模块，一起学习如何利用matlab去解决微积分问题。
对于本节内容，主要分为两个部分讲解，数值微分和数值积分，那么下面，就开始今天的学习吧！
一、数值微分
在正式开始之前，有几个新概念需要讲解一下
（1）数值差分与差商：
微积分中，任意函数 f(x) 在x0 点的导数是通过极限定义的（如下图所示）：

如果去掉极限定义中h趋向于0的极限过程，得到函数在x0 点处以h（ h>0 ）为步长的向前差分、向后差分和中心差分公式（如下图）：

当步长h充分小时，得到函数在x0 点处以h（ h>0 ）为步长的向前差商、向后差商和中心差商公式（如下图）：

函数 f(x) 在点x0的微分接近于函数在该点的差分，而f在点x的导数接近于函数在该点的差商。
（2）数值微分的实现：
MATLAB 提供了求向前差分的函数diff ，其调用格式有三种：
dx=diff(x) ：计算向量x的向前差分， dx(i)=x(i+1)- x(i) ， i=1，2， … ，n, 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) 在某点处的差商作为其导数的近似值。
例、设f(x)=ln(x)，采用上述方法绘制出[0，10]之间的其导函数的图像，求解出导函数在1处的值，对比理论值1之间的差值，看相差大不大。

>>  x=0:0.01:10;
y=log(x);
f1=diff(y)./diff(x);
plot(x(1:end-1),f1);

打开matlab的工作区查看f1在100处的值，即就是导函数在1出=处的值（如图）：

值为1.005，与理论值1相差0.005，可能是由于x=0:0.01:10;这一语句导致，发现当取点越多时，结果越接近1，再看图像，其符合1/x的函数图像。
二、数值积分
1、数值积分基本原理
在高等数学中，计算定积分依靠微积分基本定理，只要找到被积函数 f(x)的原函数F(x)，则可用牛顿 — 莱布尼兹（ Newton- – Leibniz ）公式：

在有些情况下，应用牛顿—莱布尼兹公式有困难，例如，当被积函数的原函数无法用初等函数表示，或被积函数是用离散的表格形式给出的。这时就需要用数值解法来求定积分的近似值。
求定积分的数值方法多种多样，如梯形法、辛普森（ Simpson ）法、高斯求积公式等。它们的基本思想都是将积分区间 [a ， b] 分成n个子区间 [xi，xi+1] ，i=1 ，2 ，… ，n，其中x1 =a ，xn+1=b ，这样求定积分问题就分解为下面的求和问题。

在每一个小的子区间上定积分的值可以近似求得, , 从而避免了牛顿 — 莱布尼兹公式需要寻求原函数的困难。
2、数值积分的实现
（1）基于自适应辛普森方法：
[ I,n ]=quad( filename,a,b,tol,trace)
（2）基于自适应 Gauss- – Lobatto方法：
[ I,n ]= quadl ( filename,a,b,tol,trace)
其中， filename 是被积函数名；a和b分别是定积分的下限和上限，积分限 [a ，b]必须是有限的，不能为无穷大（ Inf ）； tol 用来控制积分精度，默认时取tol =10^(- 6 )； trace 控制是否展现积分过程，若取非0则展现积分过程，取0则不展现，默认时取 trace=0 ；返回参数I即定积分的值，n为被积函数的调用次数。
例、分别用 quad 函数和 quadl 函数求定积分的近似值，并在相同的积分精度下，比较被积函数的调用次数。求解在[1，10]上，sin(x)的定积分值。

>>  format long
f=@(x)cos(x);
[ I,n ]=quad(f,1,10,1e-8)
I =-1.385492095697597
n =229
>> [ I,n ]=quadl(f,1,10,1e-8)
I =-1.385492095687089
n =138
>> sin(10)-sin(1)
ans =-1.385492095697266

可以看出，两种方式求出的定积分值与理论值sin(10)-sin(1)几乎一致，但两种函数求解时的被积函数的调用次数却有差别。
（3）基于全局自适应积分方法：
I=integral( filename,a,b )
其中，I是计算得到的积分；filename 是被积函数；a和b分别是定积分的下限和上限，积分限可以为无穷大。
（4）基于自适应高斯- – 克朗罗德方法
[ I,err ]= quadgk( filename,a,b)
其中，err返回近似误差范围，其他参数的含义和用法与quad函数相同。积分上下限可以是无穷大（−Inf 或 Inf），也可以是复数。如果积分上下限是复数，则 quadgk函数在复平面上求积分。
上述四种积分的求解方法各有不同，同时也各有优劣，一、二两种方法要求积分上下限是有限的，不可取无穷值，三、四两种方法中积分上下限可以取无穷值，且第四种方法还可以在复平面上求积分。
下面介绍一种特殊的积分方法：梯形积分法
I= trapz( x,y)
其中，向量x、y 定义函数关系y=f(x)，值得注意的是，此处的x，y表示为向量形式。
例、 5 设 x=1:10，y=[5,4,1,9,3,3,5,8,7,2], 用 trapz函数计算定积分并且画图。

>> x=1:10;y=[5,4,1,9,3,3,5,8,7,2];
>> plot(x,y);
>> I1=trapz(x,y)
I1 =43.500000000000000

三、多重定积分的数值求解
现将多重积分的求解函数总结如下，它们的用法和上述的一重积分用法基本一致
求二重积分的数值解：

I=integral2( filename,a,b,c,d)
I=quad2d( filename,a,b,c,d)
I= dblquad( filename,a,b,c,d,tol)
求三重积分的数值解：

I=integral3( filename,a,b,c,d,e,f) )
I= triplequad( filename,a,b,c,d,e,f,tol)
本节内容就讲述到这里，下节将讲述线性方程组求解，敬请期待！