MATLAB求分数阶微分的数值解，G-L定义，R-L定义，Caputo定义

分数阶微积分学是整数阶微积分学的直接拓展，将一阶导数、二阶导数、一重积分、二重积分等整数阶微积分拓展到0.75阶导数、 $\sqrt{2}$ 阶导数等实数甚至是复数阶的导数或积分。这无疑拓展了微积分学的深度。

对于整数阶微积分，一般可以具有简洁明确的物理意义，比如位移、速度和加速度可以很好地解释一个信号与其整数阶导数之间的关系。然而分数阶微积分却没有那么简洁易懂的物理解释。目前对于分数阶微积分的定义，比较应用广泛的是G-L定义，R-L定义和Caputo定义。

Grünwald-Letnikov定义：

用MATLAB语言编写出Grünwald-Letnikov分数阶微积分数值计算的函数：

function dy = glfdiff(y,t,gam)if strcmp(class(y),'function_handel')y = y(t);
endh = t(2)-t(1);
w = 1;
y = y(:);
t = t(:);for j = 2:length(t)w(j) = w(j-1)*(1-(gam+1)/(j-1));
end
for i = 1:length(t)dy(i) = w(1:i)*[y(i:-1:1)]/h^gam;
end

该函数的调用格式为：y1=glfdiff(y,t,r)。其中t为等间距的时间向量，r为阶次，y可以为函数f(t)在t各点的上的函数值，也可以为f(t)的函数句柄，若r为负值，则计算原函数的-r阶积分。

Riemann-Liouville定义：

f(t)的分数阶积分定义：

f(t)的分数阶微分定义应基于积分定义给出：

下面给出一种R-L定义的分数阶微积分的数值解法的函数：

function [dy,t] = rlfdiff(y,t,r)h = t(2)-t(1);
n = length(t);
dy1 = zeros(1,n);
y = y(:);
t = t(:);if r>-1,m = ceil(r)+1;p = m-r;y3 = t.^(p-1);
elseif r==-1n = length(t);yy = 0.5*(y(1:n-1)+y(2:n)).*diff(t);dy(1) = 0;for i=2:ndy(i) = dy(i-1)+yy(i-1);endreturn
elsem=-r;y3 = t.^(m-1);
endfor i = 1:ndy1(i) = y(i:-1:1).'*(y3(1:i));
endif r>-1dy = diff(dy1,m)/(h^(m-1))/gamma(p);t = t(1:n-m);
elsedy = dy1*h/gamma(m);
end

该函数的调用格式为[y1,t]=rlfdiff(y,t,r)。其调用格式类似刚刚的glfdiff()函数，不同的是该函数返回了两个向量y1和t。因为该函数使用了MATLAB求差分的函数diff()，向量y1的实际长度可能小于y的长度。

Caputo定义：

事实上Caputo分数阶微积分的定义与R-L定义是相同的，对于r阶微积分，记m=[r]，则二值之间的关系为：

Caputo分数阶微积分的数值解函数：

function dy = caputo(y,t,alfa,y0,L)dy = glfdiff(y,t,alfa);
if nargin<=4L = 10;
end
if alfa>0q = ceil(alfa);if alfa<=1y0 = y(1);endfor k = 0:q-1dy = dy-y0(k+1)*t.^(k-alfa)./gamma(k+1-alfa);endyy1 = interp1(t(L+1:end),dy(L+1:end),t(1:L),'spline');dy(1:L) = yy1;
end

