Matlab caputo fractional code

23 Ansichten (letzte 30 Tage)
noureddine karim
noureddine karim am 20 Apr. 2023
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

Antworten (0)

Kategorien

Mehr zu Programming finden Sie in Help Center und File Exchange

Tags

Produkte


Version

R2023a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by