Euller method to solve mackey glass system

In my understanding , dt needs to be very small so that the Euller method works.

I am solving the mackey glass equation

dx(t)/dt = 0.2x(t-tau)/(1+x(t-tau)^10) - 0.1x(t)

If I use dt = 1, the graph is very "mackey-glassy", but for any value different than 1, the system converges to a step function.

What am I doing wrong?

Thanks,

x = make_mg_euller(0.3,17);

Warning: Function Warning: Name is nonexistent or not a directory: /users/mss.system.452ATJ/./mackeyglass 
> In path (line 109)
In addpath (line 86)
In addpath (line 47)
In LiveEditorEvaluationHelperEeditorId>make_mg_euller (line 8)
In LiveEditorEvaluationHelperEeditorId (line 1)
In connector.internal.fevalMatlab
In connector.internal.fevalJSON 

figure; 
plot(x)

function [x] = make_mg_euller(x0,tau)
addpath ./mackeyglass/


% Euler's Method
% Initial conditions and setup
dt = 0.01;  % step size
n = 400000; % simulation steps
x = double(zeros(n,1));  % allocate
x(1) = x0; % initial condition

% The loop to solve the DE
    for i = 1:tau-1
    x(i+1) = x(i) + double (dt * ( -  0.1 * x(i) ) );
    end
    for i = tau:n-1
    x(i+1) = x(i) + double(dt*  ( ( 0.2 * x(i-tau+1) ) / ( 1 + x(i-tau+1)^10 )  - ( 0.1 * x(i) ) ))  ;
    end
end

[Torsten]

You try to solve a delay differential equation. Use dde23.

If you want to solve it with your own code, note that tau has to be a multiple of dt, and you must access x(i-j) to get x(t-tau) if tau = j*dt. For t <= tau, x must explicitly be prescribed as a "history" function or by an ordinary differential equation without delay terms and an initial condition at t=0.

[L]

Thanks for the tip @Torsten.

Answer 1

[LeoAiE]

The main issue with your implementation is that you're using the parameter tau as an integer representing the time delay in terms of steps instead of a continuous time value. With this, when you use a different dt, the delay in the equation isn't adjusted properly.

To fix the issue, you need to calculate the number of steps that correspond to the time delay tau based on the dt value. Here's an updated version of your code with the necessary changes:

x = make_mg_euller(0.3,17);
figure; 
plot(x)

function [x] = make_mg_euller(x0, tau)
    % Euler's Method
    % Initial conditions and setup
    dt = 0.01; % step size
    n = 400000; % simulation steps
    x = double(zeros(n, 1)); % allocate
    x(1) = x0; % initial condition

    % Calculate the number of steps corresponding to the time delay tau
    tau_steps = round(tau / dt);

    % The loop to solve the DE
    for i = 1:tau_steps - 1
        x(i + 1) = x(i) + double(dt * (-0.1 * x(i)));
    end

    for i = tau_steps:n - 1
        x(i + 1) = x(i) + double(dt * ((0.2 * x(i - tau_steps + 1)) / (1 + x(i - tau_steps + 1)^10) - (0.1 * x(i))));
    end
end

[Fares]

Hi LeoAiE, thanks for your answer. I am applying your code to my case but I am wondering why we need to ignore the delay terms in the first loop. Is it just enough to include them without the delay parameter? Many thanks!