# include # include # include # include "stiff_ode.h" /******************************************************************************/ 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 stiff_deriv ( double t, double y[], double dydt[] ) /******************************************************************************/ /* stiff_deriv() evaluates the right hand side of the stiff ODE. Discussion: y' = lambda * ( cos(t) - y ) y(t0) = y0 Licensing: This code is distributed under the MIT license. Modified: 20 June 2025 Author: John Burkardt Input: double T, Y[1]: the time and solution value. Output: double DYDT[1]: the derivative value. */ { double lambda; stiff_parameters ( NULL, NULL, NULL, NULL, &lambda, NULL, NULL, NULL ); dydt[0] = lambda * ( cos ( t ) - y[0] ); 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) - 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 - 0 ) / ( 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] ) - y[i] ); } return; } /******************************************************************************/ void stiff_euler_backward ( int n, double t[], double y[] ) /******************************************************************************/ /* Purpose: stiff_euler_backward() uses the backward Euler method on the stiff ODE. Discussion: y' = lambda * ( cos(t) - 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 - 0 ) / ( 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+1] ) ) / ( 1.0 + lambda * dt ); } return; } /******************************************************************************/ double *stiff_exact ( int n, double t[] ) /******************************************************************************/ /* Purpose: stiff_exact() evaluates the exact solution of the stiff ODE. Discussion: y' = lambda * ( cos(t) - 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 values. double T[N]: the evaluation times. Output: double STIFF_EXACT[N]: the exact solution values. */ { int i; double lambda; double *y; stiff_parameters ( NULL, NULL, NULL, NULL, &lambda, NULL, NULL, NULL ); y = ( double * ) malloc ( n * sizeof ( double ) ); for ( i = 0; i < n; i++ ) { y[i] = lambda * ( sin ( t[i] ) + lambda * cos ( t[i] ) - lambda * exp ( - lambda * t[i] ) ) / ( lambda * lambda + 1.0 ); } return y; } /******************************************************************************/ void stiff_midpoint ( int n, double t[], double y[] ) /******************************************************************************/ /* Purpose: stiff_midpoint() uses the midpoint method on the stiff ODE. Discussion: y' = lambda * ( cos(t) - 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 - 0 ) / ( 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] + lambda * 0.5 * dt * ( cos ( t[i] ) - y[i] + cos ( t[i+1] ) ) ) / ( 1.0 + lambda * 0.5 * dt ); } return; } /******************************************************************************/ void stiff_parameters ( double *lambda_in, double *t0_in, double *y0_in, double *tstop_in, double *lambda_out, double *t0_out, double *y0_out, double *tstop_out ) /******************************************************************************/ /* Purpose: stiff_parameters() returns parameters of the stiff ODE. Discussion: y' = lambda * ( cos(t) - y ) y(t0) = y0 Licensing: This code is distributed under the MIT license. Modified: 20 June 2025 Author: John Burkardt Input: double *lambda_in, a parameter. double *t0_in, double *y0_in: the initial time and value. double *tstop_in: the final time. Output: double *lambda_out, a parameter. double *t0_out, double *y0_out: the initial time and value. double *tstop_out: the final time. */ { static double lambda_default = 50.0; static double t0_default = 0.0; static double y0_default = 0.0; static double tstop_default = 1.0; /* New values, if supplied on input, overwrite the current values. */ if ( lambda_in ) { lambda_default = *lambda_in; } if ( t0_in ) { t0_default = *t0_in; } if ( y0_in ) { y0_default = *y0_in; } if ( tstop_in ) { tstop_default = *tstop_in; } /* The current values are copied to the output. */ if ( lambda_out ) { *lambda_out = lambda_default; } if ( t0_out ) { *t0_out = t0_default; } if ( y0_out ) { *y0_out = y0_default; } if ( tstop_out ) { *tstop_out = tstop_default; } return; }