ordinary_differential_equations.cpp

/****************************************************************************************************************/
/*                                                                                                              */
/*   OpenNN: Open Neural Networks Library                                                                       */
/*   www.artelnics.com/opennn                                                                                   */
/*                                                                                                              */
/*   O R D I N A R Y   D I F F E R E N T I A L   E Q U A T I O N S   C L A S S                                  */
/*                                                                                                              */
/*   Roberto Lopez                                                                                              */ 
/*   Artelnics - Making intelligent use of data                                                                 */
/*   [email protected]                                                                                 */
/*                                                                                                              */
/****************************************************************************************************************/

// OpenNN includes

#include "ordinary_differential_equations.h"

namespace OpenNN
{

// GENERAL CONSTRUCTOR


OrdinaryDifferentialEquations::OrdinaryDifferentialEquations(void) : MathematicalModel()
{                                            
   set_default();
}


// XML CONSTRUCTOR


OrdinaryDifferentialEquations::OrdinaryDifferentialEquations(const tinyxml2::XMLDocument& ordinary_differential_equations_document)
: MathematicalModel(ordinary_differential_equations_document)
{                
   set_default();
}


// FILE CONSTRUCTOR


OrdinaryDifferentialEquations::OrdinaryDifferentialEquations(const std::string& file_name)
: MathematicalModel(file_name)
{                
   set_default();
}


// COPY CONSTRUCTOR


OrdinaryDifferentialEquations::OrdinaryDifferentialEquations(const OrdinaryDifferentialEquations& other_ordinary_differential_equations) 
 : MathematicalModel()
{
   set(other_ordinary_differential_equations);
}


// DESTRUCTOR


OrdinaryDifferentialEquations::~OrdinaryDifferentialEquations(void)
{ 
}


// ASSIGNMENT OPERATOR

// OrdinaryDifferentialEquations& operator = (const OrdinaryDifferentialEquations&) method


OrdinaryDifferentialEquations& OrdinaryDifferentialEquations::operator = (const OrdinaryDifferentialEquations& other_ordinary_differential_equations)
{
   if(this != &other_ordinary_differential_equations) 
   {
      initial_independent_variable = other_ordinary_differential_equations.initial_independent_variable;
      final_independent_variable = other_ordinary_differential_equations.final_independent_variable;
      initial_dependent_variables = other_ordinary_differential_equations.initial_dependent_variables;
   }

   return(*this);
}


// EQUAL TO OPERATOR

// bool operator == (const OrdinaryDifferentialEquations&) const method


bool OrdinaryDifferentialEquations::operator == (const OrdinaryDifferentialEquations& other_ordinary_differential_equations) const
{
   if(initial_independent_variable == other_ordinary_differential_equations.initial_independent_variable
   && final_independent_variable == other_ordinary_differential_equations.final_independent_variable
   && initial_dependent_variables == other_ordinary_differential_equations.initial_dependent_variables)
   {
      return(true);
   }
   else
   {
      return(false);
   }
}


// METHODS

// const double& get_initial_independent_variable(void) const method


const double& OrdinaryDifferentialEquations::get_initial_independent_variable(void) const
{
   return(initial_independent_variable);
}


// const double& get_final_independent_variable(void) const method


const double& OrdinaryDifferentialEquations::get_final_independent_variable(void) const
{
   return(final_independent_variable);
}


// const Vector<double>& get_initial_dependent_variables(void) const method


const Vector<double>& OrdinaryDifferentialEquations::get_initial_dependent_variables(void) const
{
   return(initial_dependent_variables);
}


// const double& get_initial_dependent_variable(const size_t&) const method


const double& OrdinaryDifferentialEquations::get_initial_dependent_variable(const size_t& i) const
{
   return(initial_dependent_variables[i]);
}


// const SolutionMethod& get_solution_method(void) const method


const OrdinaryDifferentialEquations::SolutionMethod& OrdinaryDifferentialEquations::get_solution_method(void) const
{
   return(solution_method);
}


// std::string write_solution_method(void) const method


std::string OrdinaryDifferentialEquations::write_solution_method(void) const
{
   if(solution_method == RungeKutta)
   {
      return("RungeKutta");
   }
   else if(solution_method == RungeKuttaFehlberg)
   {
      return("RungeKuttaFehlberg");
   }
   else
   {
      std::ostringstream buffer;

      buffer << "OpenNN Exception: OrdinaryDifferentialEquations class.\n"
             << "std::string write_solution_method(void) const method.\n"
             << "Unknown solution method.\n";
 
      throw std::logic_error(buffer.str());
   }
}


// const size_t& get_points_number(void) const method


const size_t& OrdinaryDifferentialEquations::get_points_number(void) const
{
   return(points_number);
}


// const double& get_tolerance(void) const method


const double& OrdinaryDifferentialEquations::get_tolerance(void) const
{
   return(tolerance);
}


// const size_t& get_initial_size(void) const method


const size_t& OrdinaryDifferentialEquations::get_initial_size(void) const
{
   return(initial_size);
}


// const size_t& get_warning_size(void) const method


const size_t& OrdinaryDifferentialEquations::get_warning_size(void) const
{
   return(warning_size);
}


// const size_t& get_error_size(void) const method


const size_t& OrdinaryDifferentialEquations::get_error_size(void) const
{
   return(error_size);
}


// void set(const OrdinaryDifferentialEquations&) method


void OrdinaryDifferentialEquations::set(const OrdinaryDifferentialEquations& other_ordinary_differential_equations)
{
   independent_variables_number = other_ordinary_differential_equations.independent_variables_number;
   dependent_variables_number = other_ordinary_differential_equations.dependent_variables_number;

   initial_independent_variable = other_ordinary_differential_equations.initial_independent_variable;
   final_independent_variable = other_ordinary_differential_equations.final_independent_variable;

   initial_dependent_variables = other_ordinary_differential_equations.initial_dependent_variables;

   solution_method = other_ordinary_differential_equations.solution_method;

   points_number = other_ordinary_differential_equations.points_number;

   tolerance = other_ordinary_differential_equations.tolerance;

   initial_size = other_ordinary_differential_equations.initial_size;
   warning_size = other_ordinary_differential_equations.warning_size;
   error_size = other_ordinary_differential_equations.error_size;

   display = other_ordinary_differential_equations.display;
}


// void set_initial_independent_variable(const double&) method


void OrdinaryDifferentialEquations::set_initial_independent_variable(const double& new_initial_independent_variable)
{
   initial_independent_variable = new_initial_independent_variable;
}


// void set_final_independent_variable(const double&) method


void OrdinaryDifferentialEquations::set_final_independent_variable(const double& new_final_independent_variable)
{
    final_independent_variable = new_final_independent_variable;
}


// void set_initial_dependent_variables(const Vector<double>&) const method


void OrdinaryDifferentialEquations::set_initial_dependent_variables(const Vector<double>& new_initial_dependent_variables) 
{
   initial_dependent_variables = new_initial_dependent_variables;
}


// void set_initial_dependent_variable(const size_t&, const double&) const method


void OrdinaryDifferentialEquations::set_initial_dependent_variable(const size_t& i, const double& new_initial_dependent_variable)
{
    initial_dependent_variables[i] = new_initial_dependent_variable;
}


// void set_solution_method(const SolutionMethod&) method


void OrdinaryDifferentialEquations::set_solution_method(const SolutionMethod& new_solution_method)
{
   solution_method = new_solution_method;
}


// void set_solution_method(const std::string&) method


void OrdinaryDifferentialEquations::set_solution_method(const std::string& new_solution_method)
{
   if(new_solution_method == "RungeKutta")
   {
      set_solution_method(RungeKutta);
   }
   else if(new_solution_method == "RungeKuttaFehlberg")
   {
      set_solution_method(RungeKuttaFehlberg);
   }
   else
   {
      std::ostringstream buffer;

      buffer << "OpenNN Exception: OrdinaryDifferentialEquations class.\n"
             << "void set_solution_method(const std::string&) method.\n"
             << "Unknown solution method: " << new_solution_method << ".\n";

      throw std::logic_error(buffer.str());
   }
}


// void set_points_number(const size_t&) method


void OrdinaryDifferentialEquations::set_points_number(const size_t& new_points_number)
{
   points_number = new_points_number;
}


// void set_tolerance(const double&) method


void OrdinaryDifferentialEquations::set_tolerance(const double& new_tolerance)
{
   tolerance = new_tolerance;
}


// void set_initial_size(const size_t&) method


void OrdinaryDifferentialEquations::set_initial_size(const size_t& new_initial_size)
{
   initial_size = new_initial_size;
}


// void set_warning_size(const size_t&) method


void OrdinaryDifferentialEquations::set_warning_size(const size_t& new_warning_size)
{
   warning_size = new_warning_size;
}


// void set_error_size(const size_t&) method


void OrdinaryDifferentialEquations::set_error_size(const size_t& new_error_size)
{
   error_size = new_error_size;
}


// void set_default(void) method


void OrdinaryDifferentialEquations::set_default(void)
{
   independent_variables_number = 1;
   dependent_variables_number = 0;

   solution_method = RungeKuttaFehlberg;

   points_number = 101;

   tolerance = 1.0e-6;

   initial_size = (size_t)1.0e3;
   warning_size = (size_t)1.0e6;
   error_size = (size_t)1.0e9;

   display = true;

}


// Matrix<double> calculate_Runge_Kutta_solution(const NeuralNetwork& neural_network, const double& xi, const double& xf, const Vector<double>& yi) const method


Matrix<double> OrdinaryDifferentialEquations::calculate_Runge_Kutta_solution(const NeuralNetwork& neural_network) const
{
   const size_t variables_number = count_variables_number();

   const double h = (final_independent_variable - initial_independent_variable)/(points_number-1.0);      

   // Fourth order Runge-Kutta coefficients

   Vector< Vector<double> > c(4);

   Matrix<double> solution(points_number, variables_number);

   Vector<double> variables(variables_number);

   // Initial variables

   solution(0,0) = initial_independent_variable;

   for(size_t j = 0; j < dependent_variables_number; j++)
   {
      solution(0,1+j) = initial_dependent_variables[j];
   }

   // Main loop

   for(size_t i = 0; i < points_number-1; i++)
   {
      solution(i+1,0) = solution(i,0) + h;

      // First coefficient 

      variables[0] = solution(i,0); 

      for(size_t j = 0; j < dependent_variables_number; j++)
      {
         variables[1+j] = solution(i,1+j);
      }   

      c[0] = calculate_dependent_variables_dots(neural_network, variables);

      // Second coefficient 

      variables[0] = solution(i,0) + h/2.0; 

      for(size_t j = 0; j < dependent_variables_number; j++)
      {
         variables[1+j] = solution(i,1+j) + h*c[0][j]/2.0;
      }   

      c[1] = calculate_dependent_variables_dots(neural_network, variables);

      // Third coefficient

      variables[0] = solution(i,0) + h/2.0; 
               
      for(size_t j = 0; j < dependent_variables_number; j++)
      {
         variables[1+j] = solution(i,1+j) + h*c[1][j]/2.0;
      }
           
      c[2] = calculate_dependent_variables_dots(neural_network, variables);

      // Fourth coefficient

      variables[0] = solution(i+1,0);

      for(size_t j = 0; j < dependent_variables_number; j++)
      {
         variables[1+j] = solution(i,1+j) + h*c[2][j]/2.0;
      }
            
      c[3] = calculate_dependent_variables_dots(neural_network, variables);

      // Dependent variables

      for(size_t j = 0; j < dependent_variables_number; j++)
      {
         solution(i+1,1+j) = solution(i,1+j) + h*(c[0][j] + 2.0*c[1][j] + 2.0*c[2][j] + c[3][j])/6.0;
      }
   }

   solution(points_number-1,0) = final_independent_variable;

   return(solution);
}


// Matrix<double> calculate_Runge_Kutta_final_solution(const NeuralNetwork& neural_network) const method


Vector<double> OrdinaryDifferentialEquations::calculate_Runge_Kutta_final_solution(const NeuralNetwork& neural_network) const
{
   const size_t variables_number = count_variables_number();

   const double h = (final_independent_variable - initial_independent_variable)/(points_number-1.0);      

   // Fourth order Runge-Kutta coefficients

   Vector< Vector<double> > c(4);

   Vector<double> final_solution(variables_number);

   Vector<double> variables(variables_number);

   // Initial variables

   final_solution[0] = initial_independent_variable;

   for(size_t j = 0; j < dependent_variables_number; j++)
   {
      final_solution[0] = initial_dependent_variables[j];
   }

   // Main loop

   for(size_t i = 0; i < points_number-1; i++)
   {
      final_solution[0] = final_solution[0] + h;

      // First coefficient 

      variables[0] = final_solution[0]; 

      for(size_t j = 0; j < dependent_variables_number; j++)
      {
         variables[1+j] = final_solution[1+j]; 
      }   

      c[0] = calculate_dependent_variables_dots(neural_network, variables);

      // Second coefficient 

      variables[0] = final_solution[0] + h/2.0; 

      for(size_t j = 0; j < dependent_variables_number; j++)
      {
         variables[1+j] = final_solution[1+j] + h*c[0][j]/2.0; 
      }   

      c[1] = calculate_dependent_variables_dots(neural_network, variables);

      // Third coefficient

      variables[0] = final_solution[0] + h/2.0; 
               
      for(size_t j = 0; j < dependent_variables_number; j++)
      {
         variables[1+j] = final_solution[1+j] + h*c[1][j]/2.0; 
      }
           
      c[2] = calculate_dependent_variables_dots(neural_network, variables);

      // Fourth coefficient

      variables[0] = final_solution[0]; 

      for(size_t j = 0; j < dependent_variables_number; j++)
      {
         variables[1+j] = final_solution[1+j] + h*c[2][j]/2.0; 
      }
            
      c[3] = calculate_dependent_variables_dots(neural_network, variables);

      // Dependent variables

      for(size_t j = 0; j < dependent_variables_number; j++)
      {
         final_solution[1+j] = final_solution[1+j] + h*(c[0][j] + 2.0*c[1][j] + 2.0*c[2][j] + c[3][j])/6.0;
      }
   }

   return(final_solution);
}


// Matrix<double> calculate_Runge_Kutta_Fehlberg_solution(const NeuralNetwork&) const method


Matrix<double> OrdinaryDifferentialEquations::calculate_Runge_Kutta_Fehlberg_solution(const NeuralNetwork& neural_network) const
{
   const double epsilon = 1.0e-12;//std::numeric_limits<double>::epsilon();

   const double a2 = 1.0/5.0;
   const double a3 = 3.0/10.0;
   const double a4 = 3.0/5.0;
   const double a5 = 1.0;
   const double a6 = 7.0/8.0;

   const double b21 = 1.0/5.0;   
   const double b31 = 3.0/40.0;
   const double b32 = 9.0/40.0;   
   const double b41 = 3.0/10.0;
   const double b42 = -9.0/10.0;   
   const double b43 = 6.0/5.0; 
   const double b51 = -11.0/54.0;
   const double b52 = 5.0/2.0;
   const double b53 = -70.0/27.0;
   const double b54 = 35.0/27.0; 
   const double b61 = 1631.0/55296.0;
   const double b62 = 175.0/512.0;   
   const double b63 = 575.0/13824.0;   
   const double b64 = 44275.0/110592.0;   
   const double b65 = 253.0/4096.0;

   const double c41 = 37.0/378.0;
   const double c42 = 0.0;
   const double c43 = 250.0/621.0;
   const double c44 = 125.0/594.0;
   const double c45 = 0.0;
   const double c46 = 512.0/1771.0;

   const double c51 = 2825.0/27648.0;  
   const double c52 = 0.0;  
   const double c53 = 18575.0/48384.0;   
   const double c54 = 13525.0/55296.0;
   const double c55 = 277.0/14336.0;
   const double c56 = 1.0/4.0;

   const double d1 = c41 - c51;
   const double d2 = c42 - c52;
   const double d3 = c43 - c53;
   const double d4 = c44 - c54;
   const double d5 = c45 - c55;
   const double d6 = c46 - c56;

   const size_t variables_number = count_variables_number();

   size_t size = initial_size;

   Vector<double> errors(dependent_variables_number, 0.0);

   double error = 0.0;

   Vector< Vector<double> > c(6);

   double hmin = 0.0;
   double h = (final_independent_variable - initial_independent_variable)*1.0e-3;

   Vector<double> variables(variables_number);

   Matrix<double> solution(size, variables_number);

   size_t point_index = 0;

   // Initial values

   solution(0,0) = initial_independent_variable;

   for(size_t j = 0; j < dependent_variables_number; j++)
   {
      solution(0,1+j) = initial_dependent_variables[j];
   } 

   // Main loop

   while(solution(point_index,0) < final_independent_variable)
   {   
      // Set smallest allowable stepsize

      hmin = 32.0*epsilon*fabs(solution(point_index,0));

      if(h < hmin)
      {
         if(display)
         {
            std::cout << "OpenNN Warning: OrdinaryDifferentialEquations class.\n" 
                      << "calculate_Runge_Kutta_Fehlberg_solution() method.\n" 
                      << "Step size is less than smallest allowable." << std::endl;
         }

         h = hmin;
      }

      if(solution(point_index,0) + h > final_independent_variable)
      {        
         h = final_independent_variable - solution(point_index,0);
      }

      // First coefficient

      variables[0] = solution(point_index,0);

      for(size_t j = 0; j < dependent_variables_number; j++)
      {
         variables[1+j] = solution(point_index,1+j);
      }

      c[0] = calculate_dependent_variables_dots(neural_network, variables)*h;

      // Second coefficient

      variables[0] = solution(point_index,0) + a2*h;

      for(size_t j = 0; j < dependent_variables_number; j++)
      {
         variables[1+j] = solution(point_index,1+j)+ b21*c[0][j];
      }

      c[1] = calculate_dependent_variables_dots(neural_network, variables)*h;

      // Third coefficient

      variables[0] = solution(point_index,0) + a3*h;

      for(size_t j = 0; j < dependent_variables_number; j++)
      {
         variables[1+j] =  solution(point_index,1+j) + b31*c[0][j] + b32*c[1][j];
      }

      c[2] = calculate_dependent_variables_dots(neural_network, variables)*h;

      // Fourth coefficient

      variables[0] = solution(point_index,0) + a4*h;

      for(size_t j = 0; j < dependent_variables_number; j++)
      {
         variables[1+j] = solution(point_index,1+j) + b41*c[0][j] + b42*c[1][j] + b43*c[2][j];
      }

      c[3] = calculate_dependent_variables_dots(neural_network, variables)*h; 

      // Fifth coefficient

      variables[0] = solution(point_index,0) + a5*h;

      for(size_t j = 0; j < dependent_variables_number; j++)
      {
         variables[1+j] = solution(point_index,1+j) + b51*c[0][j] + b52*c[1][j] + b53*c[2][j] + b54*c[3][j];
      }
      
      c[4] = calculate_dependent_variables_dots(neural_network, variables)*h;
            
      // Sixth coefficient

      variables[0] = solution(point_index,0) + a6*h;
            
      for(size_t j = 0; j < dependent_variables_number; j++)
      {
         variables[1+j] = solution(point_index,1+j) +  b61*c[0][j] + b62*c[1][j] + b63*c[2][j] + b64*c[3][j] + b65*c[4][j];
      }
      
      c[5] = calculate_dependent_variables_dots(neural_network, variables)*h;

      // Error estimate

      for(size_t j = 0; j < dependent_variables_number; j++)
      {
         errors[j] = fabs(d1*c[0][j] + d2*c[1][j] + d3*c[2][j] + d4*c[3][j] + d5*c[4][j] + d6*c[5][j]);         
      }

      error = errors.calculate_maximum();
   
      if(error <= tolerance)
      {
         solution(point_index+1,0) = solution(point_index,0) + h;

         for(size_t j = 0; j < dependent_variables_number; j++)
         {
            solution(point_index+1,1+j) = solution(point_index,1+j) + c51*c[0][j] + c52*c[1][j] + c53*c[2][j] + c54*c[3][j] + c55*c[4][j] + c56*c[5][j];
         }

         point_index++;

         if(error != 0.0)
         {
            h *= 0.9*pow(fabs(tolerance/error), 0.2);
         }

         if(point_index >= size)
         {
            size *= 2;

            if(display && size > warning_size)
            {
               std::cout << "OpenNN Warning: OrdinaryDifferentialEquations class." << std::endl
                         << "calculate_Runge_Kutta_Fehlberg_solution() method." << std::endl
                         << "Solution size is " << size << std::endl;
            }
            else if(size > error_size)
            {
               std::ostringstream buffer;

               buffer << "OpenNN Exception: OrdinaryDifferentialEquations class." << std::endl
                      << "calculate_Runge_Kutta_Fehlberg_solution() method." << std::endl
                      << "Solution size is bigger than greatest allowable." << std::endl;
                            
               throw std::logic_error(buffer.str());          
            }

            // @todo Substitute resize method

            //solution.resize(size, variables_number);
         }
      }
      else
      {
         h *= 0.9*pow(fabs(tolerance/error), 0.25);
      }
   } // end while loop   

   // @todo Substitute resize method

   //solution.resize(point_index+1, variables_number);

   return(solution);
}


// Vector<double> calculate_Runge_Kutta_Fehlberg_final_solution(const NeuralNetwork&) const method


Vector<double> OrdinaryDifferentialEquations::calculate_Runge_Kutta_Fehlberg_final_solution(const NeuralNetwork& neural_network) const
{
   const double epsilon = 1.0e-12;//std::numeric_limits<double>::epsilon();

   const double a2 = 1.0/5.0;
   const double a3 = 3.0/10.0;
   const double a4 = 3.0/5.0;
   const double a5 = 1.0;
   const double a6 = 7.0/8.0;

   const double b21 = 1.0/5.0;   
   const double b31 = 3.0/40.0;
   const double b32 = 9.0/40.0;   
   const double b41 = 3.0/10.0;
   const double b42 = -9.0/10.0;   
   const double b43 = 6.0/5.0; 
   const double b51 = -11.0/54.0;
   const double b52 = 5.0/2.0;
   const double b53 = -70.0/27.0;
   const double b54 = 35.0/27.0; 
   const double b61 = 1631.0/55296.0;
   const double b62 = 175.0/512.0;   
   const double b63 = 575.0/13824.0;   
   const double b64 = 44275.0/110592.0;   
   const double b65 = 253.0/4096.0;

   const double c41 = 37.0/378.0;
   const double c42 = 0.0;
   const double c43 = 250.0/621.0;
   const double c44 = 125.0/594.0;
   const double c45 = 0.0;
   const double c46 = 512.0/1771.0;

   const double c51 = 2825.0/27648.0;  
   const double c52 = 0.0;  
   const double c53 = 18575.0/48384.0;   
   const double c54 = 13525.0/55296.0;
   const double c55 = 277.0/14336.0;
   const double c56 = 1.0/4.0;

   const double d1 = c41 - c51;
   const double d2 = c42 - c52;
   const double d3 = c43 - c53;
   const double d4 = c44 - c54;
   const double d5 = c45 - c55;
   const double d6 = c46 - c56;

   const size_t variables_number = count_variables_number();

   Vector<double> errors(dependent_variables_number, 0.0);

   double error = 0.0;

   Vector< Vector<double> > c(6);

   double hmin = 0.0;
   double h = (final_independent_variable - initial_independent_variable)*1.0e-3;

   Vector<double> variables(variables_number);

   Vector<double> final_solution(variables_number);

   // Initial values

   final_solution[0] = initial_independent_variable;

   for(size_t j = 0; j < dependent_variables_number; j++)
   {
      final_solution[1+j] = initial_dependent_variables[j];
   } 

   // Main loop

   while(final_solution[0] < final_independent_variable)
   {   
      // Set smallest allowable stepsize

      hmin = 32.0*epsilon*fabs(final_solution[0]);

      if(h < hmin)
      {
         if(display)
         {
            std::cout << "OpenNN Warning: OrdinaryDifferentialEquations class.\n" 
                      << "calculate_Runge_Kutta_Fehlberg_solution() method.\n" 
                      << "Step size is less than smallest allowable." << std::endl;
         }

         h = hmin;
      }

      if(final_solution[0] + h > final_independent_variable)
      {        
         h = final_independent_variable - final_solution[0];
      }

      // First coefficient

      variables[0] = final_solution[0]; 

      for(size_t j = 0; j < dependent_variables_number; j++)
      {
         variables[1+j] = final_solution[1+j];
      }

      c[0] = calculate_dependent_variables_dots(neural_network, variables)*h;

      // Second coefficient

      variables[0] = final_solution[0] + a2*h;

      for(size_t j = 0; j < dependent_variables_number; j++)
      {
         variables[1+j] = final_solution[1+j]+ b21*c[0][j];
      }

      c[1] = calculate_dependent_variables_dots(neural_network, variables)*h;

      // Third coefficient

      variables[0] = final_solution[0] + a3*h;

      for(size_t j = 0; j < dependent_variables_number; j++)
      {
         variables[1+j] =  final_solution[1+j] + b31*c[0][j] + b32*c[1][j]; 
      }

      c[2] = calculate_dependent_variables_dots(neural_network, variables)*h;

      // Fourth coefficient

      variables[0] = final_solution[0] + a4*h; 

      for(size_t j = 0; j < dependent_variables_number; j++)
      {
         variables[1+j] = final_solution[1+j] + b41*c[0][j] + b42*c[1][j] + b43*c[2][j];
      }

      c[3] = calculate_dependent_variables_dots(neural_network, variables)*h; 

      // Fifth coefficient

      variables[0] = final_solution[0] + a5*h;

      for(size_t j = 0; j < dependent_variables_number; j++)
      {
         variables[1+j] = final_solution[1+j] + b51*c[0][j] + b52*c[1][j] + b53*c[2][j] + b54*c[3][j];
      }
      
      c[4] = calculate_dependent_variables_dots(neural_network, variables)*h;
            
      // Sixth coefficient

      variables[0] = final_solution[0] + a6*h; 
            
      for(size_t j = 0; j < dependent_variables_number; j++)
      {
         variables[1+j] = final_solution[1+j] +  b61*c[0][j] + b62*c[1][j] + b63*c[2][j] + b64*c[3][j] + b65*c[4][j];
      }
      
      c[5] = calculate_dependent_variables_dots(neural_network, variables)*h;

      // Error estimate

      for(size_t j = 0; j < dependent_variables_number; j++)
      {
         errors[j] = fabs(d1*c[0][j] + d2*c[1][j] + d3*c[2][j] + d4*c[3][j] + d5*c[4][j] + d6*c[5][j]);         
      }

      error = errors.calculate_maximum();
   
      if(error <= tolerance)
      {
         final_solution[0] +=  h;

         for(size_t j = 0; j < dependent_variables_number; j++)
         {
            final_solution[1+j] += c51*c[0][j] + c52*c[1][j] + c53*c[2][j] + c54*c[3][j] + c55*c[4][j] + c56*c[5][j];   
         }

         if(error != 0.0)
         {
            h *= 0.9*pow(fabs(tolerance/error), 0.2);
         }
      }
      else
      {
         h *= 0.9*pow(fabs(tolerance/error), 0.25);
      }
   } // end while loop   


   return(final_solution);
}


// Matrix<double> calculate_solutions(const NeuralNetwork&) const method

Matrix<double> OrdinaryDifferentialEquations::calculate_solutions(const NeuralNetwork& neural_network) const
{
   switch(solution_method)
   {
      case RungeKutta:
      {
         return(calculate_Runge_Kutta_solution(neural_network));
      }            
      break;

      case RungeKuttaFehlberg:
      {
         return(calculate_Runge_Kutta_Fehlberg_solution(neural_network));
      }
      break;

      default:
      {
         std::ostringstream buffer;

         buffer << "OpenNN Error: OrdinaryDifferentialEquations class\n"
                << "Vector<double> calculate_solutions(const NeuralNetwork&) const method.\n"               
                << "Unknown solution method.\n";

         throw std::logic_error(buffer.str());
      }
      break;
   }
}


// Vector<double> calculate_final_solutions(const NeuralNetwork&) const method

Vector<double> OrdinaryDifferentialEquations::calculate_final_solutions(const NeuralNetwork& neural_network) const
{
   switch(solution_method)
   {
      case RungeKutta:
      {
         return(calculate_Runge_Kutta_final_solution(neural_network));
      }            
      break;

      case RungeKuttaFehlberg:
      {
         return(calculate_Runge_Kutta_Fehlberg_final_solution(neural_network));
      }
      break;

      default:
      {
         std::ostringstream buffer;

         buffer << "OpenNN Exception: ScalingLayer class\n"
                << "Vector<double> calculate_final_solutions(const NeuralNetwork&) const method.\n"               
                << "Unknown solution method.\n";

         throw std::logic_error(buffer.str());
      }
      break;
   }
}


// std::string to_string(void) const method


std::string OrdinaryDifferentialEquations::to_string(void) const
{
   std::ostringstream buffer; 

   buffer << "Mathematical model\n"
          << "Independent variables number: " << independent_variables_number << "\n" 
          << "Dependent variables number: " << dependent_variables_number << "\n"
          << "Initial independent variable: " << initial_independent_variable << "\n"
          << "Final independent variable: " << final_independent_variable << "\n"
          << "Initial dependent variables: " <<  initial_dependent_variables << "\n"
          << "Solution method: " << write_solution_method() << "\n"
          << "Points number: " << points_number << "\n"
          << "Tolerance" << tolerance << "\n"
          << "Initial size: " << initial_size << "\n"
          << "Warning size: " << warning_size << "\n"
          << "Error size: " << error_size << "\n"
          << "Display: " << display << std::endl;

   return(buffer.str());
}


// tinyxml2::XMLDocument* to_XML(void) const method


tinyxml2::XMLDocument* OrdinaryDifferentialEquations::to_XML(void) const
{    
   tinyxml2::XMLDocument* document = new tinyxml2::XMLDocument;

   std::ostringstream buffer;

   tinyxml2::XMLElement* ordinary_differential_equations_element = document->NewElement("OrdinaryDifferentialEquations");

   document->InsertFirstChild(ordinary_differential_equations_element);

   // Independent variables number 
   {
      tinyxml2::XMLElement* element = document->NewElement("IndependentVariablesNumber");
      ordinary_differential_equations_element->LinkEndChild(element);

      buffer.str("");
      buffer << independent_variables_number;

      tinyxml2::XMLText* text = document->NewText(buffer.str().c_str());
      element->LinkEndChild(text);
   }

   // Dependent variables number 
   {
      tinyxml2::XMLElement* element = document->NewElement("DependentVariablesNumber");
      ordinary_differential_equations_element->LinkEndChild(element);

      buffer.str("");
      buffer << dependent_variables_number;

      tinyxml2::XMLText* text = document->NewText(buffer.str().c_str());
      element->LinkEndChild(text);
   }

   // Initial independent variable
   {
      tinyxml2::XMLElement* element = document->NewElement("InitialIndependentVariable");
      ordinary_differential_equations_element->LinkEndChild(element);

      buffer.str("");
      buffer << initial_independent_variable;

      tinyxml2::XMLText* text = document->NewText(buffer.str().c_str());
      element->LinkEndChild(text);
   }

   // Final independent variable
   {
      tinyxml2::XMLElement* element = document->NewElement("FinalIndependentVariable");
      ordinary_differential_equations_element->LinkEndChild(element);

      buffer.str("");
      buffer << final_independent_variable;

      tinyxml2::XMLText* text = document->NewText(buffer.str().c_str());
      element->LinkEndChild(text);
   }

   // Initial dependent variables 
   {
      tinyxml2::XMLElement* element = document->NewElement("InitialDependentVariables");
      ordinary_differential_equations_element->LinkEndChild(element);

      buffer.str("");
      buffer << initial_dependent_variables;

      tinyxml2::XMLText* text = document->NewText(buffer.str().c_str());
      element->LinkEndChild(text);
   }

   // Solution method
   {
      tinyxml2::XMLElement* element = document->NewElement("SolutionMethod");
      ordinary_differential_equations_element->LinkEndChild(element);

      tinyxml2::XMLText* text = document->NewText(write_solution_method().c_str());
      element->LinkEndChild(text);
   }

   // Points number
   {
      tinyxml2::XMLElement* element = document->NewElement("PointsNumber");
      ordinary_differential_equations_element->LinkEndChild(element);

      buffer.str("");
      buffer << points_number;

      tinyxml2::XMLText* text = document->NewText(buffer.str().c_str());
      element->LinkEndChild(text);
   }

   // Tolerance
   {
      tinyxml2::XMLElement* element = document->NewElement("Tolerance");
      ordinary_differential_equations_element->LinkEndChild(element);

      buffer.str("");
      buffer << tolerance;

      tinyxml2::XMLText* text = document->NewText(buffer.str().c_str());
      element->LinkEndChild(text);
   }

   // Initial size
   {
      tinyxml2::XMLElement* element = document->NewElement("InitialSize");
      ordinary_differential_equations_element->LinkEndChild(element);

      buffer.str("");
      buffer << initial_size;

      tinyxml2::XMLText* text = document->NewText(buffer.str().c_str());
      element->LinkEndChild(text);
   }

   // Warning size
   {
      tinyxml2::XMLElement* element = document->NewElement("WarningSize");
      ordinary_differential_equations_element->LinkEndChild(element);

      buffer.str("");
      buffer << warning_size;

      tinyxml2::XMLText* text = document->NewText(buffer.str().c_str());
      element->LinkEndChild(text);
   }

   // Error size
   {
      tinyxml2::XMLElement* element = document->NewElement("ErrorSize");
      ordinary_differential_equations_element->LinkEndChild(element);

      buffer.str("");
      buffer << error_size;

      tinyxml2::XMLText* text = document->NewText(buffer.str().c_str());
      element->LinkEndChild(text);
   }

   // Display 
   {
      tinyxml2::XMLElement* element = document->NewElement("Display");
      ordinary_differential_equations_element->LinkEndChild(element);

      buffer.str("");
      buffer << display;

      tinyxml2::XMLText* text = document->NewText(buffer.str().c_str());
      element->LinkEndChild(text);
   }

   return(document);   
}


// void from_XML(const tinyxml2::XMLDocument&) method


void OrdinaryDifferentialEquations::from_XML(const tinyxml2::XMLDocument& document)
{
  // Dependent variables number
  {
     const tinyxml2::XMLElement* element = document.FirstChildElement("DependentVariablesNumber");

     if(element)
     {
        const char* text = element->GetText();

        if(text)
        {
           try
           {
              set_dependent_variables_number(atoi(text));
           }
           catch(const std::logic_error& e)
           {
              std::cout << e.what() << std::endl;
           }
        }
     }
  }

  // Initial independent variable
  {
     const tinyxml2::XMLElement* element = document.FirstChildElement("InitialIndependentVariable");

     if(element)
     {
        const char* text = element->GetText();

        if(text)
        {
           try
           {
              set_initial_independent_variable(atof(text));
           }
           catch(const std::logic_error& e)
           {
              std::cout << e.what() << std::endl;
           }
        }
     }
  }

  // Final independent variable
  {
     const tinyxml2::XMLElement* element = document.FirstChildElement("FinalIndependentVariable");

     if(element)
     {
        const char* text = element->GetText();

        if(text)
        {
           try
           {
              set_final_independent_variable(atof(text));
           }
           catch(const std::logic_error& e)
           {
              std::cout << e.what() << std::endl;
           }
        }
     }
  }

  // Initial dependent variables
  {
     const tinyxml2::XMLElement* element = document.FirstChildElement("InitialDependentVariables");

     if(element)
     {
        const char* text = element->GetText();

        if(text)
        {
           try
           {
              Vector<double> new_initial_dependent_variables;
              new_initial_dependent_variables.parse(text);

              set_initial_dependent_variables(new_initial_dependent_variables);
           }
           catch(const std::logic_error& e)
           {
              std::cout << e.what() << std::endl;
           }
        }
     }
  }

  // Solution method
  {
     const tinyxml2::XMLElement* element = document.FirstChildElement("SolutionMethod");

     if(element)
     {
        const char* text = element->GetText();

        if(text)
        {
           try
           {
              std::string new_solution_method(text);

              set_solution_method(new_solution_method);
           }
           catch(const std::logic_error& e)
           {
              std::cout << e.what() << std::endl;
           }
        }
     }
  }

  // Points number
  {
     const tinyxml2::XMLElement* element = document.FirstChildElement("PointsNumber");

     if(element)
     {
        const char* text = element->GetText();

        if(text)
        {
           try
           {
              set_points_number(atoi(text));
           }
           catch(const std::logic_error& e)
           {
              std::cout << e.what() << std::endl;
           }
        }
     }
  }

  // Tolerance
  {
     const tinyxml2::XMLElement* element = document.FirstChildElement("Tolerance");

     if(element)
     {
        const char* text = element->GetText();

        if(text)
        {
           try
           {
              set_tolerance(atof(text));
           }
           catch(const std::logic_error& e)
           {
              std::cout << e.what() << std::endl;
           }
        }
     }
  }

  // Initial size
  {
     const tinyxml2::XMLElement* element = document.FirstChildElement("InitialSize");

     if(element)
     {
        const char* text = element->GetText();

        if(text)
        {
           try
           {
              set_initial_size(atoi(text));
           }
           catch(const std::logic_error& e)
           {
              std::cout << e.what() << std::endl;
           }
        }
     }
  }

  // Warning size
  {
     const tinyxml2::XMLElement* element = document.FirstChildElement("WarningSize");

     if(element)
     {
        const char* text = element->GetText();

        if(text)
        {
           try
           {
              set_warning_size(atoi(text));
           }
           catch(const std::logic_error& e)
           {
              std::cout << e.what() << std::endl;
           }
        }
     }
  }

  // Error size
  {
     const tinyxml2::XMLElement* element = document.FirstChildElement("ErrorSize");

     if(element)
     {
        const char* text = element->GetText();

        if(text)
        {
           try
           {
              set_error_size(atoi(text));
           }
           catch(const std::logic_error& e)
           {
              std::cout << e.what() << std::endl;
           }
        }
     }
  }

  // Display
  {
     const tinyxml2::XMLElement* display_element = document.FirstChildElement("Display");

     if(display_element)
     {
        const char* display_text = display_element->GetText();

        if(display_text)
        {
           try
           {
              std::string display_string(display_text);

              set_display(display_string != "0");
           }
           catch(const std::logic_error& e)
           {
              std::cout << e.what() << std::endl;
           }
        }
     }
  }
}


// void save_data(const NeuralNetwork&, const std::string&) const method


void OrdinaryDifferentialEquations::save_data(const NeuralNetwork& neural_network, const std::string& file_name) const
{
   const Matrix<double> solution = calculate_Runge_Kutta_solution(neural_network);

   solution.save(file_name);
}

}


// OpenNN: Open Neural Networks Library.
// Copyright (c) 2005-2015 Roberto Lopez.
//
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or any later version.
//
// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.

// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA