# include <cmath>
# include <complex>
# include <cstdlib>
# include <ctime>
# include <iomanip>
# include <iostream>

using namespace std;

# include "fsolve.hpp"

int main ( );
void fsolve_test1 ( );
void f1 ( int n, double x[], double fvec[] );
void fsolve_test2 ( );
void f2 ( int n, double x[], double fvec[] );
void fsolve_test3 ( );
void f3 ( int n, double x[], double fvec[] );
void fsolve_test4 ( );
void f4 ( int n, double x[], double fvec[] );
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[], string title );
void stiff_bdf2_test ( );
void stiff_be_test ( );
void stiff_tr_test ( );
void stiff_dydt ( double t, double y[], double dy[] );
void timestamp ( );

//****************************************************************************80

int main ( )

//****************************************************************************80
//
//  Purpose:
//
//    fsolve_test() tests fsolve().
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    29 November 2023
//
//  Author:
//
//    John Burkardt
//
{
  timestamp ( );
  cout << "\n";
  cout << "fsolve_test()\n";
  cout << "  C++ version:\n";
  cout << "  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.
//
  cout << "\n";
  cout << "fsolve_test():\n";
  cout << "  Normal end of execution.\n";
  cout << "\n";
  timestamp ( );

  return 0;
}
//****************************************************************************80

void fsolve_test1 ( )

//****************************************************************************80
//
//  Purpose:
//
//    fsolve_test1() tests fsolve() with a system of 1 equation.
//
//  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 lwa;
  int n = 1;
  double tol = 0.00001;
  double *wa;
  double *x;

  lwa = ( n * ( 3 * n + 13 ) ) / 2;
  fx = new double[n];
  wa = new double[lwa];
  x = new double[n];

  cout << "\n";
  cout << "fsolve_test1():\n";
  cout << "  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, wa, lwa );

  cout << "\n";
  cout << "  Returned value of INFO = " << info << "\n";
  r8vec2_print ( n, x, fx, "  Final X, FX" );
//
//  Free memory.
//
  delete [] fx;
  delete [] wa;
  delete [] x;

  return;
}
//****************************************************************************80

void f1 ( int n, double x[], double fx[] )

//****************************************************************************80
//
//  Purpose:
//
//    f1() evaluates a system of 1 nonlinear equation.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    08 April 2010
//
//  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;
}
//****************************************************************************80

void fsolve_test2 ( )

//****************************************************************************80
//
//  Purpose:
//
//    fsolve_test2() tests fsolve() with a system of 2 equations.
//
//  Discussion:
//
//    This is an example of what your main program would look
//    like if you wanted to use MINPACK to solve N nonlinear equations
//    in N unknowns.  In this version, we avoid computing the jacobian
//    matrix, and request that MINPACK approximate it for us.
//
//    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:
//
//    08 April 2021
//
//  Author:
//
//    John Burkardt
//
{
  double *fx;
  int info;
  int lwa;
  int n = 2;
  double tol = 0.00001;
  double *wa;
  double *x;

  lwa = ( n * ( 3 * n + 13 ) ) / 2;
  fx = new double[n];
  wa = new double[lwa];
  x = new double[n];

  cout << "\n";
  cout << "fsolve_test2():\n";
  cout << "  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, wa, lwa );

  cout << "\n";
  cout << "  Returned value of INFO = " << info << "\n";
  r8vec2_print ( n, x, fx, "  Final X, FX" );
//
//  Free memory.
//
  delete [] fx;
  delete [] wa;
  delete [] x;

  return;
}
//****************************************************************************80

void f2 ( int n, double x[], double fx[] )

//****************************************************************************80
//
//  Purpose:
//
//    f2() evaluates a system of 2 nonlinear equations.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    08 April 2010
//
//  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;
}
//****************************************************************************80

void fsolve_test3 ( )

