# include # include # include # include # include # include # include "fsolve.h" # include "midpoint_adaptive.h" /******************************************************************************/ void mad_compute ( void f ( double t, double y[], double dydt[] ), double t0, double tmax, int m, double y0[], double tau0, double reltol, double abstol, int nmax, int *nstep, int *n_rejected, int *n_fsolve, double t[], double y[] ) /******************************************************************************/ /* Purpose: mad_compute() estimates the solution of an ordinary differential equation. Discussion: mad_compute() uses an adaptive time stepping algorithm for the implicit midpoint method, applied to a system of ordinary differential equations. The implicit midpoint method has local truncation error O(h^3). The time step is adaptively controlled with respect to relative and absolute tolerances applied to the estimated local truncation error (LTE). The adaptivity with respect to relative and absolute error tolerances is described on page 168 in the Hairer reference. Licensing: This code is distributed under the MIT license. Modified: 15 July 2024 Author: John Burkardt, Catalin Trenchea Reference: Catalin Trenchea, John Burkardt, Refactorization of the midpoint rule, Applied Mathematics Letters, Volume 107, September 2020. Ernst Hairer, Syvert Norsett, Gerhard Wanner, Solving ordinary differential equations, I. Nonstiff problems, Springer Series in Computational Mathematics, Number 8. Springer-Verlag, Berlin, 1987. Input: void f ( double t, double y[], double dydt[] ): evaluates the right hand side of the differential equation. double t0, tmax: the interval of integration. int m: the dimension of a single solution vector. double y9(m): the vector of initial conditions at the starting time. double tau0: the initial timestep. double reltol, abstol: the tolerances for the local truncation error. int nmax: the maximum number of time steps allowed. Output: int *nstep: the number of time steps taken. int *n_rejected: the number of rejected time steps. int *n_fsolve: the number of failed calls to fsolve(). double t(0:nmax), y(0:nmax,1:m): contains in entries 0:nstep the the sequence of times and solution estimates actually computed. Local: int count: the number of steps, whether accepted or rejected. It should be the case that n <= count. double kappa: a factor for the stepsize update. 0 < kappa <= 1. The value kappa = 0.85 is suggested. double theta: the theta-method parameter. theta = 0.0 for the backward Euler method. theta = 0.5 for the midpoint method. theta = 1.0 for the forward Euler method. */ { int count; double factor; bool fsolve_fail; int i; double kappa; int n; double *tau; double theta; double tn; double tnmax; double tnmin; double *yn; tau = ( double * ) malloc ( ( nmax + 1 ) * sizeof ( double ) ); yn = ( double * ) malloc ( m * sizeof ( double ) ); kappa = 0.85; theta = 0.5; *n_rejected = 0; *n_fsolve = 0; /* On starting this loop with a given value of n, we have t(0:n-1), y(0:n-1,1:m) and tau(0:n-1). We want to take step n, and compute t(n), y(n,1:m), tau(n). */ for ( n = 0; n <= nmax; n++ ) { /* n = 0: set n data */ if ( n == 0 ) { t[n] = t0; for ( i = 0; i < m; i++ ) { y[i+n*m] = y0[i]; } tau[n] = tau0; count = 0; } /* n = 1, 2: y1, y2 Steps 1 and 2 use fixed tau. */ else if ( n <= 2 ) { mad_step ( m, t[n-1], y+(n-1)*m, tau[n-1], f, theta, reltol, &tn, yn, &fsolve_fail ); if ( fsolve_fail ) { printf ( "\n" ); printf ( "mad_compute(): Fatal error!\n" ); printf ( " Call to fsolve() failed on first or second step.\n" ); exit ( 1 ); } t[n] = tn; for ( i = 0; i < m; i++ ) { y[i+n*m] = yn[i]; } tau[n] = tau0; count = count + 1; } /* n = 3, ..., nmax: use data from step n-1 to compute data for step n. Tau gets updated. Trying a step to t(n) = t(n-1) + theta * tau(n-1). If step fails, may need to reduce theta and try again. */ else { while ( true ) { mad_step ( m, t[n-1], y+(n-1)*m, tau[n-1], f, theta, reltol, &tn, yn, &fsolve_fail ); if ( fsolve_fail ) { *n_fsolve = *n_fsolve + 1; continue; } t[n] = tn; for ( i = 0; i < m; i++ ) { y[i+n*m] = yn[i]; } count = count + 1; /* Estimate local truncation error. */ mad_lte ( m, y+(n-3)*m, tau+(n-3), reltol, abstol, &tnmin, &tnmax ); /* If the LTE test was met, we accept T and Y, and set the next time step. */ if ( tnmin <= 1.0 ) { factor = kappa * pow ( 1.0 / tnmax, 1.0/3.0 ); factor = fmin ( factor, 1.5 ); factor = fmax ( factor, 0.02 ); tau[n] = factor * tau[n-1]; break; } /* If the step was not accepted, we try again. */ *n_rejected = *n_rejected + 1; } } /* Exit if we reached or surpassed final time b. */ if ( tmax <= t[n] ) { break; } } *nstep = n; /* Free memory. */ free ( tau ); free ( yn ); return; } /******************************************************************************/ void mad_lte ( int m, double y[], double tau[], double reltol, double abstol, double *tnmin, double *tnmax ) /******************************************************************************/ /* Purpose: mad_lte() estimates local truncation error for midpoint adaptive method. Discussion: This function is called by mad_compute() to estimate the local truncation error incurred during a single step of the implicit midpoint method. We assume the step has been taken over an interval which we will represent as going from (t1,y1) to (t2,y2). We may regard y2 as the estimated solution of the local ODE: dy/dt = f(t,y), y(t1) = y1. We suppose that a function y*() is the exact solution of this system and we define the local truncation error for this step as lte = y*(t2) - y2 We wish to accept this step as long as lte is small, that is: lte = y*(t2) - y2 < abstol + reltol * max ( y1, y2 ) Of course, to carry out this check, we need to estimate lte or y*(t2). This is done using the Milne Device, described in the Milne reference. Essentially, lte is estimated using a weighted difference of solution estimates y2 and an estimate formed by an Adams-Bashforth AB2-like method. See page 5 of the Burkardt/Trenchea reference for details. Licensing: This code is distributed under the MIT license. Modified: 29 June 2024 Author: John Burkardt, Catalin Trenchea Reference: John Burkardt, Catalin Trenchea, Refactorization of the midpoint rule, Applied Mathematics Letters, Volume 107, 106438, 2020. William Milne, Numerical Integration of Ordinary Differential Equations, American Mathematical Monthly, Volume 33, number 9, pages 455–460, 1926. Input: integer m: the spatial dimension of the solution. real y[4*m]: the last four solution vectors. real tau[3]: the last three stepsizes. real reltol, abstol: local truncation error tolerances. Output: real tnmin, tnmax: the minimum and maximum entries of the local truncation error estimator. */ { double c1; double c2; double c3; int i; double r; double tn; double *uab2; /* Evaluate the AB2-like solution. */ c1 = ( tau[2] + tau[1] ) * ( tau[2] + tau[1] + tau[0] ) / ( tau[1] * ( tau[1] + tau[0] ) ); c2 = - tau[2] * ( tau[2] + tau[1] + tau[0] ) / ( tau[1] * tau[0] ); c3 = tau[2] * ( tau[2] + tau[1] ) / ( tau[0] * ( tau[1] + tau[0] ) ); uab2 = ( double * ) malloc ( m * sizeof ( double ) ); for ( i = 0; i < m; i++ ) { uab2[i] = c1 * y[i+2*m] + c2 * y[i+1*m] + c3 * y[i+0*m]; } /* Evaluate R, the variable error coefficient in the LTE. */ r = 1.0 / 24.0 + 1.0 / 8.0 * ( 1.0 + tau[1] / tau[2] ) * ( 1.0 + 2.0 * tau[1] / tau[2] + tau[0] / tau[2] ); /* Use AB2 estimator to approximate the LTE vector. */ *tnmax = 0; *tnmin = HUGE_VAL; for ( i = 0; i < m; i++ ) { tn = ( y[i+3*m] - uab2[i] ) * 24.0 * r / ( 24.0 * r - 1.0 ) / ( abstol + fabs ( y[i+3*m] ) * reltol ); tn = fabs ( tn ); if ( *tnmax < tn ) { *tnmax = tn; } if ( tn < *tnmin ) { *tnmin = tn; } } free ( uab2 ); return; } /******************************************************************************/ void mad_phase_plot ( int n, int m, double t[], double y[], char *label ) /******************************************************************************/ /* Purpose: mad_phase_plot() writes data and command files for a phase plot. Discussion: We assume a simple phase plane plot of y(0) versus y(1). Licensing: This code is distributed under the MIT license. Modified: 29 June 2024 Author: John Burkardt Input: int n: the number of steps. int m: the size of a solution vector. double t(0:n): the sequence of times. double y[n*m]: N solutions, each of length M. char *label: a title. */ { char command_filename[255]; FILE *command; char data_filename[255]; FILE *data; int i; int j; char png_filename[255]; /* Create the data file. */ strcpy ( data_filename, label ); strcat ( data_filename, "_phase_data.txt" ); data = fopen ( data_filename, "wt" ); for ( j = 0; j < n; j++ ) { fprintf ( data, " %g", t[j] ); for ( i = 0; i < m; i++ ) { fprintf ( data, " %g", y[i+j*m] ); } fprintf ( data, "\n" ); } fclose ( data ); printf ( "\n" ); printf ( " Phase graphics data stored in \"%s\".\n", data_filename ); /* Create the command file. */ strcpy ( command_filename, label ); strcat ( command_filename, "_phase_commands.txt" ); strcpy ( png_filename, label ); strcat ( png_filename, "_phase.png" ); 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'\n", png_filename ); fprintf ( command, "set xlabel '<-- y[0][t] -->'\n" ); fprintf ( command, "set ylabel '<-- y[1][t] -->'\n" ); fprintf ( command, "set title 'MAD phase plane'\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 ( " Phase plot commands stored in \"%s\".\n", command_filename ); return; } /******************************************************************************/ void mad_solution_plot ( int n, int m, double t[], double y[], char *label ) /******************************************************************************/ /* Purpose: mad_solution_plot() writes data and command files for a solution plot. Licensing: This code is distributed under the MIT license. Modified: 29 June 2024 Author: John Burkardt Input: int n: the number of steps. int m: the size of a solution vector. double t(0:n): the sequence of time values. double y(0:n * 1:m): the sequence of solutions. char *label: a title. */ { char command_filename[255]; FILE *command; char data_filename[255]; FILE *data; int i; int j; char png_filename[255]; /* Create the data file. */ strcpy ( data_filename, label ); strcat ( data_filename, "_solution_data.txt" ); data = fopen ( data_filename, "wt" ); for ( j = 0; j < n; j++ ) { fprintf ( data, " %g", t[j] ); for ( i = 0; i < m; i++ ) { fprintf ( data, " %g", y[i+j*m] ); } fprintf ( data, "\n" ); } fclose ( data ); printf ( "\n" ); printf ( " Solution data stored in \"%s\".\n", data_filename ); /* Create the command file. */ strcpy ( command_filename, label ); strcat ( command_filename, "_solution_commands.txt" ); strcpy ( png_filename, label ); strcat ( png_filename, "_solution.png" ); 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'\n", png_filename ); fprintf ( command, "set xlabel '<-- t -->'\n" ); fprintf ( command, "set ylabel '<-- y[t] -->'\n" ); fprintf ( command, "set title 'MAD solution Y(T)'\n" ); fprintf ( command, "set grid\n" ); fprintf ( command, "set style data lines\n" ); fprintf ( command, "plot \\\n" ); for ( i = 2; i <= m; i++ ) { fprintf ( command, "'%s' using 1:%d with lines lw 3,\\\n", data_filename, i ); } fprintf ( command, "'%s' using 1:%d with lines lw 3\n", data_filename, m + 1 ); fprintf ( command, "quit\n" ); fclose ( command ); printf ( " Solution plot commands stored in \"%s\".\n", command_filename ); return; } /******************************************************************************/ void mad_stats ( int n, int n_rejected, int n_fsolve, double t[] ) /******************************************************************************/ /* Purpose: mad_stats() reports time and timestep statistics for mad_compute(). Discussion: This function prints some information about the time nodes and time steps. Licensing: This code is distributed under the MIT license. Modified: 15 July 2024 Author: John Burkardt, Catalin Trenchea Reference: Catalin Trenchea, John Burkardt, Refactorization of the midpoint rule, Applied Mathematics Letters, Volume 107, September 2020. Input: int n: the number of time steps int n_rejected: the number of rejected time steps. int n_fsolve: the number of failed calls to fsolve(). real t(0:n): the sequence of time values. */ { double dt; double dt_max; double dt_min; int i; dt_max = 0.0; dt_min = HUGE_VAL; for ( i = 0; i < n - 1; i++ ) { dt = t[i+1] - t[i]; if ( dt_max < dt ) { dt_max = dt; } if ( dt < dt_min ) { dt_min = dt; } } printf ( "\n" ); printf ( "mad_stats():\n" ); printf ( " Number of timesteps = %d\n", n ); printf ( " Rejected timesteps = %d\n", n_rejected ); printf ( " Failed fsolve calls = %d\n", n_fsolve ); printf ( " First time interval: [%g,%g]\n", t[0], t[1] ); printf ( " Final time interval: [%g,%g]\n", t[n-2], t[n-1] ); printf ( " Minimum timestep is %g\n", dt_min ); printf ( " Maximum timestep is %g\n", dt_max ); return; } /******************************************************************************/ void mad_step ( int m, double to, double yo[], double tauo, void f ( double t, double y[], double dydt[] ), double theta, double tol, double *tn, double yn[], bool *fsolve_fail ) /******************************************************************************/ /* Purpose: mad_step() extends the (t,y) solution sequence by one more step. Discussion: This code is requested to add one more pair of (t,y) data to a string of estimated solutions to the differential equation. On input, an initial condition and n solutions have already been computed and stored in the arrays t and y. A stepsize tau(n) and factor theta are provided. These define a solution ym of the implicit backward Euler method at tm = t(n) + theta * tau(n). Then ym is used to construct the implicit midpoint solution y(n+1,1:n) at time t(n+1) = t(n) + tau(n). Licensing: This code is distributed under the MIT license. Modified: 15 July 2024 Author: John Burkardt, Catalin Trenchea Reference: Catalin Trenchea, John Burkardt, Refactorization of the midpoint rule, Applied Mathematics Letters, Volume 107, September 2020. Input: int m: the size of a single solution vector y. double to, yo(1:m): the "old" time and solution. double tauo: the time step. external f: the function that evaluates the right hand side of the differential equation, of the form void f ( t, y, dydt ) double theta: the theta-method parameter. Common values include: theta = 0.0 for the backward Euler method. theta = 0.5 for the midpoint method. theta = 1.0 for the forward Euler method. Output: double tn, yn(1:m): time and solution information has been added as index n+1. bool fsolve_fail: true if the call to fsolve() failed. */ { int i; int info; double *fm; double tm; double *ym; *fsolve_fail = false; /* Make a forward solution prediction at to + theta * tau */ fm = ( double * ) malloc ( m * sizeof ( double ) ); ym = ( double * ) malloc ( m * sizeof ( double ) ); tm = to + theta * tauo; f ( tm, yo, fm ); for ( i = 0; i < m; i++ ) { ym[i] = yo[i] + theta * tauo * fm[i]; } /* Starting from the forward prediction, solve the implicit backward Euler equation for ym. */ info = fsolve_be ( f, m, to, yo, tm, ym, fm, tol ); if ( info != 1 ) { printf ( "\n" ); printf ( "mad_step(): Fatal error!\n" ); printf ( " fsolve_be() returns info = %d\n", info ); *fsolve_fail = true; return; } /* Now predict the solution at tn. */ *tn = to + tauo; for ( i = 0; i < m; i++ ) { yn[i] = ( 1.0 / theta ) * ym[i] + ( 1.0 - 1.0 / theta ) * yo[i]; } free ( ym ); free ( fm ); return; } /******************************************************************************/ void mad_timestep_plot ( int n, double t[], char *label ) /******************************************************************************/ /* Purpose: mad_timestep_plot() writes data and command files for a timestep plot. Licensing: This code is distributed under the MIT license. Modified: 29 June 2024 Author: John Burkardt Input: int n: the number of steps. double t(0:n): the sequence of time values. char *label: a title. */ { char command_filename[255]; FILE *command; char data_filename[255]; FILE *data; int j; char png_filename[255]; /* Create the data file. */ strcpy ( data_filename, label ); strcat ( data_filename, "_timestamp_data.txt" ); data = fopen ( data_filename, "wt" ); for ( j = 1; j < n; j++ ) { fprintf ( data, " %d %g\n", j - 1, t[j] - t[j-1] ); } fclose ( data ); printf ( "\n" ); printf ( " Time step data stored in \"%s\".\n", data_filename ); /* Create the command file. */ strcpy ( command_filename, label ); strcat ( command_filename, "_timestep_commands.txt" ); strcpy ( png_filename, label ); strcat ( png_filename, "_timestep.png" ); 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'\n", png_filename ); fprintf ( command, "set xlabel '<-- Index -->'\n" ); fprintf ( command, "set ylabel '<-- DT -->'\n" ); fprintf ( command, "set title 'MAD Timesteps'\n" ); fprintf ( command, "set grid\n" ); fprintf ( command, "set style data lines\n" ); fprintf ( command, "plot '%s' using 1:2 with lines lw 3\n", data_filename ); fprintf ( command, "quit\n" ); fclose ( command ); printf ( " Timestep plot commands stored in \"%s\".\n", command_filename ); return; }