im looking for the error in the following code, please help

% Fractional SIR System using Caputo Derivative

clear all; clc;

% Parameters
N = 1000;       % Population size
beta = 0.5;     % Contact rate
gamma = 0.1;    % Recovery rate
alpha = 0.7;    % Fractional order

% Initial conditions
I0 = 1;         % Initial number of infected individuals
S0 = N - I0;    % Initial number of susceptible individuals
R0 = 0;         % Initial number of recovered individuals

% Time steps
t_start = 0;
t_end = 50;
dt = 0.1;
t = t_start:dt:t_end;

% Caputo derivative
CaputoDerivative = @(f,alpha,dt) 1/gamma(alpha)*(f(1)-f(2))/dt^(alpha);

% Initialization
S = zeros(1,length(t));
I = zeros(1,length(t));
R = zeros(1,length(t));

S(1) = S0;
I(1) = I0;
R(1) = R0;

% Iterative solution
for n = 1:length(t)-1
    DalphaS = CaputoDerivative(S(1:n),alpha,dt);
    DalphaI = CaputoDerivative(I(1:n),alpha,dt);
    DalphaR = CaputoDerivative(R(1:n),alpha,dt);
    
    S(n+1) = S(n) - beta*S(n)*DalphaI*dt;
    I(n+1) = I(n) + beta*S(n)*DalphaI*dt - gamma*I(n)*DalphaI*dt;
    R(n+1) = R(n) + gamma*I(n)*DalphaR*dt;
end

Array indices must be positive integers or logical values.

Error in solution>@(f,alpha,dt)1/gamma(alpha)*(f(1)-f(2))/dt^(alpha) (line 23)
CaputoDerivative = @(f,alpha,dt) 1/gamma(alpha)*(f(1)-f(2))/dt^(alpha);

[noureddine karim]

I change the code to this following, but same probleme

% Parameters
alpha = 0.5; % Fractional order
beta = 0.1; % Infection rate
gamma = 0.05; % Recovery rate
N = 10000; % Total population
I0 = 10; % Initial number of infected individuals
R0 = 0; % Initial number of recovered individuals
S0 = N - I0 - R0; % Initial number of susceptible individuals
tend = 100; % End time
dt = 0.1; % Time step

% Initialization
t = 0:dt:tend;
S = zeros(size(t));
I = zeros(size(t));
R = zeros(size(t));
S(1) = S0;
I(1) = I0;
R(1) = R0;

% Fractional differentiation operator
Dalpha = @(t) (1./gamma(alpha)).*(t^(1-alpha)).*diff([0,t.^(alpha-1)]);

% SIR model simulation
for n = 1:length(t)-1
    S(n+1) = S(n) - beta*S(n)*Dalpha(I(n));
    I(n+1) = I(n) + beta*S(n)*Dalpha(I(n)) - gamma*I(n);
    R(n+1) = R(n) + gamma*I(n);
end

% Plotting
plot(t,S,'b',t,I,'r',t,R,'g')
legend('S','I','R')
xlabel('Time')
ylabel('Number of individuals')
title('SIR Fractional Model Simulation')

[Torsten]

I change the code to this following, but same probleme

But you still didn't change the

gamma(alpha)

in the definition of the fractional differentiation operator. It throws an error if you previously defined a scalar gamma as

gamma = 0.05; % Recovery rate

because the gamma function will not be used then.

Further

diff([0,t.^(alpha-1)]) 

equals

t.^(alpha-1) - 0 = t.^(alpha-1)

If this is what you want, there is no need to use "diff" here.

[noureddine karim]

Thanks man

Answer 1

[Image Analyst]

See the FAQ for a thorough discussion:

https://matlab.fandom.com/wiki/FAQ#%22Subscript_indices_must_either_be_real_positive_integers_or_logicals.%22

Answer 2

[Torsten]

DalphaS = CaputoDerivative(S(1:n),alpha,dt);
DalphaI = CaputoDerivative(I(1:n),alpha,dt);
DalphaR = CaputoDerivative(R(1:n),alpha,dt);

For n = 1, you call CaputoDerivative with S(1:n), I(1:n) and R(1:n) each having only a single element.

But CaputoDerivative expects that it has at least two elements: f(1) and f(2).

And - as Walter noticed - gamma(alpha) is not defined in your case because you set gamma to a value instead of using the gamma function in your formula.

Answer 3

[Walter Roberson]

gamma = 0.1;    % Recovery rate

gamma is a scalar not a function

CaputoDerivative = @(f,alpha,dt) 1/gamma(alpha)*(f(1)-f(2))/dt^(alpha);

That tries to use the gamma function, but gamma is a scalar

[Walter Roberson]

If the intention is then you would need 1./(gamma .* alpha) . (Note that could be calculated in advance instead of recalculating it each time.)