//****************************************************************************80
//
//  Purpose:
//
//    fsolve_test3() tests fsolve() with 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 lwa;
  int n = 4;
  double tol = 0.00001;
  double *wa;
  double *x;

  lwa = ( n * ( 3 * n + 13 ) ) / 2;
  fx = new double[n];
  wa = new double[lwa];
  x = new double[n];

  cout << "\n";
  cout << "fsolve_test3():\n";
  cout << "  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, wa, lwa );

  cout << "\n";
  cout << "  Returned value of INFO = " << info << "\n";
  r8vec2_print ( n, x, fx, "  Final X, FX" );
//
//  Free memory.
//
  delete [] fx;
  delete [] wa;
  delete [] x;

  return;
}
//****************************************************************************80

void f3 ( int n, double x[], double fx[] )

//****************************************************************************80
//
//  Purpose:
//
//    f3() evaluates a system of 4 nonlinear equations.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    08 April 2010
//
//  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;
}
//****************************************************************************80

void fsolve_test4 ( )

//****************************************************************************80
//
//  Purpose:
//
//    fsolve_test4() tests fsolve() with 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 lwa;
  int n = 8;
  double tol = 0.00001;
  double *wa;
  double *x;

  lwa = ( n * ( 3 * n + 13 ) ) / 2;
  fx = new double[n];
  wa = new double[lwa];
  x = new double[n];

  cout << "\n";
  cout << "fsolve_test4():\n";
  cout << "  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, wa, lwa );

  cout << "\n";
  cout << "  Returned value of INFO = " << info << "\n";
  r8vec2_print ( n, x, fx, "  Final X, FX" );
//
//  Free memory.
//
  delete [] fx;
  delete [] wa;
  delete [] x;

  return;
}
//****************************************************************************80

void f4 ( int n, double x[], double fx[] )

//****************************************************************************80
//
//  Purpose:
//
//    f4() evaluates a system of 8 nonlinear equations.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    08 April 2010
//
//  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;
}
//****************************************************************************80

void predator_prey_be_test ( )

//****************************************************************************80
//
//  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];

  cout << "\n";
  cout << "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 );
  cout << "\n";
  cout << "  Initial ||ym-yo-(tm-to)*dydt(tm,ym)||:\n";
  cout << enorm ( m, fm ) << "\n";

  info = fsolve_be ( predator_prey_dydt, m, to, yo, tm, ym, fm, tol );
  if ( info != 1 )
  {
    cout << "\n";
    cout << "  fsolve_be() returned error flag info = " << info << "\n";
  }

  backward_euler_residual ( predator_prey_dydt, m, to, yo, tm, ym, fm );
  cout << "\n";
  cout << "  Final ||ym-yo-(tm-to)*dydt(tm,ym)||:\n";
  cout << enorm ( m, fm ) << "\n";

  return;
}
//****************************************************************************80

void predator_prey_tr_test ( )

//****************************************************************************80
//
//  Purpose:
//
//    predator_prey_tr_test() tests fsolve_tr() on the predator prey ODE.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    08 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];

  cout << "\n";
  cout << "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 );
  cout << "\n";
  cout << "  Initial residual:\n";
  cout << enorm ( m, ft ) << "\n";

  info = fsolve_tr ( predator_prey_dydt, m, to, yo, tn, yn, ft, tol );
  if ( info != 1 )
  {
    cout << "\n";
    cout << "  fsolve_tr() returned error flag info = " << info << "\n";
  }

  trapezoidal_residual ( predator_prey_dydt, m, to, yo, tn, yn, ft );
  cout << "\n";
  cout << "  Final residual:\n";
  cout << enorm ( m, ft ) << "\n";

  return;
}
//****************************************************************************80

void predator_prey_dydt ( double t, double y[], double dy[] )

//****************************************************************************80
//
//  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;
}
//****************************************************************************80

void r8vec2_print ( int n, double a1[], double a2[], string title )

