# include # include # include # include # include # include "stiff_exact.h" int main ( ); void stiff_value ( ); void stiff_euler_test ( int n ); void stiff_euler ( ); void plot2 ( int n1, double t1[], double y1[], int n2, double t2[], double y2[], char *header, char *title ); void plot4 ( int n1, double t1[], double y1[], int n2, double t2[], double y2[], int n3, double t3[], double y3[], int n4, double t4[], double y4[], char *header, char *title ); double *r8vec_linspace_new ( int n, double a, double b ); void timestamp ( ); /******************************************************************************/ int main ( ) /******************************************************************************/ /* Purpose: stiff_exact_test() tests stiff_exact(). Licensing: This code is distributed under the MIT license. Modified: 23 June 2025 Author: John Burkardt */ { int n; timestamp ( ); printf ( "\n" ); printf ( "stiff_exact_test():\n" ); printf ( " C version\n" ); printf ( " Test stiff_exact(), which evaluates exact solutions\n" ); printf ( " of a stiff ODE.\n" ); stiff_value ( ); n = 27; stiff_euler_test ( n ); /* Terminate. */ printf ( "\n" ); printf ( "stiff_exact_test():\n" ); printf ( " Normal end of execution.\n" ); timestamp ( ); return 0; } /******************************************************************************/ void stiff_value ( ) /******************************************************************************/ /* Purpose: stiff_value() evaluates the exact solution for various initial conditions. Discussion: y' = lambda * ( cos(t-t0) - y ) y(t0) = y0 Licensing: This code is distributed under the MIT license. Modified: 23 June 2025 Author: John Burkardt */ { double lambda; int nt; double t0; double *t1; double *t2; double *t3; double *t4; double tstop; double y0; double *y1; double *y2; double *y3; double *y4; nt = 101; /* Set initial parameter values. */ lambda = 50.0; t0 = 0.0; y0 = 0.0; tstop = 1.0; stiff_parameters ( &lambda, &t0, &y0, &tstop, NULL, NULL, NULL, NULL ); /* Report initial parameter values. */ printf ( "\n" ); printf ( " Initial parameters:\n" ); printf ( " lambda = %g\n", lambda ); printf ( " t0 = %g\n", t0 ); printf ( " y0 = %g\n", y0 ); printf ( " tstop = %g\n", tstop ); /* Default with t0 = 0, y0 = 0. */ t0 = 0.0; y0 = 0.0; stiff_parameters ( &lambda, &t0, &y0, &tstop, NULL, NULL, NULL, NULL ); t1 = r8vec_linspace_new ( nt, t0, tstop ); y1 = stiff_exact ( nt, t1 ); /* Try t0=0, y0=0.9. */ t0 = 0.0; y0 = 0.9; stiff_parameters ( &lambda, &t0, &y0, &tstop, NULL, NULL, NULL, NULL ); t2 = r8vec_linspace_new ( nt, t0, tstop ); y2 = stiff_exact ( nt, t2 ); /* Try t0=0, y0=1.5. */ t0 = 0.0; y0 = 1.5; stiff_parameters ( &lambda, &t0, &y0, &tstop, NULL, NULL, NULL, NULL ); t3 = r8vec_linspace_new ( nt, t0, tstop ); y3 = stiff_exact ( nt, t3 ); /* t0 = 0.25, y0 = 0.5 */ t0 = 0.25; y0 = 0.5; stiff_parameters ( &lambda, &t0, &y0, &tstop, NULL, NULL, NULL, NULL ); t4 = r8vec_linspace_new ( nt, t0, tstop ); y4 = stiff_exact ( nt, t4 ); /* Plot them. */ plot4 ( nt, t1, y1, nt, t2, y2, nt, t3, y3, nt, t4, y4, "stiff_value", "Stiff ODE exact solutions" ); /* Free memory. */ free ( t1 ); free ( t2 ); free ( t3 ); free ( t4 ); free ( y1 ); free ( y2 ); free ( y3 ); free ( y4 ); return; } /******************************************************************************/ void stiff_euler_test ( int n ) /******************************************************************************/ /* Purpose: stiff_euler_test() tests stiff_euler(). Licensing: This code is distributed under the MIT license. Modified: 23 June 2025 Author: John Burkardt Input: int N: the number of steps to take. */ { double lambda; const int n2 = 101; double t0; double *t1; double *t2; double tstop; double y0; double *y1; double *y2; printf ( "\n" ); printf ( "stiff_euler_test():\n" ); printf ( " Solve stiff ODE using the Euler method.\n" ); /* Set initial parameter values. */ lambda = 50; t0 = 0.0; y0 = 0.0; tstop = 1.0; stiff_parameters ( &lambda, &t0, &y0, &tstop, NULL, NULL, NULL, NULL ); printf ( "\n" ); printf ( " Parameters:\n" ); printf ( " lambda = %g\n", lambda ); printf ( " t0 = %g\n", t0 ); printf ( " y0 = %g\n", y0 ); printf ( " tstop = %g\n", tstop ); stiff_parameters ( NULL, NULL, NULL, NULL, NULL, &t0, NULL, &tstop ); t1 = ( double * ) malloc ( ( n + 1 ) * sizeof ( double ) ); y1 = ( double * ) malloc ( ( n + 1 ) * sizeof ( double ) ); stiff_euler ( n, t1, y1 ); t2 = r8vec_linspace_new ( n2, t0, tstop ); y2 = stiff_exact ( n2, t2 ); plot2 ( n+1, t1, y1, n2, t2, y2, "stiff_euler", "Stiff ODE solved using euler method" ); free ( t1 ); free ( t2 ); free ( y1 ); free ( y2 ); return; } /******************************************************************************/ void stiff_euler ( int n, double t[], double y[] ) /******************************************************************************/ /* Purpose: stiff_euler() uses the Euler method on the stiff equation. Discussion: y' = lambda * ( cos(t-t0) - y ) y(t0) = y0 Licensing: This code is distributed under the MIT license. Modified: 20 June 2025 Author: John Burkardt Input: int N: the number of steps to take. Output: double T[N+1], Y[N+1]: the times and estimated solutions. */ { double dt; int i; double lambda; double t0; double tstop; double y0; stiff_parameters ( NULL, NULL, NULL, NULL, &lambda, &t0, &y0, &tstop ); dt = ( tstop - t0 ) / ( double ) ( n ); t[0] = t0; y[0] = y0; for ( i = 0; i < n; i++ ) { t[i+1] = t[i] + dt; y[i+1] = y[i] + dt * lambda * ( cos ( t[i] - t0 ) - y[i] ); } return; } /******************************************************************************/ 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; } /******************************************************************************/ void plot4 ( int n1, double t1[], double y1[], int n2, double t2[], double y2[], int n3, double t3[], double y3[], int n4, double t4[], double y4[], char *header, char *title ) /******************************************************************************/ /* Purpose: plot4() plots four curves together. Licensing: This code is distributed under the MIT license. Modified: 23 June 2025 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 second dataset. int N3: the size of the data set 3. double T3(N3), Y3(N3), dataset 3. int N3: the size of the data set 4. double T4(N4), Y4(N4), dataset 4. 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]; char data3_filename[80]; char data4_filename[80]; FILE *data; int i; /* Create the data files. */ 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 ); strcpy ( data3_filename, header ); strcat ( data3_filename, "_data3.txt" ); data = fopen ( data3_filename, "wt" ); for ( i = 0; i < n3; i++ ) { fprintf ( data, " %g %g\n", t3[i], y3[i] ); } fclose ( data ); strcpy ( data4_filename, header ); strcat ( data4_filename, "_data4.txt" ); data = fopen ( data4_filename, "wt" ); for ( i = 0; i < n4; i++ ) { fprintf ( data, " %g %g\n", t4[i], y4[i] ); } fclose ( data ); printf ( "\n" ); printf ( " stiff_ode_test: data stored in '%s', '%s', '%s', and '%s'.\n", data1_filename, data2_filename, data3_filename, data3_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, " '%s' using 1:2 with lines lw 3 lt rgb 'green',\\\n", data3_filename ); fprintf ( command, " '%s' using 1:2 with lines lw 3 lt rgb 'magenta'\n", data4_filename ); fprintf ( command, "quit\n" ); fclose ( command ); printf ( " stiff_exact_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 Input: int N, the number of entries in the vector. 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: 01 May 2021 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 }