[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);

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值