Matlab caputo fractional code

조회 수: 23 (최근 30일)
noureddine karim
noureddine karim 2023년 4월 20일
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

답변 (0개)

카테고리

Help CenterFile Exchange에서 Programming에 대해 자세히 알아보기

태그

제품


릴리스

R2023a

Community Treasure Hunt

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

Start Hunting!

Translated by