点击蓝字

关注我们

01

数值微分与数值积分

一、数值微分

(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-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)=sin x，在[0,2π]范围内随机采样，计算f’(x)的近似值，并与理论值f’(x)=cos x进行比较。

>> 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)d =    0.0456

二、数值积分

(1)数值积分基本原理

在高等数学中，计算定积分依靠微积分基本原理，只要找到被积函数f(x)的原函数F(x)，则可用牛顿—莱布尼兹(Newton-Leibinz)公式：

在有些情况下，应用牛顿—莱布尼兹公式有困难，例如，当被积函数的原函数无法用初等函数表示，或被积函数是用离散的表格形式给出的。这是就需要用数值解法来求定积分的近似值。

求定积分的数值方法多种多样，如梯形法、辛普森法，高斯求积公式等。它们的基本思想都是将积分区间[a,b]分成n个子区间[xi,xi+1],i=1,2,...,n，其中x1=a,xn+1=b，这样求定积分问题就分解为下面的求和问题：

在每一个小的子区间上定积分的值可以近似求得，从而避免了牛顿—莱布尼兹公式需要寻求原函数的困难。

(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为被积函数的调用次数。

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

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

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

I=intergral(filename,a,b)

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

【例】求定积分。

被积函数文件fe.m：

function f=fe(x)f=1./(x.*sqrt(1-log(x).^2));

>> I=intergral(@fe,1,exp(1))>> I=integral(@fe,1,exp(1))I =    1.5708

④基于自适应高斯—克朗罗德方法

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

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

【例】求定积分。

被积函数文件fsx.m：

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

>> I=quadgk(@fsx,2/pi,+Inf)I =    1.0000

⑤基于梯形积分方法

已知(xi,yi)(i=1,2,...,n)，且a=x1，求近似值。

I=trapz(x,y)

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

trapz函数采用梯形积分法则，积分的近似值为：

其中，hi=xi+1-xi

可以用以下语句实现：

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

【例】设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.5745>> 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

参考课程：中南大学《科学计算与MATLAB语言》

如有侵权，请联系删除

微信公众号｜技术人的文学小天地

微信号｜15035869237