草庐IT

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

clear_lantern 2023-12-19 原文

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

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

h = 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==-1
    n = length(t);
    yy = 0.5*(y(1:n-1)+y(2:n)).*diff(t);
    dy(1) = 0;
    for i=2:n
        dy(i) = dy(i-1)+yy(i-1);
    end
    return
else
    m=-r;
    y3 = t.^(m-1);
end

for i = 1:n
    dy1(i) = y(i:-1:1).'*(y3(1:i));
end

if r>-1
    dy = diff(dy1,m)/(h^(m-1))/gamma(p);
    t = t(1:n-m);
else
    dy = 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<=4
    L = 10;
end
if alfa>0
    q = ceil(alfa);
    if alfa<=1
        y0 = y(1);
    end
    for k = 0:q-1
        dy = dy-y0(k+1)*t.^(k-alfa)./gamma(k+1-alfa);
    end
    yy1 = 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定义的更多相关文章

  1. ruby - Facter::Util::Uptime:Module 的未定义方法 get_uptime (NoMethodError) - 2

    我正在尝试设置一个puppet节点,但ruby​​gems似乎不正常。如果我通过它自己的二进制文件(/usr/lib/ruby/gems/1.8/gems/facter-1.5.8/bin/facter)在cli上运行facter,它工作正常,但如果我通过由ruby​​gems(/usr/bin/facter)安装的二进制文件,它抛出:/usr/lib/ruby/1.8/facter/uptime.rb:11:undefinedmethod`get_uptime'forFacter::Util::Uptime:Module(NoMethodError)from/usr/lib/ruby

  2. ruby-on-rails - 如果为空或不验证数值,则使属性默认为 0 - 2

    我希望我的UserPrice模型的属性在它们为空或不验证数值时默认为0。这些属性是tax_rate、shipping_cost和price。classCreateUserPrices8,:scale=>2t.decimal:tax_rate,:precision=>8,:scale=>2t.decimal:shipping_cost,:precision=>8,:scale=>2endendend起初,我将所有3列的:default=>0放在表格中,但我不想要这样,因为它已经填充了字段,我想使用占位符。这是我的UserPrice模型:classUserPrice回答before_val

  3. ruby-on-rails - Rails 3.2.1 中 ActionMailer 中的未定义方法 'default_content_type=' - 2

    我在我的项目中添加了一个系统来重置用户密码并通过电子邮件将密码发送给他,以防他忘记密码。昨天它运行良好(当我实现它时)。当我今天尝试启动服务器时,出现以下错误。=>BootingWEBrick=>Rails3.2.1applicationstartingindevelopmentonhttp://0.0.0.0:3000=>Callwith-dtodetach=>Ctrl-CtoshutdownserverExiting/Users/vinayshenoy/.rvm/gems/ruby-1.9.3-p0/gems/actionmailer-3.2.1/lib/action_mailer

  4. ruby-on-rails - form_for 中不在模型中的自定义字段 - 2

    我想向我的Controller传递一个参数,它是一个简单的复选框,但我不知道如何在模型的form_for中引入它,这是我的观点:{:id=>'go_finance'}do|f|%>Transferirde:para:Entrada:"input",:placeholder=>"Quantofoiganho?"%>Saída:"output",:placeholder=>"Quantofoigasto?"%>Nota:我想做一个额外的复选框,但我该怎么做,模型中没有一个对象,而是一个要检查的对象,以便在Controller中创建一个ifelse,如果没有检查,请帮助我,非常感谢,谢谢

  5. ruby - 主要 :Object when running build from sublime 的未定义方法 `require_relative' - 2

    我已经从我的命令行中获得了一切,所以我可以运行rubymyfile并且它可以正常工作。但是当我尝试从sublime中运行它时,我得到了undefinedmethod`require_relative'formain:Object有人知道我的sublime设置中缺少什么吗?我正在使用OSX并安装了rvm。 最佳答案 或者,您可以只使用“require”,它应该可以正常工作。我认为“require_relative”仅适用于ruby​​1.9+ 关于ruby-主要:Objectwhenrun

  6. ruby - 在 Ruby 中有条件地定义函数 - 2

    我有一些代码在几个不同的位置之一运行:作为具有调试输出的命令行工具,作为不接受任何输出的更大程序的一部分,以及在Rails环境中。有时我需要根据代码的位置对代码进行细微的更改,我意识到以下样式似乎可行:print"Testingnestedfunctionsdefined\n"CLI=trueifCLIdeftest_printprint"CommandLineVersion\n"endelsedeftest_printprint"ReleaseVersion\n"endendtest_print()这导致:TestingnestedfunctionsdefinedCommandLin

  7. ruby - 定义方法参数的条件 - 2

    我有一个只接受一个参数的方法:defmy_method(number)end如果使用number调用方法,我该如何引发错误??通常,我如何定义方法参数的条件?比如我想在调用的时候报错:my_method(1) 最佳答案 您可以添加guard在函数的开头,如果参数无效则引发异常。例如:defmy_method(number)failArgumentError,"Inputshouldbegreaterthanorequalto2"ifnumbereputse.messageend#=>Inputshouldbegreaterthano

  8. ruby - 如何在 Grape 中定义哈希数组? - 2

    我使用Ember作为我的前端和GrapeAPI来为我的API提供服务。前端发送类似:{"service"=>{"name"=>"Name","duration"=>"30","user"=>nil,"organization"=>"org","category"=>nil,"description"=>"description","disabled"=>true,"color"=>nil,"availabilities"=>[{"day"=>"Saturday","enabled"=>false,"timeSlots"=>[{"startAt"=>"09:00AM","endAt"=>

  9. ruby - 获取模块中定义的所有常量的值 - 2

    我想获取模块中定义的所有常量的值:moduleLettersA='apple'.freezeB='boy'.freezeendconstants给了我常量的名字:Letters.constants(false)#=>[:A,:B]如何获取它们的值的数组,即["apple","boy"]? 最佳答案 为了做到这一点,请使用mapLetters.constants(false).map&Letters.method(:const_get)这将返回["a","b"]第二种方式:Letters.constants(false).map{|c

  10. ruby - 这两个 Ruby 类初始化定义有什么区别? - 2

    我正在阅读一本关于Ruby的书,作者在编写类初始化定义时使用的形式与他在本书前几节中使用的形式略有不同。它看起来像这样:classTicketattr_accessor:venue,:datedefinitialize(venue,date)self.venue=venueself.date=dateendend在本书的前几节中,它的定义如下:classTicketattr_accessor:venue,:datedefinitialize(venue,date)@venue=venue@date=dateendend在第一个示例中使用setter方法与在第二个示例中使用实例变量之间是

随机推荐