/* NAME * mg - plot the phase space of the Mackey-Glass system * MISCELLANY * The plot region is determined by the points that are initially * skipped. If this number is too small (i.e., it is not very * representative of the range of the plotted values), then you * may need to increase the number specified by the -skip option. * Alternatively, you can adjust the value given to -factor, which * simply fractionally increases the border of the plot. * * The program uses a second-order Euler's method to perform the * numerical integration, which is sufficient for simple tasks such * as this. * BUGS * No sanity checks are performed to make sure that any of the * options make sense. * AUTHOR * Copyright (c) 1997, Gary William Flake. * * Permission granted for any use according to the standard GNU * ``copyleft'' agreement provided that the author's comments are * neither modified nor removed. No warranty is given or implied. */ #include #include #include #include "misc.h" int width = 480, height = 480, skip = 2000, points = 10000; int delta = 6, data = 0, invert = 0, mag = 1, tau = 17; double a = 0.2, b = 0.1, dt = 0.1, factor = 0.2, xx0 = 1.23456789; char *term = NULL; char help_string[] = "\ The phase space of the Mackey-Glass system, which is described by the \ delay differential equation \ \ dx(t)/dt = A * x(t-Tau) / (1 + x(t-Tau)^10) - B * x(t), \ \ is plotted according to the specified parameters. The x-coordinate of \ the plot is determined by x(t) while the y-coordinate is determined by \ x(t-delta).\ "; OPTION options[] = { { "-width", OPT_INT, &width, "Width of the plot in pixels." }, { "-height", OPT_INT, &height, "Height of the plot in pixels." }, { "-skip", OPT_INT, &skip, "Number of initial points to skip." }, { "-points", OPT_INT, &points, "Number of points to plot." }, { "-delta", OPT_INT, &delta, "Time steps to delay for." }, { "-tau", OPT_INT, &tau, "Value of the Tau parameter." }, { "-A", OPT_DOUBLE, &a, "Value of the A parameter." }, { "-B", OPT_DOUBLE, &b, "Value of the B parameter." }, { "-dt", OPT_DOUBLE, &dt, "Time step size." }, { "-x0", OPT_DOUBLE, &xx0, "Initial X value." }, { "-factor", OPT_DOUBLE, &factor, "Auto-scale expansion factor." }, { "-data", OPT_SWITCH, &data, "Don't plot, but print points." }, { "-inv", OPT_SWITCH, &invert, "Invert all colors?" }, { "-mag", OPT_INT, &mag, "Magnification factor." }, { "-term", OPT_STRING, &term, "How to plot points." }, { NULL, OPT_NULL, NULL, NULL } }; /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ /* Perform a second-order Euler step according to the MG system equation. */ double euler(double xt, double xtau, double a, double b, double dt) { double x1, x2; x1 = (a * xtau) / (1 + pow(xtau, 10.0)) - b * xt; x2 = (a * xtau) / (1 + pow(xtau, 10.0)) - b * (dt * x1 + xt); return(xt + 0.5 * dt * (x1 + x2)); } /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ int main(int argc, char **argv) { extern int plot_mag; extern int plot_inverse; int i, hsz, h, ttau, tdelta; /* Variables used to compute auto-scaling. */ double xmin, xmax, ymin, ymax; double x, y, xx = 0, yy = 0, xt, xtau, temp, *hold; get_options(argc, argv, options, help_string); if(!data) { plot_mag = mag; plot_inverse = invert; plot_init(width, height, 2, term); plot_set_all(0); } /* Recalculate Tau and delta to reflect time step size. */ ttau = 1 / dt * tau + 0.5; tdelta = 1 / dt * delta + 0.5; /* Initialize the buffer space. */ hsz = MAX(tdelta + 1, ttau + 1); hold = xmalloc(hsz * sizeof(double)); /* Initialize the system. */ h = 0; for(i = 0; i < hsz; i++) hold[i] = xx0; /* Set insane minimum and maximum guesses. */ xmin = ymin = 10e10; xmax = ymax = -10e10; /* For all points (plus the skip and delay values. */ for(i = 0; i < points + skip + tdelta + 1; i++) { /* Compute the time evolution of the system. */ xtau = hold[(h + hsz - (ttau + 1)) % hsz]; xt = hold[(h + hsz - 1) % hsz]; x = euler(xt, xtau, a, b, dt); /* Save the state so that we can remember delayed values. */ hold[h] = x; h++; h %= hsz; /* Get the delayed value. */ y = hold[(h + hsz - tdelta) % hsz]; if(data) { if(i >= skip + delta + 1) printf("%f\t%f\n", x, y); } else { /* If we are still skipping points, adjust the best guesses for * the minimum and maximums. */ if(i < skip + delta + 1) { xmin = (x < xmin) ? x : xmin; xmax = (x > xmax) ? x : xmax; ymin = (y < ymin) ? y : ymin; ymax = (y > ymax) ? y : ymax; } /* If this is the last point to be skipped, reset the plotting * range based on the minimum and maximums. */ if(i == skip + delta + 1) { temp = (xmax - xmin) * factor; xmin -= temp; xmax += temp; temp = (ymax - ymin) * factor; ymin -= temp; ymax += temp; plot_set_range(xmin, xmax, ymin, ymax); } /* Plot a line from the last point to the current point. */ if(i >= skip + delta + 1) plot_line(xx, yy, x, y, 1); /* Save the last point. */ xx = x; yy = y; } } if(!data) plot_finish(); exit(0); } /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */