Hi Muhammad,
Transitioning from solving a single delay differential equation (DDE) to handling a system of DDEs with MATLAB's "dde23" function involves a few key adjustments to ensure that your code correctly manages vector inputs and outputs, as DDE systems involve multiple state variables. Here are some points and suggestions for modifying your existing code to handle a system of DDEs:
- Ensure g Function Handles Vector Inputs: The function g(t,y,Z,par) that you pass to "dde23" must be capable of handling vector inputs for 'y' and 'Z'.
- Initial History Function phi as a Vector Function: The initial history function phi should return a vector of values corresponding to the initial states of all equations in your system for times 't' <= 0.
- Adjustment for Multiple State Variables and Delays: While dealing with multiple state variables and delays, ensure that your code correctly handles indexing into these vectors/matrices to compute the current and delayed states for each variable.
- Vectorizing the Computation for Efficiency: Try to vectorize the computations as much as possible instead of using nested loops, which can significantly slow down MATLAB code, especially for large systems or long time spans.
The optimized MATLAB code with demo equations for illustration is as follows:
% Define parameters (if any)
par = []; % Placeholder for parameters
% Define delays
tau = [1]; % Example delay
% Time span
tspan = [0, 10];
% Initial history function (for t <= 0)
% It should return a vector for a system of DDEs
phi = @(t) [1; 0.5]; % Example initial values for two state variables
% Define the DDE system
% y and Z are vectors. Z(:,i) is the state at t - tau(i)
g = @(t, y, Z, par) [
-y(1) + Z(1,1); % Example equation for the first state variable
-0.5 * y(2) + 0.5 * Z(1,1) % Example equation for the second state variable
];
% Solve the system of DDEs
sol = dde23(@(t, y, Z) g(t, y, Z, par), tau, phi, tspan);
% User input parameters (for example purposes)
dist = 1; % Distribution type
t_r = 0.5; % Time ratio
N = 100; % Number of points
n_a = 0.1; % Noise amplitude
% Select data points based on user input
if dist == 1
t_ae = linspace(tspan(1), tspan(2), N);
else
t_ae = tspan(1) + (tspan(2) - tspan(1)) * rand(1, N);
end
% Find current state terms for all variables
x_t = deval(sol, t_ae); % x_t is now a matrix with rows corresponding to variables
% Add noise to the current state terms
x_tn = x_t + n_a * randn(size(x_t));
% Compute delayed state terms for all variables and delays
x_d = zeros(size(x_t, 1), length(tau), length(t_ae)); % Adjusted for multiple delays
for k = 1:length(tau)
t_d = t_ae - tau(k);
for i = 1:length(t_d)
if t_d(i) <= 0
x_d(:, k, i) = phi(t_d(i)); % Adjusted to handle vector output from phi
else
for j = 1:size(x_t, 1) % For each state variable
x_d(j, k, i) = interp1(t_ae, x_tn(j, :), t_d(i), 'linear', 'extrap');
end
end
end
end
% At this point, x_t and x_d contain the current and delayed states of the system,
% which you can use for further analysis or visualization.
This code considers a hypothetical system of two DDEs with one delay. You may need adjustments based on the specifics of your system of DDEs.
You may refer to the following documentation links to have a better understanding on working with systems of delay differential equations (DDEs) in MATLAB:
