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]);```

Answer 1

[Jacob Mathew]

Hey Naman,

I understand that you tried to plot a system of fractional order differential equations. To simplify solving and plotting fractional order differential equation, you can install and use the “FOTF Toolbox” Add On. You can find the Add On at the link below:

https://www.mathworks.com/matlabcentral/fileexchange/60874-fotf-toolbox

Refer to the following MATLAB Answer on how to install toolboxes:

https://www.mathworks.com/matlabcentral/answers/101885-how-do-i-install-additional-toolboxes-into-an-existing-installation-of-matlab

We will specifically use the “fode_sol” function from the toolbox. Its implementation can be found here:

https://viewer.mathworks.com/?viewer=plain_code&url=https%3A%2F%2Fwww.mathworks.com%2Fmatlabcentral%2Fmlc-downloads%2Fdownloads%2Fe34b8d0e-decc-4329-9417-d5735ec524b9%2F4b5a4f0f-b0db-4286-be15-ef698df894ae%2Ffiles%2Ffode_sol.m&embed=web

An example of using the function to plot a system of fractional order equations is as follows:

% Define the coefficients and orders for the first equation
a1 = [1, 0.5];         % Coefficients of the numerator for x(t)
na1 = [0.9, 0.9];      % Orders of the numerator (fractional order)
b1 = [1, -0.5];        % Coefficients of the denominator for x(t)
nb1 = [0, 0.9];        % Orders of the denominator (fractional order)

% Define the coefficients and orders for the second equation
a2 = [0.5, 1];         % Coefficients of the numerator for y(t)
na2 = [0.8, 0.8];      % Orders of the numerator (fractional order)
b2 = [1, -0.5];        % Coefficients of the denominator for y(t)
nb2 = [0, 0.8];        % Orders of the denominator (fractional order)

% Define the time vector
t = linspace(0, 10, 1000)';

% Define the input signals (e.g., step functions)
u1 = ones(size(t));
u2 = ones(size(t));

% Solve the fractional-order differential equations using fode_sol
x = fode_sol(a1, na1, b1, nb1, u1, t);
y = fode_sol(a2, na2, b2, nb2, u2, t);

% Plot the results
figure;
plot(t, x, 'r', 'DisplayName', 'x(t)');
hold on;
plot(t, y, 'b', 'DisplayName', 'y(t)');
xlabel('Time t');
ylabel('State variables');
legend;
title('Solution of System of Fractional-Order Differential Equations');
grid on;

Depending on your requirements, you can refer to other functions available in the FOTF Toolbox

[Naman]

@Jacob.—Thankyou for the response. I see these links, I will check and try them out. I had found fde12.m file in some other discussion, and that has helped me make some plots. Certainly, perhaps this toolbox might give better results if it uses other numerical techniques than what fde12.m does. I will try it out. Thankyou...