Plotting a system of fractional order differential equations
Hi. I have been trying to plot fractional order differential equations with MATLAB. I have tried to use this technique called "Predictor-Corrector Method" and got a code listing using AI. The code seems to be giving the plots, though they do not match with the expected results. I will attach the code listing below. But my question is how to actually plot such systems? Is there a standard method in-built, or a widely accepted add-on/package? Thanks. “`% Parameters
alpha = 0.9;
LAMBDA = 3;
mu = 0.4;
beta = 0.004;
gamma = 0.0005;
m = 0.03;
n = 0.006;
delta = 0.02;
rho = 0.003;
h = 0.1; % Time step
T = 300; % Total time
% Time steps
time = 0:h:T;
time_steps = length(time);
% Initial conditions
s = zeros(1, time_steps);
i = zeros(1, time_steps);
r = zeros(1, time_steps);
s(1) = 85;
i(1) = 3;
r(1) = 0;
% Predictor-Corrector Method
for k = 1:(time_steps – 1)
% Predictor step (Adams-Bashforth)
if k == 1
ds_pred = h * (LAMBDA – mu * s(k) – (beta * s(k) * i(k)) / (1 + gamma * i(k)));
di_pred = h * ((beta * s(k) * i(k)) / (1 + gamma * i(k)) – (m * i(k)) / (1 + n * i(k)) – (mu + delta + rho) * i(k));
dr_pred = h * ((m * i(k)) / (1 + n * i(k)) + delta * i(k) – mu * r(k));
else
ds_pred = h / 2 * ((LAMBDA – mu * s(k) – (beta * s(k) * i(k)) / (1 + gamma * i(k))) + …
(LAMBDA – mu * s(k-1) – (beta * s(k-1) * i(k-1)) / (1 + gamma * i(k-1))));
di_pred = h / 2 * (((beta * s(k) * i(k)) / (1 + gamma * i(k)) – (m * i(k)) / (1 + n * i(k)) – (mu + delta + rho) * i(k)) + …
((beta * s(k-1) * i(k-1)) / (1 + gamma * i(k-1)) – (m * i(k-1)) / (1 + n * i(k-1)) – (mu + delta + rho) * i(k-1)));
dr_pred = h / 2 * (((m * i(k)) / (1 + n * i(k)) + delta * i(k) – mu * r(k)) + …
((m * i(k-1)) / (1 + n * i(k-1)) + delta * i(k-1) – mu * r(k-1)));
end
% Update predicted values
s(k + 1) = s(k) + ds_pred;
i(k + 1) = i(k) + di_pred;
r(k + 1) = r(k) + dr_pred;
% Corrector step (Adams-Moulton)
ds_corr = h / 2 * ((LAMBDA – mu * s(k + 1) – (beta * s(k + 1) * i(k + 1)) / (1 + gamma * i(k + 1))) + …
(LAMBDA – mu * s(k) – (beta * s(k) * i(k)) / (1 + gamma * i(k))));
di_corr = h / 2 * (((beta * s(k + 1) * i(k + 1)) / (1 + gamma * i(k + 1)) – (m * i(k + 1)) / (1 + n * i(k + 1)) – (mu + delta + rho) * i(k + 1)) + …
((beta * s(k) * i(k)) / (1 + gamma * i(k)) – (m * i(k)) / (1 + n * i(k)) – (mu + delta + rho) * i(k)));
dr_corr = h / 2 * (((m * i(k + 1)) / (1 + n * i(k + 1)) + delta * i(k + 1) – mu * r(k + 1)) + …
((m * i(k)) / (1 + n * i(k)) + delta * i(k) – mu * r(k)));
% Update corrected values
s(k + 1) = s(k) + ds_corr;
i(k + 1) = i(k) + di_corr;
r(k + 1) = r(k) + dr_corr;
end
% Plot the results
figure;
plot(time, s, ‘-o’, ‘DisplayName’, ‘s(t)’);
hold on;
plot(time, i, ‘-s’, ‘DisplayName’, ‘i(t)’);
plot(time, r, ‘-d’, ‘DisplayName’, ‘r(t)’);
xlabel(‘Time’);
ylabel(‘Population’);
title(‘Solution of the System of Fractional Differential Equations’);
legend;
grid on;
% Set plot range
xlim([0 30]);
ylim([0 100]);“`
<</matlabcentral/answers/uploaded_files/1745286/fo_M.png>>Hi. I have been trying to plot fractional order differential equations with MATLAB. I have tried to use this technique called "Predictor-Corrector Method" and got a code listing using AI. The code seems to be giving the plots, though they do not match with the expected results. I will attach the code listing below. But my question is how to actually plot such systems? Is there a standard method in-built, or a widely accepted add-on/package? Thanks. “`% Parameters
alpha = 0.9;
LAMBDA = 3;
mu = 0.4;
beta = 0.004;
gamma = 0.0005;
m = 0.03;
n = 0.006;
delta = 0.02;
rho = 0.003;
h = 0.1; % Time step
T = 300; % Total time
% Time steps
time = 0:h:T;
time_steps = length(time);
% Initial conditions
s = zeros(1, time_steps);
i = zeros(1, time_steps);
r = zeros(1, time_steps);
s(1) = 85;
i(1) = 3;
r(1) = 0;
% Predictor-Corrector Method
for k = 1:(time_steps – 1)
% Predictor step (Adams-Bashforth)
if k == 1
ds_pred = h * (LAMBDA – mu * s(k) – (beta * s(k) * i(k)) / (1 + gamma * i(k)));
di_pred = h * ((beta * s(k) * i(k)) / (1 + gamma * i(k)) – (m * i(k)) / (1 + n * i(k)) – (mu + delta + rho) * i(k));
dr_pred = h * ((m * i(k)) / (1 + n * i(k)) + delta * i(k) – mu * r(k));
else
ds_pred = h / 2 * ((LAMBDA – mu * s(k) – (beta * s(k) * i(k)) / (1 + gamma * i(k))) + …
(LAMBDA – mu * s(k-1) – (beta * s(k-1) * i(k-1)) / (1 + gamma * i(k-1))));
di_pred = h / 2 * (((beta * s(k) * i(k)) / (1 + gamma * i(k)) – (m * i(k)) / (1 + n * i(k)) – (mu + delta + rho) * i(k)) + …
((beta * s(k-1) * i(k-1)) / (1 + gamma * i(k-1)) – (m * i(k-1)) / (1 + n * i(k-1)) – (mu + delta + rho) * i(k-1)));
dr_pred = h / 2 * (((m * i(k)) / (1 + n * i(k)) + delta * i(k) – mu * r(k)) + …
((m * i(k-1)) / (1 + n * i(k-1)) + delta * i(k-1) – mu * r(k-1)));
end
% Update predicted values
s(k + 1) = s(k) + ds_pred;
i(k + 1) = i(k) + di_pred;
r(k + 1) = r(k) + dr_pred;
% Corrector step (Adams-Moulton)
ds_corr = h / 2 * ((LAMBDA – mu * s(k + 1) – (beta * s(k + 1) * i(k + 1)) / (1 + gamma * i(k + 1))) + …
(LAMBDA – mu * s(k) – (beta * s(k) * i(k)) / (1 + gamma * i(k))));
di_corr = h / 2 * (((beta * s(k + 1) * i(k + 1)) / (1 + gamma * i(k + 1)) – (m * i(k + 1)) / (1 + n * i(k + 1)) – (mu + delta + rho) * i(k + 1)) + …
((beta * s(k) * i(k)) / (1 + gamma * i(k)) – (m * i(k)) / (1 + n * i(k)) – (mu + delta + rho) * i(k)));
dr_corr = h / 2 * (((m * i(k + 1)) / (1 + n * i(k + 1)) + delta * i(k + 1) – mu * r(k + 1)) + …
((m * i(k)) / (1 + n * i(k)) + delta * i(k) – mu * r(k)));
% Update corrected values
s(k + 1) = s(k) + ds_corr;
i(k + 1) = i(k) + di_corr;
r(k + 1) = r(k) + dr_corr;
end
% Plot the results
figure;
plot(time, s, ‘-o’, ‘DisplayName’, ‘s(t)’);
hold on;
plot(time, i, ‘-s’, ‘DisplayName’, ‘i(t)’);
plot(time, r, ‘-d’, ‘DisplayName’, ‘r(t)’);
xlabel(‘Time’);
ylabel(‘Population’);
title(‘Solution of the System of Fractional Differential Equations’);
legend;
grid on;
% Set plot range
xlim([0 30]);
ylim([0 100]);“`
<</matlabcentral/answers/uploaded_files/1745286/fo_M.png>> Hi. I have been trying to plot fractional order differential equations with MATLAB. I have tried to use this technique called "Predictor-Corrector Method" and got a code listing using AI. The code seems to be giving the plots, though they do not match with the expected results. I will attach the code listing below. But my question is how to actually plot such systems? Is there a standard method in-built, or a widely accepted add-on/package? Thanks. “`% Parameters
alpha = 0.9;
LAMBDA = 3;
mu = 0.4;
beta = 0.004;
gamma = 0.0005;
m = 0.03;
n = 0.006;
delta = 0.02;
rho = 0.003;
h = 0.1; % Time step
T = 300; % Total time
% Time steps
time = 0:h:T;
time_steps = length(time);
% Initial conditions
s = zeros(1, time_steps);
i = zeros(1, time_steps);
r = zeros(1, time_steps);
s(1) = 85;
i(1) = 3;
r(1) = 0;
% Predictor-Corrector Method
for k = 1:(time_steps – 1)
% Predictor step (Adams-Bashforth)
if k == 1
ds_pred = h * (LAMBDA – mu * s(k) – (beta * s(k) * i(k)) / (1 + gamma * i(k)));
di_pred = h * ((beta * s(k) * i(k)) / (1 + gamma * i(k)) – (m * i(k)) / (1 + n * i(k)) – (mu + delta + rho) * i(k));
dr_pred = h * ((m * i(k)) / (1 + n * i(k)) + delta * i(k) – mu * r(k));
else
ds_pred = h / 2 * ((LAMBDA – mu * s(k) – (beta * s(k) * i(k)) / (1 + gamma * i(k))) + …
(LAMBDA – mu * s(k-1) – (beta * s(k-1) * i(k-1)) / (1 + gamma * i(k-1))));
di_pred = h / 2 * (((beta * s(k) * i(k)) / (1 + gamma * i(k)) – (m * i(k)) / (1 + n * i(k)) – (mu + delta + rho) * i(k)) + …
((beta * s(k-1) * i(k-1)) / (1 + gamma * i(k-1)) – (m * i(k-1)) / (1 + n * i(k-1)) – (mu + delta + rho) * i(k-1)));
dr_pred = h / 2 * (((m * i(k)) / (1 + n * i(k)) + delta * i(k) – mu * r(k)) + …
((m * i(k-1)) / (1 + n * i(k-1)) + delta * i(k-1) – mu * r(k-1)));
end
% Update predicted values
s(k + 1) = s(k) + ds_pred;
i(k + 1) = i(k) + di_pred;
r(k + 1) = r(k) + dr_pred;
% Corrector step (Adams-Moulton)
ds_corr = h / 2 * ((LAMBDA – mu * s(k + 1) – (beta * s(k + 1) * i(k + 1)) / (1 + gamma * i(k + 1))) + …
(LAMBDA – mu * s(k) – (beta * s(k) * i(k)) / (1 + gamma * i(k))));
di_corr = h / 2 * (((beta * s(k + 1) * i(k + 1)) / (1 + gamma * i(k + 1)) – (m * i(k + 1)) / (1 + n * i(k + 1)) – (mu + delta + rho) * i(k + 1)) + …
((beta * s(k) * i(k)) / (1 + gamma * i(k)) – (m * i(k)) / (1 + n * i(k)) – (mu + delta + rho) * i(k)));
dr_corr = h / 2 * (((m * i(k + 1)) / (1 + n * i(k + 1)) + delta * i(k + 1) – mu * r(k + 1)) + …
((m * i(k)) / (1 + n * i(k)) + delta * i(k) – mu * r(k)));
% Update corrected values
s(k + 1) = s(k) + ds_corr;
i(k + 1) = i(k) + di_corr;
r(k + 1) = r(k) + dr_corr;
end
% Plot the results
figure;
plot(time, s, ‘-o’, ‘DisplayName’, ‘s(t)’);
hold on;
plot(time, i, ‘-s’, ‘DisplayName’, ‘i(t)’);
plot(time, r, ‘-d’, ‘DisplayName’, ‘r(t)’);
xlabel(‘Time’);
ylabel(‘Population’);
title(‘Solution of the System of Fractional Differential Equations’);
legend;
grid on;
% Set plot range
xlim([0 30]);
ylim([0 100]);“`
<</matlabcentral/answers/uploaded_files/1745286/fo_M.png>> fde MATLAB Answers — New Questions
