mackeyglass
This script generates a Mackey-Glass time series using the 4th order Runge-Kutta method. The code is a straighforward translation in Matlab of C source code provided by Roger Jang, which is available here
Contents
The theory
Mackey-Glass time series refers to the following, delayed differential equation:
It can be numerically solved using, for example, the 4th order Runge-Kutta method, at discrete, equally spaced time steps:
where the function mackeyglass_rk4 numerically solves the Mackey-Glass delayed differential equation using the 4-th order Runge Kutta. This is the RK4 method:
where mackeyglass_eq is the function which return the value of the Mackey-Glass delayed differential equation in (1) once its inputs and its parameters (a,b) are provided.
Here is an example:
Input parameters
a = 0.2; % value for a in eq (1) b = 0.1; % value for b in eq (1) tau = 17; % delay constant in eq (1) x0 = 1.2; % initial condition: x(t=0)=x0 deltat = 0.1; % time step size (which coincides with the integration step) sample_n = 12000; % total no. of samples, excluding the given initial condition interval = 1; % output is printed at every 'interval' time steps
Main algorithm
- x_t : x at instant t , i.e. x(t) (current value of x)
- x_t_minus_tau : x at instant (t-tau) , i.e. x(t-tau)
- x_t_plus_deltat : x at instant (t+deltat), i.e. x(t+deltat) (next value of x)
- X : the (sample_n+1)-dimensional vector containing x0 plus all other computed values of x
- T : the (sample_n+1)-dimensional vector containing time samples
- x_history : a circular vector storing all computed samples within x(t-tau) and x(t)
time = 0; index = 1; history_length = floor(tau/deltat); x_history = zeros(history_length, 1); % here we assume x(t)=0 for -tau <= t < 0 x_t = x0; X = zeros(sample_n+1, 1); % vector of all generated x samples T = zeros(sample_n+1, 1); % vector of time samples for i = 1:sample_n+1, X(i) = x_t; if (mod(i-1, interval) == 0), disp(sprintf('%4d %f', (i-1)/interval, x_t)); end if tau == 0, x_t_minus_tau = 0.0; else x_t_minus_tau = x_history(index); end x_t_plus_deltat = mackeyglass_rk4(x_t, x_t_minus_tau, deltat, a, b); if (tau ~= 0), x_history(index) = x_t_plus_deltat; index = mod(index, history_length)+1; end time = time + deltat; T(i) = time; x_t = x_t_plus_deltat; end figure plot(T, X); set(gca,'xlim',[0, T(end)]); xlabel('t'); ylabel('x(t)'); title(sprintf('A Mackey-Glass time serie (tau=%d)', tau));
0 1.200000 1 1.188060 2 1.176238 3 1.164535 4 1.152947 5 1.141475 6 1.130117 7 1.118873 8 1.107740 9 1.096717 10 1.085805 11 1.075001 12 1.064305 13 1.053715 14 1.043230 15 1.032850 16 1.022573 17 1.012398 18 1.002324 19 0.992351 20 0.982477 ...