该函数的调用格式为y1=caputo(y,t,r,y0,L)。当r<=0时，则可以返回Caputo分数阶微积分；如果r<1，则y(t0)由向量y直接提取出来，就可以计算Caputo分数阶导数；如果r>1，则初值向量y0实际上的内容y0=[y(t0),y'(t0),y"(t0),...]直到[r]-1阶。

举两个例子观察一下不同定义下求得的分数阶导数：

对于阶跃函数，其0.75阶导数的表达式为：t^(-0.75)/Gama(0.25)。比较G-L定义和R-L定义的微积分的数值解与精确解：

t0 = 0:0.01:3;
y0 = t0.^(-0.75)/gamma(0.25);
h = 0.01;
t1 = 0:h:3;
y = ones(size(t1));
y1 = glfdiff(y,t1,0.75);
[y2,t2] = rlfdiff(y,t1,0.75);
plot(t0,y0,t1,y1,'--',t2,y2,':');
ylim([0 6]);

运行结果：

根据运行结果，G-L定义的分数阶微分求得的数值解更接近于其精确解，这也可能是函数求解的方法的精度的差异，并不能依此说明G-L定义的分数阶积分更准确。

对于函数y=sin(3t+1)，对比其G-L定义和Caputo定义下的0.3阶微分的数值解：

t = 0:0.01:pi;
y = sin(3*t+1);
d = t.^(-0.3)*sin(1)/gamma(0.7);y1=glfdiff(y,t,0.3);
y2=caputo(y,t,0.3,0);plot(t,y1,t,y2,'--',t,d,':')
axis([0 pi -3 3])

运行结果：

由于函数在t0时刻的初值为sin1，故两个定义间的数值解求取出现了偏差。也就是说Caputo定义的分数阶积分在初值不为零时是与其他定义不等价的。

事实上，分数阶微积分求数值解还有更多高精度的算法，本文只给出了比较基础的代码。但是无论什么方法，都很难避免会引入新的误差，分数阶微积分的深度还亟待更加的深入。

MATLAB求分数阶微分的数值解，G-L定义，R-L定义，Caputo定义相关推荐

图像增强-分数阶微分(vc++)
分数阶微分,Grumwald-Letnikov定义在图像的数值实现中更为准确: Gamma函数: 若s(t)的持续期t [a,t],将函数持续期间[a,t]按单位间隔h=1进行等分,得到: 推到一元函 ...
使用matlab求高阶累积量
本文介绍如何使用matlab自带的高阶累积量函数求取一个随机过程的高阶累积量,运行demo之前确保matlab中已经安装了HOSA的工具箱. 安装过程中经常会碰到的两个问题: 1.info.xml出问 ...
利用2阶分数阶微分掩模的边缘检测（Matlab代码实现）
matlab求偏微分方法解析解,偏微分数值解(2,MATLAB求解方法)学案.ppt
这部分主要讨论如何用MATLAB实现对偏微分方程的数值仿真求解．MATLAB的偏微分方程工具箱(PDE Toolbox)的出现,为偏微分方程的求解以及定性研究提供了捷径．主要步骤为: 2.1 用偏微分 ...
matlab中动力学方程,Matlab求动力学的微分方程组拟合
!using["XSLSF"]; //使用命名空间XSLSF //函数定义,用于计算微分方程组中各方程右端函数值,连分式法对微分方程组积分一步函数pb ...
matlab 求高阶累积量,高阶累积量matlab源码
[实例简介] 对数字调制信号ASK.FSK.PSK类利用高阶累积量特征进行识别的matlab程序 [实例截图] [核心代码] 高阶累积量识别信号 └── 高阶累积量 ├── 识别率 │ ├── a ...
Caputo 分数阶一维问题基于 L1 逼近的空间二阶方法(附Matlab代码)
Caputo 分数阶一维问题基于 L1 逼近的空间二阶方法 Caputo 分数阶一维问题基于 L1 逼近的快速差分方法(附Matlab程序) 文章目录 Caputo 分数阶一维问题基于 L1 逼近的空 ...
matlab 分数阶混沌系统的完全同步控制
1.内容简介略 625-可以交流.咨询.答疑 2.内容说明分数阶微积分这一重要的数学分支,其诞生在1695年,几乎和经典微积分同时出现.那一年,德国数学家Leibniz 和法国数学家L'Hopit ...
用matlab求上三角矩阵的逆,现代科学运算—MATLAB语言与应用-中国大学mooc-题库零氪...
第1章绪论 01-01 本课程的主要内容随堂测验 1. 2. 01-02 为什么学习计算机数学语言随堂测验 1. 2. 01-03 解析解与数值解随堂测验 1. A. B. C. D. 2. 01- ...

MATLAB求分数阶微分的数值解，G-L定义，R-L定义，Caputo定义

MATLAB求分数阶微分的数值解，G-L定义，R-L定义，Caputo定义相关推荐

最新文章

热门文章