# include # include # include # include # include # include "fsolve.h" int main ( ); void fsolve_test1 ( ); void f1 ( int n, double x[], double fx[] ); void fsolve_test2 ( ); void f2 ( int n, double x[], double fx[] ); void fsolve_test3 ( ); void f3 ( int n, double x[], double fx[] ); void fsolve_test4 ( ); void f4 ( int n, double x[], double fx[] ); void predator_prey_be_test ( ); void predator_prey_tr_test ( ); void predator_prey_dydt ( double t, double y[], double dy[] ); void r8vec2_print ( int n, double a1[], double a2[], char *title ); void stiff_bdf2_test ( ); void stiff_be_test ( ); void stiff_tr_test ( ); void stiff_dydt ( double t, double y[], double dy[] ); void timestamp ( ); /******************************************************************************/ int main ( ) /******************************************************************************/ /* Purpose: fsolve_test() tests fsolve(). Licensing: This code is distributed under the MIT license. Modified: 29 November 2023 Author: John Burkardt */ { timestamp ( ); printf ( "\n" ); printf ( "fsolve_test():\n" ); printf ( " C version:\n" ); printf ( " Test fsolve().\n" ); fsolve_test1 ( ); fsolve_test2 ( ); fsolve_test3 ( ); fsolve_test4 ( ); predator_prey_be_test ( ); predator_prey_tr_test ( ); stiff_bdf2_test ( ); stiff_be_test ( ); stiff_tr_test ( ); /* Terminate. */ printf ( "\n" ); printf ( "fsolve_test()\n" ); printf ( " Normal end of execution.\n" ); printf ( "\n" ); timestamp ( ); return 0; } /******************************************************************************/ void fsolve_test1 ( ) /******************************************************************************/ /* Purpose: fsolve_test1() tests fsolve on a system of 1 equations. Discussion: This is a version of Kepler's equation. Licensing: This code is distributed under the MIT license. Modified: 08 April 2021 Author: John Burkardt */ { double *fx; int info; int n = 1; double tol = 0.00001; double *x; fx = ( double * ) malloc ( n * sizeof ( double ) ); x = ( double * ) malloc ( n * sizeof ( double ) ); printf ( "\n" ); printf ( "fsolve_test1():\n" ); printf ( " fsolve() solves a nonlinear system of 1 equation.\n" ); x[0] = 0.0; f1 ( n, x, fx ); r8vec2_print ( n, x, fx, " Initial X, F(X)" ); info = fsolve ( f1, n, x, fx, tol ); printf ( "\n" ); printf ( " Returned value of INFO = %d", info ); r8vec2_print ( n, x, fx, " Final X, FX" ); /* Free memory. */ free ( fx ); free ( x ); return; } /******************************************************************************/ void f1 ( int n, double x[], double fx[] ) /******************************************************************************/ /* Purpose: f1() evaluates a nonlinear system of 1 equation. Licensing: This code is distributed under the MIT license. Modified: 08 April 2021 Author: John Burkardt Input: int N, the number of variables. double X[N], the variable values. Output: double FX[N], the function values. */ { double e; double m; e = 0.8; m = 5.0; fx[0] = x[0] - m - e * sin ( x[0] ); return; } /******************************************************************************/ void fsolve_test2 ( ) /******************************************************************************/ /* Purpose: fsolve_test2() tests fsolve on a system of 2 equations. Discussion: The set of nonlinear equations is: x1 * x1 - 10 * x1 + x2 * x2 + 8 = 0 x1 * x2 * x2 + x1 - 10 * x2 + 8 = 0 with solution x1 = x2 = 1 Licensing: This code is distributed under the MIT license. Modified: 07 April 2021 Author: John Burkardt */ { double *fx; int info; int n = 2; double tol = 0.00001; double *x; fx = ( double * ) malloc ( n * sizeof ( double ) ); x = ( double * ) malloc ( n * sizeof ( double ) ); printf ( "\n" ); printf ( "fsolve_test2():\n" ); printf ( " fsolve() solves a nonlinear system of 2 equations.\n" ); x[0] = 3.0; x[1] = 0.0; f2 ( n, x, fx ); r8vec2_print ( n, x, fx, " Initial X, F(X)" ); info = fsolve ( f2, n, x, fx, tol ); printf ( "\n" ); printf ( " Returned value of INFO = %d", info ); r8vec2_print ( n, x, fx, " Final X, FX" ); /* Free memory. */ free ( fx ); free ( x ); return; } /******************************************************************************/ void f2 ( int n, double x[], double fx[] ) /******************************************************************************/ /* Purpose: f2() evaluates a nonlinear system of 2 equations. Licensing: This code is distributed under the MIT license. Modified: 07 April 2021 Author: John Burkardt Input: int N, the number of variables. double X[N], the variable values. Output: double FX[N], the function values. */ { fx[0] = x[0] * x[0] - 10.0 * x[0] + x[1] * x[1] + 8.0; fx[1] = x[0] * x[1] * x[1] + x[0] - 10.0 * x[1] + 8.0; return; } /******************************************************************************/ void fsolve_test3 ( ) /******************************************************************************/ /* Purpose: fsolve_test3() tests fsolve on a system of 4 equations. Licensing: This code is distributed under the MIT license. Modified: 08 April 2021 Author: John Burkardt */ { double *fx; int i; int info; int n = 4; double tol = 0.00001; double *x; fx = ( double * ) malloc ( n * sizeof ( double ) ); x = ( double * ) malloc ( n * sizeof ( double ) ); printf ( "\n" ); printf ( "fsolve_test3():\n" ); printf ( " fsolve() solves a nonlinear system of 4 equations.\n" ); for ( i = 0; i < n; i++ ) { x[i] = 0.0; } f3 ( n, x, fx ); r8vec2_print ( n, x, fx, " Initial X, F(X)" ); info = fsolve ( f3, n, x, fx, tol ); printf ( "\n" ); printf ( " Returned value of INFO = %d", info ); r8vec2_print ( n, x, fx, " Final X, FX" ); /* Free memory. */ free ( fx ); free ( x ); return; } /******************************************************************************/ void f3 ( int n, double x[], double fx[] ) /******************************************************************************/ /* Purpose: f3() evaluates a nonlinear system of 4 equations. Licensing: This code is distributed under the MIT license. Modified: 08 April 2021 Author: John Burkardt Input: int N, the number of variables. double X[N], the variable values. Output: double FX[N], the function values. */ { int i; for ( i = 0; i < n; i++ ) { fx[i] = pow ( x[i] - ( i + 1 ), 2 ); } return; } /******************************************************************************/ void fsolve_test4 ( ) /******************************************************************************/ /* Purpose: fsolve_test4() tests fsolve on a system of 8 equations. Licensing: This code is distributed under the MIT license. Modified: 08 April 2021 Author: John Burkardt */ { double *fx; int i; int info; int n = 8; double tol = 0.00001; double *x; fx = ( double * ) malloc ( n * sizeof ( double ) ); x = ( double * ) malloc ( n * sizeof ( double ) ); printf ( "\n" ); printf ( "fsolve_test4():\n" ); printf ( " fsolve() solves a nonlinear system of 8 equations.\n" ); for ( i = 0; i < n; i++ ) { x[i] = 0.0; } f4 ( n, x, fx ); r8vec2_print ( n, x, fx, " Initial X, F(X)" ); info = fsolve ( f4, n, x, fx, tol ); printf ( "\n" ); printf ( " Returned value of INFO = %d", info ); r8vec2_print ( n, x, fx, " Final X, FX" ); /* Free memory. */ free ( fx ); free ( x ); return; } /******************************************************************************/ void f4 ( int n, double x[], double fx[] ) /******************************************************************************/ /* Purpose: f4() evaluates a nonlinear system of 8 equations. Licensing: This code is distributed under the MIT license. Modified: 08 April 2021 Author: John Burkardt Input: int N, the number of variables. double X[N], the variable values. Output: double FX[N], the function values. */ { int i; for ( i = 0; i < n; i++ ) { fx[i] = ( 3.0 - 2.0 * x[i] ) * x[i] + 1.0; if ( 0 < i ) { fx[i] = fx[i] - x[i-1]; } if ( i < n - 1 ) { fx[i] = fx[i] - 2.0 * x[i+1]; } } return; } /******************************************************************************/ void predator_prey_be_test ( ) /******************************************************************************/ /* Purpose: predator_prey_be_test() tests fsolve_be() on the predator prey ODE. Licensing: This code is distributed under the MIT license. Modified: 08 November 2023 Author: John Burkardt */ { double fm[2]; int info; int m = 2; double tm; double to; double tol; double ym[2]; double yo[2]; printf ( "\n" ); printf ( "predator_prey_be_test():\n" ); to = 0.0; yo[0] = 5000.0; yo[1] = 100.0; tm = 1.0; ym[0] = yo[0]; ym[1] = yo[1]; tol = 1.0E-05; backward_euler_residual ( predator_prey_dydt, m, to, yo, tm, ym, fm ); printf ( "\n" ); printf ( " Initial ||ym-yo-(tm-to)*dydt(tm,ym)||:\n" ); printf ( "%g\n", enorm ( m, fm ) ); info = fsolve_be ( predator_prey_dydt, m, to, yo, tm, ym, fm, tol ); if ( info != 1 ) { printf ( "\n" ); printf ( " fsolve_be() returned error flag info = %d\n", info ); } backward_euler_residual ( predator_prey_dydt, m, to, yo, tm, ym, fm ); printf ( "\n" ); printf ( " Final ||ym-yo-(tm-to)*dydt(tm,ym)||:\n" ); printf ( "%g\n", enorm ( m, fm ) ); return; } /******************************************************************************/ void predator_prey_tr_test ( ) /******************************************************************************/ /* Purpose: predator_prey_tr_test() tests fsolve_tr() on the predator prey ODE. Licensing: This code is distributed under the MIT license. Modified: 15 November 2023 Author: John Burkardt */ { double ft[2]; int info; int m = 2; double tn; double to; double tol; double yn[2]; double yo[2]; printf ( "\n" ); printf ( "predator_prey_tr_test():\n" ); to = 0.0; yo[0] = 5000.0; yo[1] = 100.0; tn = 1.0; yn[0] = yo[0]; yn[1] = yo[1]; tol = 1.0E-05; trapezoidal_residual ( predator_prey_dydt, m, to, yo, tn, yn, ft ); printf ( "\n" ); printf ( " Initial residual:\n" ); printf ( "%g\n", enorm ( m, ft ) ); info = fsolve_tr ( predator_prey_dydt, m, to, yo, tn, yn, ft, tol ); if ( info != 1 ) { printf ( "\n" ); printf ( " fsolve_tr() returned error flag info = %d\n", info ); } trapezoidal_residual ( predator_prey_dydt, m, to, yo, tn, yn, ft ); printf ( "\n" ); printf ( " Final residual:\n" ); printf ( "%g\n", enorm ( m, ft ) ); return; } /******************************************************************************/ void predator_prey_dydt ( double t, double y[], double dy[] ) /******************************************************************************/ /* 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 DY[2]: the right hand side of the 2 ODE's. */ { dy[0] = 2.0 * y[0] - 0.001 * y[0] * y[1]; dy[1] = - 10.0 * y[1] + 0.002 * y[0] * y[1]; return; } /******************************************************************************/ void r8vec2_print ( int n, double a1[], double a2[], char *title ) /******************************************************************************/ /* Purpose: r8vec2_print() prints an R8VEC2. Discussion: An R8VEC2 is a dataset consisting of N pairs of real values, stored as two separate vectors A1 and A2. Licensing: This code is distributed under the MIT license. Modified: 26 March 2009 Author: John Burkardt Input: int N, the number of components of the vector. double A1[N], double A2[N], the vectors to be printed. char *TITLE, a title. */ { int i; fprintf ( stdout, "\n" ); fprintf ( stdout, "%s\n", title ); fprintf ( stdout, "\n" ); for ( i = 0; i < n; i++ ) { fprintf ( stdout, " %4d: %14g %14g\n", i, a1[i], a2[i] ); } return; } /******************************************************************************/ void stiff_bdf2_test ( ) /******************************************************************************/ /* Purpose: stiff_bdf2_test() tests fsolve_bdf2() on the stiff ODE. Licensing: This code is distributed under the MIT license. Modified: 29 November 2023 Author: John Burkardt */ { double fm[1]; int info; int m = 1; double t1; double t2; double t3; double tol; double y1[1]; double y2[1]; double y3[1]; printf ( "\n" ); printf ( "stiff_bdf2_test():\n" ); t1 = 0.0; y1[0] = 0.0; t2 = 0.25; y2[0] = 50.0 * ( sin ( t2 ) + 50.0 * cos ( t2 ) - 50.0 * exp ( - 50.0 * t2 ) ) / 2501.0; t3 = 0.5; y3[0] = y2[0]; tol = 1.0E-05; bdf2_residual ( stiff_dydt, m, t1, y1, t2, y2, t3, y3, fm ); printf ( "\n" ); printf ( " Initial residual:\n" ); printf ( "%g\n", enorm ( m, fm ) ); info = fsolve_bdf2 ( stiff_dydt, m, t1, y1, t2, y2, t3, y3, fm, tol ); if ( info != 1 ) { printf ( "\n" ); printf ( " fsolve_bdf2() returned error flag info = %d\n", info ); } bdf2_residual ( stiff_dydt, m, t1, y1, t2, y2, t3, y3, fm ); printf ( "\n" ); printf ( " Final residual:\n" ); printf ( "%g\n", enorm ( m, fm ) ); return; } /******************************************************************************/ void stiff_be_test ( ) /******************************************************************************/ /* Purpose: stiff_be_test() tests fsolve_be() on the stiff ODE. Licensing: This code is distributed under the MIT license. Modified: 29 November 2023 Author: John Burkardt */ { double fm[1]; int info; int m = 1; double tm; double to; double tol; double ym[1]; double yo[1]; printf ( "\n" ); printf ( "stiff_be_test():\n" ); to = 0.0; yo[0] = 0.0; tm = 0.25; ym[0] = yo[0]; tol = 1.0E-05; backward_euler_residual ( stiff_dydt, m, to, yo, tm, ym, fm ); printf ( "\n" ); printf ( " Initial ||ym-yo-(tm-to)*dydt(tm,ym)||:\n" ); printf ( "%g\n", enorm ( m, fm ) ); info = fsolve_be ( stiff_dydt, m, to, yo, tm, ym, fm, tol ); if ( info != 1 ) { printf ( "\n" ); printf ( " fsolve_be() returned error flag info = %d\n", info ); } backward_euler_residual ( stiff_dydt, m, to, yo, tm, ym, fm ); printf ( "\n" ); printf ( " Final ||ym-yo-(tm-to)*dydt(tm,ym)||:\n" ); printf ( "%g\n", enorm ( m, fm ) ); return; } /******************************************************************************/ void stiff_tr_test ( ) /******************************************************************************/ /* Purpose: stiff_tr_test() tests fsolve_tr() on the stiff ODE. Licensing: This code is distributed under the MIT license. Modified: 29 November 2023 Author: John Burkardt */ { double ft[1]; int info; int m = 1; double tn; double to; double tol; double yn[1]; double yo[1]; printf ( "\n" ); printf ( "stiff_tr_test():\n" ); to = 0.0; yo[0] = 0.0; tn = 0.25; yn[0] = yo[0]; tol = 1.0E-05; trapezoidal_residual ( stiff_dydt, m, to, yo, tn, yn, ft ); printf ( "\n" ); printf ( " Initial residual:\n" ); printf ( "%g\n", enorm ( m, ft ) ); info = fsolve_tr ( stiff_dydt, m, to, yo, tn, yn, ft, tol ); if ( info != 1 ) { printf ( "\n" ); printf ( " fsolve_tr() returned error flag info = %d\n", info ); } trapezoidal_residual ( stiff_dydt, m, to, yo, tn, yn, ft ); printf ( "\n" ); printf ( " Final residual:\n" ); printf ( "%g\n", enorm ( m, ft ) ); 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 DY[1]: the derivative value. */ { dy[0] = 50.0 * ( cos ( t ) - y[0] ); return; } /******************************************************************************/ void timestamp ( ) /******************************************************************************/ /* Purpose: timestamp() prints the current YMDHMS date as a time stamp. Example: 31 May 2001 09:45:54 AM Licensing: This code is distributed under the MIT license. Modified: 24 September 2003 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 ); fprintf ( stdout, "%s\n", time_buffer ); return; # undef TIME_SIZE }