Why am I not getting a periodic plot?

So, this is my code:

mainfunc()

function mainfunc()

%Parameters

D0 = 0.21829;

betaP = 8.6;

Vm = (4*pi*(21)^3)/24;

kP = 0.0072419;

Qstar = 1.1;

gamaP = 0.05;

alphaP = 212.95;

bP = 308;

nP = 2;

Tprod = 61.6;

gamaT = 0.01;

alphaT = 0.00075*144.87;

kS = 2/3;

kT = 0.003*3180;

nT = 2;

etammin = 0.38874;

etammax = 2.6828;

bm = 706;

etaemin = 0.41022;

etaemax = 0.69335;

be = 92.1;

% Time-delays

taum = 8.09;

taue = 5;

% Implementation of the dde23 solver to solve system of DDEs

tspan = [0 30];

lags = [taum taue (taue + taum)]; % constant time-delays of the system

sol = dde23(@ddefunc,lags,@Xhist,tspan);

t = linspace(0,30,1000);

P = deval(sol,t,1);

T = deval(sol,t,2);

%Plots

plot(t,P);

grid on

figure

plot(t,T);

grid on

% Definion of our DDE system

function dXdt = ddefunc(t,X,Z)

t

lag1 = Z(:,1); % lag1 = taum

lag2 = Z(:,2); % lag2 = taue

lag3 = Z(:,3); % lag3 = taue + taum

dXdt = [D0*Vm*kP*Qstar*X(4)*X(5)/betaP - gamaP*X(1) - alphaP*(X(1)^nP/(bP^nP + X(1)^nP)); % X(1) = P

Tprod - gamaT*X(2) - alphaT*(X(6) + kS*betaP*X(1))*(X(2)^nT/(kT^nT + X(2)^nT)); % X(2) = T

X(3)*((etammax - etammin)*(X(2)/(bm +X(2)) - lag1(2)/(bm + lag1(2)))); % X(3) = phi1

X(4)*((etammax - etammin)*(lag2(2)/(bm + lag2(2)) - lag3(2)/(bm + lag3(2)))); % X(4) = phi2

X(5)*((etaemax - etaemin)*(X(2)/(be + X(2)) - lag2(2)/(be + lag2(2)))); % X(5) = psi

Vm*kP*Qstar*(X(3) - X(4)*X(5)) + X(6)*(etaemin + (etaemax - etaemin)*X(2)/(be + X(2)))]; % X(6) = Me

end

% Definition of the history functions (how the system acts for t<0)

function h = Xhist(~)

Pstar = 31.071;

Tstar = 100;

P0 = Pstar;

T0 = Tstar + 100;

phi10 = exp(taum*(etammin + (etammax - etammin)*Tstar/(bm + Tstar))); % = 232.178224

psi0 = exp(taue*(etaemin + (etaemax - etaemin)*Tstar/(be + Tstar))); % = 16.249041

Me0 = Vm*kP*Qstar*phi10*(psi0-1)/(etaemin + (etaemax - etaemin)*Tstar/(be + Tstar)); % = 245,266.1919

h = [P0;

T0;

phi10;

phi10;

psi0;

Me0];

end

end

I should be getting periodic plots for P and T (I am trying to reproduce someone else's plots) but I'm not. Does anyone have any idea why that might be the case?

Answer 1

[Nipun]

Hi Sofia,

I understand that you are trying to reproduce periodic plots for (P) and (T) using a system of delay differential equations (DDEs) solved by MATLAB's dde23 function but are not observing the expected behavior. To address this, ensure the accuracy of your model parameters, initial conditions, and history function. Additionally, adjusting the numerical solver's settings may help capture the system's dynamics more accurately. Here's a concise approach to potentially resolve the issue:

1. Verify Model Parameters and Functions: Double-check that all parameters and the mathematical model accurately reflect the intended system.
2. Adjust Solver Settings: Use odeset to adjust solver tolerances, which might improve the accuracy of the solution, revealing the expected periodic behavior.
3. Ensure Sufficient Simulation Time and Sampling: Make sure your tspan and the number of points sampled (e.g., in linspace) are sufficient to capture the periodic dynamics.
4. Review Initial and History Conditions: Confirm that the history function Xhist accurately represents the system's state for (t < 0).

options = odeset('RelTol',1e-6, 'AbsTol',1e-8);
sol = dde23(@ddefunc, lags, @Xhist, tspan, options);

Incorporating these considerations, here's how you could adjust your MATLAB code:

function mainfunc()
    % [Your existing parameters here]

    % Time-delays
    taum = 8.09; 
    taue = 5;
    lags = [taum, taue, (taue + taum)]; % constant time-delays

    % Solver options for improved accuracy
    options = odeset('RelTol',1e-6, 'AbsTol',1e-8);

    % Implementation of the dde23 solver
    tspan = [0, 30];
    sol = dde23(@ddefunc, lags, @Xhist, tspan, options);

    % [Continue with your plotting code]
end

This approach, especially the solver settings adjustment, will help in capturing the expected periodic behavior of your system.

Hope this helps.

Regards,

Nipun