如何用python计算和表示这两个公式?
第一个是fractional integrals (分数阶积分)
第二个是fractional derivatives(分数导数)
同学,你好,我是matlab专栏下的专家。你的这种算法在matlab上早已实现,而且代码也很简单,你可以看看:
Fractional Derivative and Integral - File Exchange - MATLAB Central
https://www.mathworks.com/matlabcentral/fileexchange/52587-fractional-derivative-and-integral
可以把代码放在这里
积分
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}'])