6.1 matlab数值微分与数值积分

数值微积分适合求解没有或很准求出微分或积分表达式的问题的计算。
1、数值微分
(1）数值差分与差商
任意函数f(x)在x0点的导数是通过极限定义的:

如果去掉极限定义中h趋向于0的极限过程，得到函数在x0点处以h为步长的向前差分。
当步长h充分小时得到函数在x0点处，以h为步长的向前差商
函数f(x)在x0点的微分接近于函数在该点的差分。而f(x)在x0点的导数接近于函数在该点的差商。

(2）数值微分的实现
MATLAB提供了求 向前差分的函数diff，其调用格式有三种:
（1）dx=diff(x): 计算向量x的一阶向前差分,dx(i)=x(i+1)-x(i)，i=1，
2，…，n-1。n是向量x的元素的个数

（2）dx=diff(x,n): 计算向量x的n阶向前差分。例如，diff(x,2)=diff(diff(x)）：x的二阶差分等于x的一阶差分再求一阶差分

（3）dx=diff(A,n,dim): 计算矩阵A的n阶差分，dim=1时（默认状态)按列计算差分;dim=2，按行计算差分。

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

例1设f(x)=sin x，在[0，2π]范围内随机采样,计算f’(x)的近似值，并与
理论值f’(x)=cos x进行比较。

x = [0,sort(2*pi*rand(1,5000)),2*pi];   %确定x向量，其首尾元素分别是0和2π，中闻是0到2π开区间的随机数并按照从小到大排序
y = sin(x); 
f1 = diff(y)./diff(x);   %利用diff函数求一阶向前差分，y的差分除以x的差分，得到差商向量f1。差商作为导数的近似值
f2 = cos(x(1:end-1));   %求各点导数的理论值向量f2
plot(x(1:end-1),f1,'--',x(1:end-1),f2,'r:');   %绘制近似值和理论值曲线
legend('近似值曲线','理论值曲线','location','southeast')

d = norm(f1-f2)   %求近似值和理论值之差的范数，其值较小，说明近似值和理论值比较接近

%输出结果
d =

    0.0439

2、数值积分

(1）数值积分基本原理
牛顿—莱布尼兹(Newton-Leibniz）公式:
当被积函数的原函数无法用初等函数表示或被积函数是用表格形式给出的，这时，就需要用数值解法来求定积分的近似值。求定积分的数值方法多种多样：如梯形法、辛普森法、高斯求积公式等。它们的基本思想都是将积分区间[a b]分成n个子区间[xi xi+1]，这里i从1变化到n，其中x1等于a ，xn+1等于b，这样求定积分问题就分解为下面的求和问题:
(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;返回参数l即定积分的值,n为被积函数的调用次数。

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

④基于自适应高斯-克朗罗德方法 （求震荡函数的积分）
[l,err]=quadgk(filename,a,b)
其中，err返回近似误差范围，其他参数的含义和用法与quad函数相同。积分上下限可以是无穷大(-lnf或Inf)，也可以是复数。如果积分上下限是复数，则quadgk函数在复平面上求积分。
⑤基于梯形积分法
l=trapz(x,y)
其中，向量x、y定义函数关系y=f(x)。x，y是两个等长的向量，从x1变化到xn，y对应从y1变化到yn，并且x1小于x2一直到小于xn，积分区间是从x1到xn

例2：分别用quad函数和quadl函数求定积分的近似值，并在相同的积分精度下，比较被积函数的调用次数。

>>format long;
f = @(x)4./(1+x.^2);   %定义被积函数
[I1,n1] = quad(f,0,1,1e-8)  %调用quad函数求定积分,精度取10的-8次方
[I2,n2] = quadl(f,0,1,1e-8)  %调用quadl函数求定积分
y = (atan(1)-atan(0))*4   %输出积分的理论值
format short;



%输出结果
I1 =
   3.141592653733437

n1 =
    61

I2 =
   3.141592653589806

n2 =
    48

y =
   3.141592653589793

例3：求定积分。

>>f = @(x)1./(x.*sqrt(1-log(x).^2));   %定义被积函数
I = integral(f,1,exp(1))     %调用integral函数求定积分

%输出结果
I =
    1.5708

例4：求定积分。

>>f = @(x)sin(1./x)./(x.^2);   %定义被积函数
I = quadgk(f,2/pi,+Inf)   %调用函数quadgk求定积分


%输出结果
I =
    1.0000

(3)多重定积分的数值求解
①求二重积分的数值解:

l=integral2(filename,a,b,c,d)
l=quad2d(filename,a,b,c,d)
l=dblquad(filename,a,b,c,d,tol)

②求三重积分的数值解:

l=integral3(filename,a,b,c,d,e,f)
l=triplequad(filename,a,b,c,d,e,f,tol)

例6：分别求二重积分和三重积分。

>> f1 = @(x,y)sin(x.^2+y).*exp(-x.^2/2);   %定义被积函数f1
f2 = @(x,y,z)4*x.*z.*exp(-z.^2.*y-x.^2);  %定义被积函数f2
I1 = quad2d(f1,-2,2,-1,1)   %调用函数quad2d求二重定积分
I2 = integral3(f2,0,pi,0,pi,0,1)   %调用函数integral3求三重定积分


%输出结果
I1 =
    1.5745

I2 =
    1.7328