MATLAB-6-1数值微分与数值积分

1. 数值微分

1） 求向前差分的函数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的ｎ阶差分，dim=1时（默认状态），按列计算差分；dim=2时，按行计算差分。
  e.g.设f(x)=sinx，在[0,2pi]范围内随机采样，计算f’(x)的近似值，并与理论值f’(x)=cos(x)进行比较。

>> x=[0,sort(2*pi*rand(1,5000)),2*pi];
>> y=sin(x);
>> f1=diff(y)./diff(x);
%f'(x)的表示方法
>> f2=cos(x(1:end-1));
>> plot(x(1:end-1),f1,x(1:end-1),f2);
>> d=norm(f1-f2)
%范数norm(A)返回A的最大奇异值，即max(svd(A))
d =
    0.0456

2. 数值积分

1）数值积分基本原理

2）数值积分的实现

i. 基于全局自适应积分方法

I=integral(filename,a,b)
I：定积分值
filename：被积函数
a，b：定积分下限、上限，积分限可以无穷大
e.g.求定积分

%先定义被积函数文件fe.m
function f = fe(x)
f = 1./(x.*sqrt(1-log(x).^2));
%注意使用点除，点乘方
end
%命令行
>> I=integral(@fe,1,exp(1))   %文件名
I =
1.5708

ii. 基于自适应辛普森方法

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

iii. 基于自适应Gauss-Lobatto方法

[I,n]=quadl(filename,a,b,tol,trace)
filename：被积函数名；
a，b为定积分的下限和上限，积分限[a,b]必须是有限的，不能为无穷大（inf）
tol：控制积分精度，默认取 tol=10^(-6)
trace：控制是否展现积分过程，若取非0则展现积分过程，取0则不展现，默认取0
I：返回参数，即定积分值
n：被积函数的调用次数
e.g.分别用quad函数和quadl函数求定积分的近似值，并在相同的积分精度下，比较被积函数的调用次数。

>> 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
>> format short

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

[I,err]=quadgk(filename,a,b)
err：近似误差范围
积分上下限可以是无穷大（-Inf或Inf），也可以是复数
若积分上下限为复数，则quadgk( )函数在复平面上求积分
e.g.求定积分
被积函数文件fe.m

function f = fe(x)
f = sin(1./x)./(x.^2);
end

命令行

>> I=quadgk(@fe,2/pi,+Inf)    %文件名
I =
	1.0000

v. 基于梯形积分法

已知$(x_i,y_i)$（i=1,2,…,n），且$a=x_1<x_2<\dots <x_n=b$，求$I={\int }_a^{b}f(x)d(x)$近似值。I=trapz(x,y)，其中，向量x，y定义函数关系y=f(x)。
原理：trape( )函数采用梯形积分法则，积分近似值为：

可用以下语句实现：sum(diff(x).*(y(1:end-1)+y(2:end))/2)
e.g.设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]);
%设值x轴数值范围为1-6，y轴数值范围为0-11
>> I1=trapz(x,y)
I1 =
    35
>> I2=sum(diff(x).*(y(1:end-1)+y(2:end))./2)
I2 =
    35

3. 多重定积分的数值求解

1）求二重积分的数值解${\int }_c^d{\int }_a^{b}f(x,y)dxdy$

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

2）求三重积分的数值解${\int }_e^f{\int }_c^d{\int }_a^{b}f(x,y,z)dxdydz$

I=integral3(filename,a,b,c,d,e,f)
I=triplequad(filename,a,b,c,d,e,f,tol)
e.g.分别求二重积分和三重积分

>> f1=@(x,y) exp(-x.^2/2).*sin(x.^2+y);
>> I1=quad2d(f1,-2,2,-1,1)
I1 =
    1.5745
>> f2=@(x,y,z) 4*x.*z.*exp(-z.^2.*y-x.^2);
>> I2=integral3(f2,0,pi,0,pi,0,1)
I2 =
    1.7328
%注意：1）函数定义方式；2）函数参数个数；3）积分下限、上限书写顺序