//****************************************************************************80
//
//  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:
//
//    19 November 2002
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    int N, the number of components of the vector.
//
//    double A1[N], double A2[N], the vectors to be printed.
//
//    string TITLE, a title.
//
{
  int i;

  cout << "\n";
  cout << title << "\n";
  cout << "\n";
  for ( i = 0; i <= n - 1; i++ )
  {
    cout << setw(6)  << i
         << ": " << setw(14) << a1[i]
         << "  " << setw(14) << a2[i] << "\n";
  }

  return;
}
//****************************************************************************80

void stiff_bdf2_test ( )

//****************************************************************************80
//
//  Purpose:
//
//    stiff_bdf2_test() tests fsolve_bdf2() on the stiff ODE.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    18 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];

  cout << "\n";
  cout << "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.50;
  y3[0] = y2[0];
  tol = 1.0E-05;

  bdf2_residual ( stiff_dydt, m, t1, y1, t2, y2, t3, y3, fm );
  cout << "\n";
  cout << "  Initial residual:\n";
  cout << enorm ( m, fm ) << "\n";

  info = fsolve_bdf2 ( stiff_dydt, m, t1, y1, t2, y2, t3, y3, fm, tol );
  if ( info != 1 )
  {
    cout << "\n";
    cout << "  fsolve_bdf2() returned error flag info = " << info << "\n";
  }

  bdf2_residual ( stiff_dydt, m, t1, y1, t2, y2, t3, y3, fm );
  cout << "\n";
  cout << "  Final residual:\n";
  cout << enorm ( m, fm ) << "\n";

  return;
}
//****************************************************************************80

void stiff_be_test ( )

//****************************************************************************80
//
//  Purpose:
//
//    stiff_be_test() tests fsolve_be() on the stiff ODE.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    08 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];

  cout << "\n";
  cout << "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 );
  cout << "\n";
  cout << "  Initial ||ym-yo-(tm-to)*dydt(tm,ym)||:\n";
  cout << enorm ( m, fm ) << "\n";

  info = fsolve_be ( stiff_dydt, m, to, yo, tm, ym, fm, tol );
  if ( info != 1 )
  {
    cout << "\n";
    cout << "  fsolve_be() returned error flag info = " << info << "\n";
  }

  backward_euler_residual ( stiff_dydt, m, to, yo, tm, ym, fm );
  cout << "\n";
  cout << "  Final ||ym-yo-(tm-to)*dydt(tm,ym)||:\n";
  cout << enorm ( m, fm ) << "\n";

  return;
}
//****************************************************************************80

void stiff_tr_test ( )

//****************************************************************************80
//
//  Purpose:
//
//    stiff_tr_test() tests fsolve_tr() on the stiff ODE.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    16 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];

  cout << "\n";
  cout << "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 );
  cout << "\n";
  cout << "  Initial residual:\n";
  cout << enorm ( m, ft ) << "\n";

  info = fsolve_tr ( stiff_dydt, m, to, yo, tn, yn, ft, tol );
  if ( info != 1 )
  {
    cout << "\n";
    cout << "  fsolve_tr() returned error flag info = " << info << "\n";
  }

  trapezoidal_residual ( stiff_dydt, m, to, yo, tn, yn, ft );
  cout << "\n";
  cout << "  Final residual:\n";
  cout << enorm ( m, ft ) << "\n";

  return;
}
//****************************************************************************80

void stiff_dydt ( double t, double y[], double dy[] )

//****************************************************************************80
//
//  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;
}
//****************************************************************************80

void timestamp ( )

//****************************************************************************80
//
//  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:
//
//    08 July 2009
//
//  Author:
//
//    John Burkardt
//
{
# define TIME_SIZE 40

  static char time_buffer[TIME_SIZE];
  const struct std::tm *tm_ptr;
  std::time_t now;

  now = std::time ( NULL );
  tm_ptr = std::localtime ( &now );

  std::strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm_ptr );

  std::cout << time_buffer << "\n";

  return;
# undef TIME_SIZE
}