数值微分与数值积分

数值微分

注意：diff函数计算的是向量元素间的差分，故差分向量元素的个数比原向量少了
一个。同样，对于矩阵来说，差分后的矩阵比原矩阵少了一行或一列。 
另外，计算差分之后，可以用f(x)在某点处的差商作为其导数的近似值。 
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，按行计算差分。

数值积分

 基于自适应辛普森方法 
[I,n]=quad(filename,a,b,tol,trace) 
 基于自适应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为被积函数的调用次数。

>> format long 
>> f=@(x) 4./(1+x.^2); 
>> [I,n]=quad(f,0,1,1e-8) 
I = 
 3.141592653733437 
n = 
 61 
>> [I,n]=quadl(f,0,1,1e-8) 
I = 
 3.141592653589806 
n = 
 48 
>> (atan(1)-atan(0))*4 
ans = 
 3.141592653589793 
>> format short

 基于全局自适应积分方法 
I=integral(filename,a,b) 
其中，I是计算得到的积分；filename是被积函数；a和b分别是定积分的下限和
上限，积分限可以为无穷大。

>> f=@(x) 1./(x.*sqrt(1-log(x).^2));
>> I=integral(f,1,exp(1))

I =

    1.5708

 基于自适应高斯-克朗罗德方法 
[I,err]=quadgk(filename,a,b) 
其中，err返回近似误差范围，其他参数的含义和用法与quad函数相同。积分
上下限可以是无穷大（−Inf或Inf），也可以是复数。如果积分上下限是复数，
则quadgk函数在复平面上求积分。

例5 设x=1:6，y=[6,8,11,7,5,2],用trapz函数计算定积分。

x=1:6;
y=[6,8,11,7,5,2];
plot(x,y,’-ko’);
grid on
axis([1,6,0,11]);
I1=trapz(x,y)
I1 =
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); 
>> I1=quad2d(f1,-2,2,-1,1) 
I1 = 
  1.574 
>> f2=@(x,y,z) 4*x.*z.*exp(-z.*z.*y-x.*x); 
>> I2=integral3(f2,0,pi,0,pi,0,1) 
I2 = 
  1.7328