# include # include # include # include # include # include "midpoint.h" int main ( ); void predator_prey_dydt ( double t, double y[], double dy[] ); void predator_prey_test ( double tspan[2], double y0[2], int n ); void stiff_dydt ( double t, double y[], double dy[] ); double *stiff_exact ( int n, double t[] ); void stiff_test ( double tspan[2], double y0[1], int n ); void plot2 ( int n1, double t1[], double y1[], int n2, double t2[], double y2[], char *header, char *title ); double *r8vec_linspace_new ( int n, double a, double b ); void timestamp ( ); /******************************************************************************/ int main ( ) /******************************************************************************/ /* Purpose: midpoint_test() tests midpoint(). Licensing: This code is distributed under the MIT license. Modified: 09 November 2023 Author: John Burkardt */ { int n; double tspan[2]; double y0[2]; timestamp ( ); printf ( "\n" ); printf ( "midpoint_test():\n" ); printf ( " C version\n" ); printf ( " Test midpoint() on several ODE's.\n" ); tspan[0] = 0.0; tspan[1] = 5.0; y0[0] = 5000.0; y0[1] = 100.0; n = 200; predator_prey_test ( tspan, y0, n ); tspan[0] = 0.0; tspan[1] = 1.0; y0[0] = 0.0; n = 27; stiff_test ( tspan, y0, n ); /* Terminate. */ printf ( "\n" ); printf ( "midpoint_test():\n" ); printf ( " Normal end of execution.\n" ); printf ( "\n" ); timestamp ( ); return ( 0 ); } /******************************************************************************/ void predator_prey_test ( double tspan[2], double y0[2], int n ) /******************************************************************************/ /* Purpose: predator_prey_test() solves predator_prey_ode() using midpoint(). Discussion: The physical system under consideration is a pair of animal populations. The PREY reproduce rapidly for each animal alive at the beginning of the year, two more will be born by the end of the year. The prey do not have a natural death rate instead, they only die by being eaten by the predator. Every prey animal has 1 chance in 1000 of being eaten in a given year by a given predator. The PREDATORS only die of starvation, but this happens very quickly. If unfed, a predator will tend to starve in about 1/10 of a year. On the other hand, the predator reproduction rate is dependent on eating prey, and the chances of this depend on the number of available prey. The resulting differential equations can be written: PREY(0) = 5000 PRED(0) = 100 d PREY / dT = 2 * PREY(T) - 0.001 * PREY(T) * PRED(T) d PRED / dT = - 10 * PRED(T) + 0.002 * PREY(T) * PRED(T) Here, the initial values (5000,100) are a somewhat arbitrary starting point. The pair of ordinary differential equations that result have an interesting behavior. For certain choices of the interaction coefficients (such as those given here), the populations of predator and prey will tend to a periodic oscillation. The two populations will be out of phase the number of prey will rise, then after a delay, the predators will rise as the prey begins to fall, causing the predator population to crash again. There is a conserved quantity, which here would be: E(r,f) = 0.002 r + 0.001 f - 10 ln(r) - 2 ln(f) Licensing: This code is distributed under the MIT license. Modified: 09 November 2023 Author: John Burkardt Reference: George Lindfield, John Penny, Numerical Methods Using MATLAB, Second Edition, Prentice Hall, 1999, ISBN: 0-13-012641-1, LC: QA297.P45. Input: double tspan[2] = [ T0, TMAX ], the initial and final times. A reasonable value might be [ 0, 5 ]. double y0[2]: = [ PREY, PRED ], the initial number of prey and predators. A reasonable value might be [ 5000, 100 ]. int n: the number of time steps. */ { char command_filename[80]; FILE *command; char data_filename[80]; FILE *data; char header[] = "predator_prey"; int i; const int m = 2; double *t; double *yout; printf ( "\n" ); printf ( "predator_prey_test():\n" ); printf ( " A pair of ordinary differential equations for a population\n" ); printf ( " of predators and prey are solved using midpoint().\n" ); t = ( double * ) malloc ( ( n + 1 ) * sizeof ( double ) ); yout = ( double * ) malloc ( ( n + 1 ) * m * sizeof ( double ) ); midpoint ( predator_prey_dydt, tspan, y0, n, m, t, yout ); /* Create the data file. */ strcpy ( data_filename, header ); strcat ( data_filename, "_data.txt" ); data = fopen ( data_filename, "wt" ); for ( i = 0; i <= n; i++ ) { fprintf ( data, " %g %g %g\n", t[i], yout[0+i*m], yout[1+i*m] ); } fclose ( data ); printf ( "\n" ); printf ( " predator_prey_test: data stored in \"%s\".\n", data_filename ); /* Create the command file. */ strcpy ( command_filename, header ); strcat ( command_filename, "_commands.txt" ); command = fopen ( command_filename, "wt" ); fprintf ( command, "# %s\n", command_filename ); fprintf ( command, "#\n" ); fprintf ( command, "# Usage:\n" ); fprintf ( command, "# gnuplot < %s\n", command_filename ); fprintf ( command, "#\n" ); fprintf ( command, "set term png\n" ); fprintf ( command, "set output '%s.png'\n", header ); fprintf ( command, "set xlabel '<-- Predator -->'\n" ); fprintf ( command, "set ylabel '<-- Prey -->'\n" ); fprintf ( command, "set title 'predator prey: midpoint method'\n" ); fprintf ( command, "set grid\n" ); fprintf ( command, "set style data lines\n" ); fprintf ( command, "plot '%s' using 2:3 with lines lw 3\n", data_filename ); fprintf ( command, "quit\n" ); fclose ( command ); printf ( " predator_prey_test: plot commands stored in \"%s\".\n", command_filename ); free ( yout ); free ( t ); return; } /******************************************************************************/ void predator_prey_dydt ( double t, double y[], double dydt[] ) /******************************************************************************/ /* Purpose: predator_prey_dydt() evaluates the right hand side of predator_prey_ode(). Licensing: This code is distributed under the MIT license. Modified: 27 April 2020 Author: John Burkardt Reference: George Lindfield, John Penny, Numerical Methods Using MATLAB, Second Edition, Prentice Hall, 1999, ISBN: 0-13-012641-1, LC: QA297.P45. Input: double t: the current time. double y[2]: the current solution variables, rabbits and foxes. Output: double dydt[2]: the right hand side of the 2 ODE's. */ { dydt[0] = 2.0 * y[0] - 0.001 * y[0] * y[1]; dydt[1] = - 10.0 * y[1] + 0.002 * y[0] * y[1]; return; } /******************************************************************************/ void stiff_test ( double tspan[2], double y0[1], int n ) /******************************************************************************/ /* Purpose: stiff_test() uses midpoint() on stiff_ode(). Licensing: This code is distributed under the MIT license. Modified: 07 April 2021 Author: John Burkardt Input: double TSPAN(2): the first and last times. double Y0: the initial condition. int N: the number of steps to take. */ { const int m = 1; const int n2 = 101; double *t1; double *t2; double *y1; double *y2; printf ( "\n" ); printf ( "stiff_test():\n" ); printf ( " Solve stiff_ode() using midpoint().\n" ); t1 = ( double * ) malloc ( ( n + 1 ) * sizeof ( double ) ); y1 = ( double * ) malloc ( ( n + 1 ) * sizeof ( double ) ); midpoint ( stiff_dydt, tspan, y0, n, m, t1, y1 ); t2 = r8vec_linspace_new ( n2, tspan[0], tspan[1] ); y2 = stiff_exact ( n2, t2 ); plot2 ( n+1, t1, y1, n2, t2, y2, "stiff", "Stiff ODE: midpoint method" ); free ( t1 ); free ( t2 ); free ( y1 ); free ( y2 ); return; } /******************************************************************************/ void stiff_dydt ( double t, double y[], double dy[] ) /******************************************************************************/ /* stiff_deriv() evaluates the right hand side of stiff_ode(). Discussion: y' = 50 * ( cos(t) - y ) y(0) = 0 Licensing: This code is distributed under the MIT license. Modified: 28 April 2020 Author: John Burkardt Input: double T, Y[1]: the time and solution value. Output: double DYDT[1]: the derivative value. */ { dy[0] = 50.0 * ( cos ( t ) - y[0] ); return; } /******************************************************************************/ double *stiff_exact ( int n, double t[] ) /******************************************************************************/ /* Purpose: stiff_exact() evaluates the exact solution of stiff_ode(). Discussion: y' = 50 * ( cos(t) - y ) y(0) = 0 Licensing: This code is distributed under the MIT license. Modified: 28 April 2020 Author: John Burkardt Input: int N: the number of values. double T[N]: the evaluation times. Output: double STIFF_EXACT[N]: the exact solution values. */ { int i; double *y; y = ( double * ) malloc ( n * sizeof ( double ) ); for ( i = 0; i < n; i++ ) { y[i] = 50.0 * ( sin ( t[i] ) + 50.0 * cos ( t[i] ) - 50.0 * exp ( - 50.0 * t[i] ) ) / 2501.0; } return y; } /******************************************************************************/ void plot2 ( int n1, double t1[], double y1[], int n2, double t2[], double y2[], char *header, char *title ) /******************************************************************************/ /* Purpose: plot2() plots two curves together. Licensing: This code is distributed under the MIT license. Modified: 28 April 2020 Author: John Burkardt Input: int N1: the size of the first data set. double T1(N1), Y1(N1), the first dataset. int N2: the size of the second data set. double T2(N2), Y2(N2), the secod dataset. char *HEADER: an identifier for the data. char *TITLE: a title to appear in the plot. */ { char command_filename[80]; FILE *command; char data1_filename[80]; char data2_filename[80]; FILE *data; int i; /* Create the data file. */ strcpy ( data1_filename, header ); strcat ( data1_filename, "_data1.txt" ); data = fopen ( data1_filename, "wt" ); for ( i = 0; i < n1; i++ ) { fprintf ( data, " %g %g\n", t1[i], y1[i] ); } fclose ( data ); strcpy ( data2_filename, header ); strcat ( data2_filename, "_data2.txt" ); data = fopen ( data2_filename, "wt" ); for ( i = 0; i < n2; i++ ) { fprintf ( data, " %g %g\n", t2[i], y2[i] ); } fclose ( data ); printf ( "\n" ); printf ( " stiff_ode_test: data stored in '%s' and '%s'.\n", data1_filename, data2_filename ); /* Create the command file. */ strcpy ( command_filename, header ); strcat ( command_filename, "_commands.txt" ); command = fopen ( command_filename, "wt" ); fprintf ( command, "# %s\n", command_filename ); fprintf ( command, "#\n" ); fprintf ( command, "# Usage:\n" ); fprintf ( command, "# gnuplot < %s\n", command_filename ); fprintf ( command, "#\n" ); fprintf ( command, "set term png\n" ); fprintf ( command, "set output '%s.png'\n", header ); fprintf ( command, "set xlabel '<-- T -->'\n" ); fprintf ( command, "set ylabel '<-- Y(T) -->'\n" ); fprintf ( command, "set title '%s'\n", title ); fprintf ( command, "set grid\n" ); fprintf ( command, "set style data lines\n" ); fprintf ( command, "plot '%s' using 1:2 with lines lw 3 lt rgb 'red',\\\n", data1_filename ); fprintf ( command, " '%s' using 1:2 with lines lw 3 lt rgb 'blue'\n", data2_filename ); fprintf ( command, "quit\n" ); fclose ( command ); printf ( " stiff_ode_test: plot commands stored in \"%s\".\n", command_filename ); return; } /******************************************************************************/ double *r8vec_linspace_new ( int n, double a, double b ) /******************************************************************************/ /* Purpose: r8vec_linspace_new() creates a vector of linearly spaced values. Discussion: An R8VEC is a vector of R8's. 4 points evenly spaced between 0 and 12 will yield 0, 4, 8, 12. In other words, the interval is divided into N-1 even subintervals, and the endpoints of intervals are used as the points. Licensing: This code is distributed under the MIT license. Modified: 29 March 2011 Author: John Burkardt Parameters: Input, int N, the number of entries in the vector. Input, double A, B, the first and last entries. Output, double R8VEC_LINSPACE_NEW[N], a vector of linearly spaced data. */ { int i; double *x; x = ( double * ) malloc ( n * sizeof ( double ) ); if ( n == 1 ) { x[0] = ( a + b ) / 2.0; } else { for ( i = 0; i < n; i++ ) { x[i] = ( ( double ) ( n - 1 - i ) * a + ( double ) ( i ) * b ) / ( double ) ( n - 1 ); } } return x; } /******************************************************************************/ void timestamp ( ) /******************************************************************************/ /* Purpose: timestamp() prints the current YMDHMS date as a time stamp. Example: 17 June 2014 09:45:54 AM Licensing: This code is distributed under the MIT license. Modified: 17 June 2014 Author: John Burkardt */ { # define TIME_SIZE 40 static char time_buffer[TIME_SIZE]; const struct tm *tm; time_t now; now = time ( NULL ); tm = localtime ( &now ); strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm ); printf ( "%s\n", time_buffer ); return; # undef TIME_SIZE }