Matlab caputo fractional code

Question

Help please, im in need for a code matlab fractional i found this code but it dosnt work, if you have in idea please chare with me 

% SIR fractional system using Caputo fractional derivative

clear all; close all; clc;

% Parameters
alpha = 0.8;     % order of the fractional derivative
beta = 0.4;      % infection rate
gamma = 0.2;     % recovery rate
N = 1000;        % total population
I0 = 1;          % initial number of infected individuals
R0 = 0;          % initial number of recovered individuals
S0 = N - I0 - R0;% initial number of susceptible individuals
tend = 10;       % end time
dt = 0.01;       % time step

% Initialization
t = 0:dt:tend;
Nt = length(t);
S = zeros(1, Nt);
I = zeros(1, Nt);
R = zeros(1, Nt);
S(1) = S0;
I(1) = I0;
R(1) = R0;

% Main loop
for n = 1:Nt-1
    % Compute fractional derivatives using the Caputo derivative
    DalphaS = CaputoDerivative(S(n),alpha,t(n),dt);
    DalphaI = CaputoDerivative(I(n),alpha,t(n),dt);
    DalphaR = CaputoDerivative(R(n),alpha,t(n),dt);
    
    % Compute derivatives using standard differential equations
    dSdt = -beta*S(n)*DalphaI/N;
    dIdt = beta*S(n)*DalphaI/N - gamma*DalphaI;
    dRdt = gamma*DalphaI;
    
    % Update variables
    S(n+1) = S(n) + dt*dSdt;
    I(n+1) = I(n) + dt*dIdt;
    R(n+1) = R(n) + dt*dRdt;
end

% Plot results
plot(t, S, 'r', t, I, 'g', t, R, 'b');
legend('Susceptible', 'Infected', 'Recovered');
xlabel('Time');
ylabel('Population');
title('Fractional SIR System using Caputo Derivative');
grid on;

% Caputo derivative function
function y = CaputoDerivative(f,alpha,t,dt)
    % Computes the Caputo fractional derivative of order alpha
    % f: function to differentiate
    % alpha: order of the fractional derivative
    % t: current time
    % dt: time step
    
    % Compute the intermediate values for the Caputo derivative
    N = length(f);
    m = floor(alpha);
    A = zeros(N,N-m);
    for k = m:N-1
        A(k+1-m,k+1-m:k) = dt.^(-(0:k-m)).*gamma(k-m+1)./gamma(k+2-m).*(-1).^(0:k-m);
    end
    
    % Compute the Caputo derivative
    y = A*f(m+1:N)';
end

Matlab caputo fractional code

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

Matlab caputo fractional code

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论