zhangguangli 发表于 2011-1-6 17:06

分数阶常微分方程的求解matlab程序

求 求解分数阶常微分方程的matlab程序代码,谢谢

固有频率 发表于 2011-1-6 22:13

下面的程序仅作参考% This function solves linear fractional-order differential equation (fode)
% with constant coefficients of the form:
%
%      n/v            (n-1)/v                1/v
% c1(1)*D   y(t)+c1(2)*D       y(t)+...+c1(n)*D   y(t)+c1(n+1)*y(t) =
%
%                m/v            (m-1)/v                1/v         
%         c2(1)*D   r(t)+c2(2)*D       r(t)+...+c2(m)*D   r(t)+c2(m+1)*r(t)      
%
%
% Inputs:
%
%         v: common denominator
%         c1 : vector of output coefficients (1x(n+1))
%         c2 : vector of input coefficients (1x(m+1))
%         r: samples of the input signal r(t)
%         h: sampling period
%
% Outputs:
%         y: vector of output samples
%         t: time vector (corresponding to y)
%
% Author: Farshad Merrikh Bayat (fbayat@ee.sharif.edu)
% URL: http:// ee.sharif.edu/~fbayat/
%
% Note: this function calls the function "fderiv.m" which is also
%       downloadable from MathWorks-File Exchange. The parameter "h" can
%       easily be tuned; it must be as small as approximate the input signal
%       r(t). The short memory principle has not neen used here, so the
%       length of input signal is limited to few hundred samples.
%
% Copyright (c), 2007.
%



function = fode2(v,c1,c2,r,h)

n = length(c1)-1;
% r = k*ones(1,100); % if the input signal is unit step
temp1 = zeros(size(r));
for i=1:length(c2)
    r_new = fderiv((length(c2)-i)/v,r,h);
    temp1 = temp1+c2(i)*r_new;
end

r = temp1;

t = *h;
y = zeros(1,length(r));

temp = zeros(1,n);

a = 0;
for i=1:length(c1)
    a = a+c1(i)/h^((n-i+1)/v);
end

for i=1:length(r)
    for k=n:-1:1
      for j=1:i-1
            temp(n-k+1) = temp(n-k+1)+(-1)^j*gamma(k/v+1)*y(i-j)/(gamma(j+1)*gamma(k/v-j+1));
      end
      temp(n-k+1) = -c1(n-k+1)*temp(n-k+1)/h^(k/v);
    end
    y(i) = (sum(temp)+r(i))/a;
    temp = zeros(1,n);
end
y = ;
%plot(t,y)
页: [1]
查看完整版本: 分数阶常微分方程的求解matlab程序