如何使用python计算fractional integrals and derivative

img

img


如何用python计算和表示这两个公式？
第一个是fractional integrals （分数阶积分）
第二个是fractional derivatives（分数导数）

同学，你好，我是matlab专栏下的专家。你的这种算法在matlab上早已实现，而且代码也很简单，你可以看看：

可以把代码放在这里
积分

function output = fint(x, alpha, n, h)
% FINT Calculate the integral of x
%    x      Sampled function
%    alpha  Fractional order
%    n      Size of the memory
%    h      Sampling rate
%    Please refer to http://www3.nd.edu/~msen/Teaching/UnderRes/FracCalc1.pdf
%    See equ. 36
%    Read section 4 p. 19: since gamma(z) where z > 172, MATLAB gives inf
%    Therefore, the size of memory is limited by  0 < n < 171
N = length(x);
output = zeros(1, N);
h_a = h^alpha;
for i = 1 : N
    sigma = 0;
    for m = 0 : 1 : n
    
    if i - m < 1
        break;
    end
        
        sigma = sigma +  gamma(alpha + m) / (factorial(m) * ...
                gamma(alpha)) * x(i - m);
    end
    output(i) = sigma * h_a;
end
end

微分

function output = fderiv(x, alpha, n, h)
% FDERIVE Calculate the derivative of x
%    x      Sampled function
%    alpha  Fractional order
%    n      Size of the memory
%    h      Sampling rate
%
%    Please refer to http://www3.nd.edu/~msen/Teaching/UnderRes/FracCalc1.pdf
%    See equ. 32
N = length(x);
output = zeros(1, N);
h_a = h^alpha;
for i = 1 : N
    sigma = 0;
    for m = 0 : 1 : n
    
    if i - m < 1
        break;
    end
        
        sigma = sigma + (-1)^m * gamma(alpha + 1) / (factorial(m) * ...
                gamma(alpha - m + 1)) * x(i - m);
    end
    
    output(i) = sigma / h_a;
end
end

计算微分和积分

ts = 0.01;
n = 170;
t = 0:ts:2;
x = t;
figure;
hold on;
plot(t, x)
x_d = fderiv(x, 1, n, ts);
x_i = fint(x, 1,   n, ts);
plot(t, x_d)
plot(t, x_i)
legend('x', 'xd', 'xint')
xlabel('t');
ylabel('x(t)');

当然转化成python也很容易
转化后的例子

import numpy as np
from math import gamma, factorial
def fint(x, alpha, n, h):
    N = len(x)
    output = np.zeros(N)
    h_a = h**alpha
    for i in range(N):
        sigma = 0
        for m in range(n+1):
            if(i-m<0):
                break 
            sigma = sigma + gamma(alpha + m) / (factorial(m)*gamma(alpha)) * x[i-m]
        output[i] = sigma*h_a
    return output

def fderiv(x, alpha, n, h):
    N = len(x);
    output = np.zeros(N)
    h_a = h**alpha
    for i in range(N):
        sigma = 0
        for m in range(n+1):
            if (i - m < 0 or alpha - m + 1<=0):
                break
            sigma = sigma + ((-1)**m) * gamma(alpha + 1) / (factorial(m) * gamma(alpha - m + 1)) * x[i-m]
        output[i] = sigma / h_a
    return output

ts = 0.01
n = 170
t = np.arange(0,2,ts)
x = t

from matplotlib import pyplot as plt

x_i = fint(x, 1,  n, ts);
x_d = fderiv(x, 1,  n, ts);
plt.plot(t, x, 'r-', t, x_i, 'b--', t, x_d,'c-.')
plt.legend(['x', 'x_{int}', 'x_{deriv}'])

img