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

zhangguangli

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

固有频率

下面的程序仅作参考

1. % This function solves linear fractional-order differential equation (fode)
   
2. % with constant coefficients of the form:
   
3. %
   
4. %        n/v            (n-1)/v                1/v
   
5. % c1(1)*D   y(t)+c1(2)*D       y(t)+...+c1(n)*D   y(t)+c1(n+1)*y(t) =
   
6. %
   
7. %                m/v            (m-1)/v                1/v           
   
8. %         c2(1)*D   r(t)+c2(2)*D       r(t)+...+c2(m)*D   r(t)+c2(m+1)*r(t)      
   
9. %
   
10. %
    
11. % Inputs:
    
12. %
    
13. %           v  : common denominator
    
14. %           c1 : vector of output coefficients (1x(n+1))
    
15. %           c2 : vector of input coefficients (1x(m+1))
    
16. %           r  : samples of the input signal r(t)
    
17. %           h  : sampling period
    
18. %
    
19. % Outputs:
    
20. %           y  : vector of output samples
    
21. %           t  : time vector (corresponding to y)
    
22. %
    
23. % Author: Farshad Merrikh Bayat ([email]fbayat@ee.sharif.edu[/email])
    
24. % URL: http:// ee.sharif.edu/~fbayat/
    
25. %
    
26. % Note: this function calls the function "fderiv.m" which is also
    
27. %       downloadable from MathWorks-File Exchange. The parameter "h" can
    
28. %       easily be tuned; it must be as small as approximate the input signal
    
29. %       r(t). The short memory principle has not neen used here, so the
    
30. %       length of input signal is limited to few hundred samples.
    
31. %
    
32. % Copyright (c), 2007.
    
33. %
    
34. 
    
35. 
    
36. 
    
37. function [t,y] = fode2(v,c1,c2,r,h)
    
38. 
    
39. n = length(c1)-1;
    
40. % r = k*ones(1,100); % if the input signal is unit step
    
41. temp1 = zeros(size(r));
    
42. for i=1:length(c2)
    
43.     r_new = fderiv((length(c2)-i)/v,r,h);
    
44.     temp1 = temp1+c2(i)*r_new;
    
45. end
    
46. 
    
47. r = temp1;
    
48. 
    
49. t = [0:1:length(r)-1]*h;
    
50. y = zeros(1,length(r));
    
51. 
    
52. temp = zeros(1,n);
    
53. 
    
54. a = 0;
    
55. for i=1:length(c1)
    
56.     a = a+c1(i)/h^((n-i+1)/v);
    
57. end
    
58. 
    
59. for i=1:length(r)
    
60.     for k=n:-1:1
    
61.         for j=1:i-1
    
62.             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));
    
63.         end
    
64.         temp(n-k+1) = -c1(n-k+1)*temp(n-k+1)/h^(k/v);
    
65.     end
    
66.     y(i) = (sum(temp)+r(i))/a;
    
67.     temp = zeros(1,n);
    
68. end
    
69. y = [0 y(1:length(y)-1)];
    
70. %plot(t,y)

复制代码