[MATLAB]数值微分与数值积分(diff/quad)

数值微分
(1)数值差分与差商

如果去掉x->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,按行计算差分。

例子：设f(x)=sinx,在[0,2Π]范围内随机采样，计算f’(x)的近似值，并与理论值f’(x)=cosx进行比较.

>> x=[0,sort(2*pi*rand(1,5000)),2*pi];
>> y=sin(x);
>> f1=diff(y)./diff(x);
>> f2=cos(x(1:end-1));
>> plot(x(1:end-1),f1,x(1:end-1),f2);
>> d=norm(f1-f2)

二、数值积分
(1)数值积分基本原理
    牛顿–莱布尼兹公式：


(2)数值积分的实现

• 基于自适应辛普森方法

 [l,n]=quad(filename,a,b,tol,trace)

• 基于自适应Gauss-Lobatto方法

 [l,n]=quadl(filename,a,b,tol,trace)

其中，filename是被积函数；a和b分别是定积分下限和上限，积分限[a,b]必须是有限的，不能为无穷大(inf);tol用来控制积分精度，默认时取tol=10^(-6);trace控制是否展现积分过程，若取非0则展现积分过程，取0则不展现，默认时取trace=0；返回参数i即定积分的值，n为被积函数的调用次数

>> format long
>> f=@(x) 4./(1+x.^2);
>> [I,n]=quad(f,0,1,1e-8);
>> [I,n]=quad(f,0,1,1e-8)

I =

   3.14159265373344


n =

    61

>> [I,n]=quadl(f,0,1,1e-8)

I =

   3.14159265358981


n =

    48

>> (atan(1)-atan(0))*4

ans =

   3.14159265358979
>> 

(2)数值积分的实现

• 基于全局自适应积分方法

I=integral(filename,a,b)

其中，I是计算得到的积分；filename是被积函数；a和b分别是定积分的上限和下限，积分限可以为无穷大。

• 基于自适应高斯-克朗罗德方法

[I,err]=quadgk(filename,a,b)

其中，err返回近似误差范围，其他参数的含义和用法与quad函数相同。积分上下限可以是无穷大（-lnf或lnf），也可以是复数。如果积分上下限是复数，则quadgk函数在复平面上求积分。

>> f=@(x) sin(1./x)./x.^2;
>> I=quadgk(f,2/pi,+Inf)

• 基于梯形积分法
  

I=trapz(x,y)

其中，向量x、y定义函数关系y=f(x).

可用以下语句实现

sum(diff(x).*(y(1:end-1)+y(2:end))/2)

>> x=1:6;
>> y=[6,8,11,7,5,2];
>> plot(x,y,'-ko');
>> grid on
>> axis([1,6,0,11]);
>> ll=trapz(x,y)

ll =

    35

>> I2=sum(diff(x).*(y(1:end-1)+y(2:end))/2)

I2 =

    35

（3）多重定积分的数值求解

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)

>> f1=@(x,y) exp(-x.^2/2).*sin(x.^2+y);
>> l1=quad2d(f1,-2,2,-1,1);
>> f2=@(x,y,z) 4*x.*z.*exp(-z.*z.*y.*x);
>> l1=integral3(f2,0,pi,0,pi,0,1);