C++ OdeSolution::rGetSolutions方法代码示例

    void TestGarysWntOdeSystemApc2Hit()
    {
#ifdef CHASTE_CVODE
        double wnt_level = 0.5;
        boost::shared_ptr<AbstractCellMutationState> p_apc2(new ApcTwoHitCellMutationState);
        Mirams2010WntOdeSystem wnt_system(wnt_level, p_apc2);

        // Solve system using CVODE solver
        // Matlab's strictest bit uses 0.01 below and relaxes it on flatter bits
        double h_value = 0.01;
        CvodeAdaptor cvode_solver;
        OdeSolution solutions;
        std::vector<double> initial_conditions = wnt_system.GetInitialConditions();

        Timer::Reset();
        solutions = cvode_solver.Solve(&wnt_system, initial_conditions, 0.0, 100.0, h_value, h_value);
        Timer::Print("1. Cvode");

        // Test solutions are OK for a small time increase
        int end = solutions.rGetSolutions().size() - 1;

        // Test the simulation is ending at the right time (going into S phase at 7.8 hours)
        TS_ASSERT_DELTA(solutions.rGetTimes()[end], 100, 1e-2);

        // Check results are correct
        TS_ASSERT_DELTA(solutions.rGetSolutions()[end][0], 433.114, 2e-3); // Tolerances relaxed for
        TS_ASSERT_DELTA(solutions.rGetSolutions()[end][1], 433.114, 2e-3); // different CVODE versions.
        TS_ASSERT_DELTA(solutions.rGetSolutions()[end][2], wnt_level, 1e-4);
#else
        std::cout << "CVODE is not enabled. " << std::endl;
        std::cout << "If required please install and alter your hostconfig settings to switch on chaste support." << std::endl;
#endif //CHASTE_CVODE
    }

    void TestBackwardEulerSystemOf3EquationsWithEvents()
    {
        OdeThirdOrderWithEvents ode_system_with_events;

        double h_value = 0.01;

        // Euler solver solution worked out
        BackwardEulerIvpOdeSolver backward_euler_solver(ode_system_with_events.GetNumberOfStateVariables());
        OdeSolution solutions;

        std::vector<double> state_variables = ode_system_with_events.GetInitialConditions();
        solutions = backward_euler_solver.Solve(&ode_system_with_events, state_variables, 0.0, 2.0, h_value, h_value);
        unsigned last = solutions.GetNumberOfTimeSteps();

        // Final time should be pi/6 (?)
        TS_ASSERT_DELTA( solutions.rGetTimes()[last], 0.5236, 0.01);

        // Penultimate y0 should be greater than -0.5
        TS_ASSERT_LESS_THAN(-0.5,solutions.rGetSolutions()[last-1][0]);

        // Final y0 should be less than -0.5
        TS_ASSERT_LESS_THAN( solutions.rGetSolutions()[last][0], -0.5);

        // Solver should correctly state the stopping event occurred
        TS_ASSERT_EQUALS(backward_euler_solver.StoppingEventOccurred(), true);
    }

    void TestRKFehlbergSystemOf3Equations() throw(Exception)
    {
        OdeThirdOrder ode_system;

        double h_value = 0.1;

        // Euler solver solution worked out
        RungeKuttaFehlbergIvpOdeSolver rkf_solver;

        OdeSolution solutions;

        std::vector<double> state_variables = ode_system.GetInitialConditions();
        solutions = rkf_solver.Solve(&ode_system, state_variables, 0.0, 2.0, 0.25, 1e-5);
        unsigned last = solutions.GetNumberOfTimeSteps();
        double numerical_solution[3];
        numerical_solution[0] = solutions.rGetSolutions()[last][0];
        numerical_solution[1] = solutions.rGetSolutions()[last][1];
        numerical_solution[2] = solutions.rGetSolutions()[last][2];

        // The tests
        double analytical_solution[3];
        analytical_solution[0] = -sin(2.0);
        analytical_solution[1] = sin(2.0)+cos(2.0);
        analytical_solution[2] = 2*sin(2.0);
        double global_error_rkf = 0.5*2*(exp(2.0)-1)*h_value;
        TS_ASSERT_DELTA(numerical_solution[0],analytical_solution[0],global_error_rkf);
        TS_ASSERT_DELTA(numerical_solution[1],analytical_solution[1],global_error_rkf);
        TS_ASSERT_DELTA(numerical_solution[2],analytical_solution[2],global_error_rkf);
    }

  // Test a given solver on an ODE which has a stopping event defined
  void MyTestSolverOnOdesWithEvents(AbstractIvpOdeSolver& rSolver)
  {
        // ODE which has solution y0 = cos(t), and stopping event y0<0,
        // ie should stop when t = pi/2;
        OdeSecondOrderWithEvents ode_with_events;

        OdeSolution solutions;
        std::vector<double> state_variables =
            ode_with_events.GetInitialConditions();
        solutions = rSolver.Solve(&ode_with_events, state_variables, 0.0, 2.0,
            0.001, 0.001);

        unsigned num_timesteps = solutions.GetNumberOfTimeSteps();

        // Final time should be around pi/2
        TS_ASSERT_DELTA( solutions.rGetTimes()[num_timesteps], M_PI_2, 0.01);

        // Penultimate y0 should be greater than zero
        TS_ASSERT_LESS_THAN( 0, solutions.rGetSolutions()[num_timesteps-1][0]);

        // Final y0 should be less than zero
        TS_ASSERT_LESS_THAN( solutions.rGetSolutions()[num_timesteps][0], 0);

        // Solver should correctly state the stopping event occurred
        TS_ASSERT_EQUALS(rSolver.StoppingEventOccurred(), true);

        // This is to cover the exception when a stopping event occurs before the first timestep.
        TS_ASSERT_THROWS_ANYTHING(rSolver.Solve(&ode_with_events, state_variables, 2.0, 3.0, 0.001));

        ///////////////////////////////////////////////
        // Repeat with sampling time larger than dt
        ///////////////////////////////////////////////

        state_variables = ode_with_events.GetInitialConditions();
        solutions = rSolver.Solve(&ode_with_events, state_variables, 0.0, 2.0,
            0.001, 0.01);

        num_timesteps = solutions.GetNumberOfTimeSteps();

        // Final time should be around pi/2
        TS_ASSERT_DELTA( solutions.rGetTimes()[num_timesteps], M_PI_2, 0.01);

        // Penultimate y0 should be greater than zero
        TS_ASSERT_LESS_THAN( 0, solutions.rGetSolutions()[num_timesteps-1][0]);

        // Final y0 should be less than zero
        TS_ASSERT_LESS_THAN( solutions.rGetSolutions()[num_timesteps][0], 0);

        // Solver should correctly state the stopping event occurred
        TS_ASSERT_EQUALS(rSolver.StoppingEventOccurred(), true);

        // Cover the check event isn't initially true exception
        std::vector<double> bad_init_cond;
        bad_init_cond.push_back(-1); //y0 < 0 so stopping event true
        bad_init_cond.push_back(0.0);
        TS_ASSERT_THROWS_ANYTHING(rSolver.Solve(&ode_with_events, bad_init_cond, 0.0, 2.0, 0.001, 0.01));
    }

    void TestWithThreeVariablesTwoSystems()
    {
        SimpleOde3 ode_for_x; // dx/dt = x -y +z
        SimpleOde6 ode_for_yz; // dy/dt = y-z  and dz/dt = 2y-z

        std::vector<AbstractOdeSystem*> ode_systems;
        ode_systems.push_back(&ode_for_x);
        ode_systems.push_back(&ode_for_yz);

        // Create combined ODE system
        CombinedOdeSystem combined_ode_system(ode_systems);

        // Tell the combined ODE system which state variables in the first ODE system
        // correspond to which parameters in the second ODE system...
        std::map<unsigned, unsigned> variable_parameter_map;
        variable_parameter_map[0] = 0; //y in the yz-ODE appears in the x-ODE
        variable_parameter_map[1] = 1; //z in the yz-ODE appears in the x-ODE
        combined_ode_system.Configure(variable_parameter_map, &ode_for_yz, &ode_for_x);

        // Test number of state variables
        unsigned num_variables = combined_ode_system.GetNumberOfStateVariables();
        TS_ASSERT_EQUALS(num_variables, 3u);

        // Combined system has no parameters
        TS_ASSERT_EQUALS(combined_ode_system.GetNumberOfParameters(), 0u);
        TS_ASSERT_EQUALS(combined_ode_system.rGetParameterNames().size(), 0u);

        // Test initial conditions
        std::vector<double> initial_conditions = combined_ode_system.GetInitialConditions();
        TS_ASSERT_DELTA(initial_conditions[0], 0.0, 1e-12);
        TS_ASSERT_DELTA(initial_conditions[1], 1.0, 1e-12);
        TS_ASSERT_DELTA(initial_conditions[2], 0.0, 1e-12);
        // Test variable names & units
        const std::vector<std::string>& r_names = combined_ode_system.rGetStateVariableNames();
        TS_ASSERT_EQUALS(r_names[0], ode_for_x.rGetStateVariableNames()[0]);
        TS_ASSERT_EQUALS(r_names[1], ode_for_yz.rGetStateVariableNames()[0]);
        TS_ASSERT_EQUALS(r_names[2], ode_for_yz.rGetStateVariableNames()[1]);

        // x'=x-y+z, y'=y-z, z'=2y-z
        // starting at (x,y,z)=(0,1,0)
        // Analytic solution is x=-sin(t), y=sin(t)+cos(t), z=2sin(t)
        EulerIvpOdeSolver solver;
        OdeSolution solutions;
        double h = 0.01;
        std::vector<double> inits = combined_ode_system.GetInitialConditions();
        solutions = solver.Solve(&combined_ode_system, inits, 0.0, 2.0, h, h);
        double global_error = 0.5*(exp(2.0)-1)*h;
        TS_ASSERT_DELTA(solutions.rGetSolutions().back()[0], -sin(2.0), global_error);
        TS_ASSERT_DELTA(solutions.rGetSolutions().back()[1], sin(2.0)+cos(2.0), global_error);
        TS_ASSERT_DELTA(solutions.rGetSolutions().back()[2], 2.0*sin(2.0), global_error);
    }

OdeSolution AbstractRushLarsenCardiacCell::Compute(double tStart, double tEnd, double tSamp)
{
    // In this method, we iterate over timesteps, doing the following for each:
    //   - update V using a forward Euler step
    //   - do as in ComputeExceptVoltage(t) to update the remaining state variables
    //     using Rush Larsen method or forward Euler as appropriate

    // Check length of time interval
    if (tSamp < mDt)
    {
        tSamp = mDt;
    }
    const unsigned n_steps = (unsigned) floor((tEnd - tStart)/tSamp + 0.5);
    assert(fabs(tStart+n_steps*tSamp - tEnd) < 1e-12);
    const unsigned n_small_steps = (unsigned) floor(tSamp/mDt+0.5);
    assert(fabs(mDt*n_small_steps - tSamp) < 1e-12);

    // Initialise solution store
    OdeSolution solutions;
    solutions.SetNumberOfTimeSteps(n_steps);
    solutions.rGetSolutions().push_back(rGetStateVariables());
    solutions.rGetTimes().push_back(tStart);
    solutions.SetOdeSystemInformation(this->mpSystemInfo);

    std::vector<double> dy(mNumberOfStateVariables, 0);
    std::vector<double> alpha(mNumberOfStateVariables, 0);
    std::vector<double> beta(mNumberOfStateVariables, 0);

    // Loop over time
    for (unsigned i=0; i<n_steps; i++)
    {
        double curr_time = tStart;
        for (unsigned j=0; j<n_small_steps; j++)
        {
            curr_time = tStart + i*tSamp + j*mDt;
            EvaluateEquations(curr_time, dy, alpha, beta);
            UpdateTransmembranePotential(dy);
            ComputeOneStepExceptVoltage(dy, alpha, beta);
            VerifyStateVariables();
        }

        // Update solutions
        solutions.rGetSolutions().push_back(rGetStateVariables());
        solutions.rGetTimes().push_back(curr_time+mDt);
    }

    return solutions;
}

    void TestBackwardEulerVanDerPolOde()
    {
        VanDerPolOde ode_system;

        double h_value = 0.01;
        double end_time = 100.0;

        // Euler solver solution worked out
        BackwardEulerIvpOdeSolver backward_euler_solver(ode_system.GetNumberOfStateVariables());
        backward_euler_solver.ForceUseOfNumericalJacobian(); // coverage
        OdeSolution solutions;

        std::vector<double> state_variables = ode_system.GetInitialConditions();

        solutions = backward_euler_solver.Solve(&ode_system, state_variables, 0.0, end_time, h_value, 5*h_value);
        unsigned last = solutions.GetNumberOfTimeSteps();

//        OutputFileHandler handler("");
//        out_stream rabbit_file=handler.OpenOutputFile("foxrabbit.dat");
//
//        for (unsigned i=0; i<last; i++)
//        {
//            (*rabbit_file) << solutions.rGetSolutions()[i][0] << "\t" << solutions.rGetSolutions()[i][1] << "\n" << std::flush;
//        }
//        rabbit_file->close();

        // assert that we are within a [-2,2] in x and [-2,2] in y (on limit cycle)
        TS_ASSERT_DELTA(solutions.rGetSolutions()[last][0], 0, 2);
        TS_ASSERT_DELTA(solutions.rGetSolutions()[last][1], 0, 2);
    }

    /**
     * Test two ODE solvers with this ODE system (correct values calculated using the Matlab solver ode15s).
     *
     */
    void TestAlarcon2004Solver()
    {
        // Set up
        double oxygen_concentration = 1.0;
        Alarcon2004OxygenBasedCellCycleOdeSystem alarcon_system(oxygen_concentration, false);

        // Create ODE solvers
        RungeKutta4IvpOdeSolver rk4_solver;
        RungeKuttaFehlbergIvpOdeSolver rkf_solver;
        BackwardEulerIvpOdeSolver back_solver(6);

        // Set up for solver
        OdeSolution solutions;
        std::vector<double> initial_conditions = alarcon_system.GetInitialConditions();
        double h_value = 1e-4; // maximum tolerance

        // Solve the ODE system using a Runge Kutta fourth order solver
        Timer::Reset();
        solutions = rk4_solver.Solve(&alarcon_system, initial_conditions, 0.0, 10.0, h_value, h_value);
        Timer::Print("1. Runge-Kutta");

        // Reset maximum tolerance for Runge Kutta Fehlber solver
        h_value = 1e-1;

        // Solve the ODE system using a Runge Kutta Fehlber solver
        initial_conditions = alarcon_system.GetInitialConditions();
        Timer::Reset();
        solutions = rkf_solver.Solve(&alarcon_system, initial_conditions, 0.0, 10.0, h_value, 1e-4);
        Timer::Print("1. Runge-Kutta-Fehlberg");

        // Test that solutions are accurate for a small time increase
        int end = solutions.rGetSolutions().size() - 1;

        // Test that the solver stops at the right time
        TS_ASSERT_DELTA(solutions.rGetTimes()[end], 9.286356375, 1e-2);

        // Test solution - note the high tolerances
        TS_ASSERT_DELTA(solutions.rGetSolutions()[end][0], 0.004000000000000, 1e-3);
        TS_ASSERT_DELTA(solutions.rGetSolutions()[end][1], 0.379221366479055, 1e-3);
        TS_ASSERT_DELTA(solutions.rGetSolutions()[end][2], 0.190488726735972, 1e-3);
        TS_ASSERT_DELTA(solutions.rGetSolutions()[end][3], 9.962110289977730, 1e-3);
        TS_ASSERT_DELTA(solutions.rGetSolutions()[end][4], 0.096476600742599, 1e-3);
        TS_ASSERT_DELTA(solutions.rGetSolutions()[end][5], 1.000000000000000, 1e-3);
    }

    void TestSolvingOdes() throw(Exception)
    {
        /* First, create an instance of the ODE class to be solved. */
        MyOde my_ode;
        /* Next, create a solver. */
        EulerIvpOdeSolver euler_solver;
        /* We will need to provide an initial condition, which needs to
         * be a {{{std::vector}}}.*/
        std::vector<double> initial_condition;
        initial_condition.push_back(1.0);
        /* Then, just call `Solve`, passing in a pointer to the ODE, the
         * initial condition, the start time, end time, the solving timestep,
         * and sampling timestep (how often we want the solution stored in the returned `OdeSolution` object).
         * Here we solve from 0 to 1, with a timestep of 0.01 but a ''sampling
         * timestep'' of 0.1. The return value is an object of type {{{OdeSolution}}}
         * (which is basically just a list of times and solutions).
         */
        OdeSolution solutions = euler_solver.Solve(&my_ode, initial_condition, 0, 1, 0.01, 0.1);
        /* Let's look at the results, which can be obtained from the {{{OdeSolution}}}
         * object using the methods {{{rGetTimes()}}} and {{{rGetSolutions()}}}, which
         * return a {{{std::vector}}} and a {{{std::vector}}} of {{{std::vector}}}s
         * respectively. */
        for (unsigned i=0; i<solutions.rGetTimes().size(); i++)
        {
            /* The {{{[0]}}} here is because we are getting the zeroth component of y (a 1-dimensional vector). */
            std::cout << solutions.rGetTimes()[i] << " " << solutions.rGetSolutions()[i][0] << "\n";
        }

        /* Alternatively, we can print the solution directly to a file, using the {{{WriteToFile}}}
         * method on the {{{OdeSolution}}} class. */
        solutions.WriteToFile("SolvingOdesTutorial", "my_ode_solution", "sec");
        /* Two files are written
         * * {{{my_ode_solution.dat}}} contains the results (a header line, then one column for time and one column per variable)
         * * {{{my_ode_solution.info}}} contains information for reading the data back, a line about the ODE solver ("{{{ODE SOLVER: EulerIvpOdeSolver}}}") and provenance information.
         */

        /* We can see from the printed out results that y goes above 2.5 somewhere just
         * before 0.6. To solve only up until y=2.5, we can solve the ODE that has the
         * stopping event defined, using the same solver as before. */
        MyOdeWithStoppingEvent my_ode_stopping;

        /* '''Note:''' ''when a {{{std::vector}}} is passed in as an initial condition
         * to a {{{Solve}}} call, it gets updated as the solve takes place''. Therefore, if
         * we want to use the same initial condition again, we have to reset it back to 1.0. */
        initial_condition[0] = 1.0;
        solutions = euler_solver.Solve(&my_ode_stopping, initial_condition, 0, 1, 0.01, 0.1);
        /* We can check with the solver that it stopped because of the stopping event, rather than because
         * it reached to end time. */
        TS_ASSERT(euler_solver.StoppingEventOccurred());
        /* Finally, let's print the time of the stopping event (to the nearest dt or so). */
        std::cout << "Stopping event occurred at t="<<solutions.rGetTimes().back()<<"\n";
    }

    void TestArchivingRkfSolver() throw(Exception)
    {
        OutputFileHandler handler("archive",false);
        std::string archive_filename;
        archive_filename = handler.GetOutputDirectoryFullPath() + "rkf_solver.arch";

        Ode5Jacobian ode_system;
        OdeSolution solutions;
        double h_value = 0.1;
        double end_time = 1.0;

        // Create and archive simulation time
        {
            std::ofstream ofs(archive_filename.c_str());
            boost::archive::text_oarchive output_arch(ofs);

            // Set up a solver
            AbstractIvpOdeSolver* const p_rkf_ode_solver = new RungeKuttaFehlbergIvpOdeSolver;


            // Should always archive a pointer
            output_arch << p_rkf_ode_solver;

            // Change stimulus a bit
            delete p_rkf_ode_solver;
        }

        // Restore
        {
            std::ifstream ifs(archive_filename.c_str(), std::ios::binary);
            boost::archive::text_iarchive input_arch(ifs);

            // Create a pointer
            AbstractIvpOdeSolver* p_rkf;
            input_arch >> p_rkf;

            std::vector<double> state_variables = ode_system.GetInitialConditions();

            solutions = p_rkf->Solve(&ode_system, state_variables, 0.0, end_time, h_value, 1e-5);
            unsigned last = solutions.GetNumberOfTimeSteps();

            double numerical_solution;
            numerical_solution = solutions.rGetSolutions()[last][0];

            // The tests
            double analytical_solution = 1.0/(1.0+4.0*exp(-100.0*end_time));

            TS_ASSERT_DELTA(numerical_solution,analytical_solution,1.0e-3);

            delete p_rkf;
        }
    }

    void TestGarysWntOdeSystemBetaCatenin1Hit() throw(Exception)
    {
#ifdef CHASTE_CVODE
        double wnt_level = 0.5;
        boost::shared_ptr<AbstractCellMutationState> p_bcat1(new BetaCateninOneHitCellMutationState);
        Mirams2010WntOdeSystem wnt_system(wnt_level, p_bcat1);

        // Solve system using CVODE solver
        // Matlab's strictest bit uses 0.01 below and relaxes it on flatter bits
        double h_value = 0.1;
        CvodeAdaptor cvode_solver;
        OdeSolution solutions;
        std::vector<double> initial_conditions = wnt_system.GetInitialConditions();

        double start_time, end_time, elapsed_time = 0.0;
        start_time = (double) std::clock();
        solutions = cvode_solver.Solve(&wnt_system, initial_conditions, 0.0, 100.0, h_value, h_value);
        end_time = (double) std::clock();
        elapsed_time = (end_time - start_time)/(CLOCKS_PER_SEC);
        std::cout << "1. Cvode Elapsed time = " << elapsed_time << " secs for 100 hours\n";

        // Test solutions are OK for a small time increase
        int end = solutions.rGetSolutions().size() - 1;

        // Tests the simulation is ending at the right time (going into S phase at 7.8 hours)
        TS_ASSERT_DELTA(solutions.rGetTimes()[end], 100, 1e-2);

        // Check results are correct
        TS_ASSERT_DELTA(solutions.rGetSolutions()[end][0], 67.5011, 1e-4);
        TS_ASSERT_DELTA(solutions.rGetSolutions()[end][1], 824.0259, 1e-4);
        TS_ASSERT_DELTA(solutions.rGetSolutions()[end][2], wnt_level, 1e-4);
#else
        std::cout << "CVODE is not enabled. " << std::endl;
        std::cout << "If required please install and alter your hostconfig settings to switch on chaste support." << std::endl;
#endif //CHASTE_CVODE
    }

    void TestArchivingSolver() throw(Exception)
    {
        OutputFileHandler handler("archive", false);
        ArchiveLocationInfo::SetArchiveDirectory(handler.FindFile(""));
        std::string archive_filename = ArchiveLocationInfo::GetProcessUniqueFilePath("backward_euler_solver.arch");

        VanDerPolOde ode_system;

        double h_value = 0.01;
        double end_time = 100.0;

        // Create and archive simulation time
        {
            std::ofstream ofs(archive_filename.c_str());
            boost::archive::text_oarchive output_arch(ofs);

            // Set up a solver
            AbstractIvpOdeSolver* const p_backward_euler_solver = new BackwardEulerIvpOdeSolver(ode_system.GetNumberOfStateVariables());

            // Should always archive a pointer
            output_arch << p_backward_euler_solver;

            // Change stimulus a bit
            delete p_backward_euler_solver;
        }

        // Restore
        {
            std::ifstream ifs(archive_filename.c_str(), std::ios::binary);
            boost::archive::text_iarchive input_arch(ifs);

            // Create a pointer
            AbstractIvpOdeSolver* p_backward_euler;
            input_arch >> p_backward_euler;
            OdeSolution solutions;

            std::vector<double> state_variables = ode_system.GetInitialConditions();

            solutions = p_backward_euler->Solve(&ode_system, state_variables, 0.0, end_time, h_value, 5*h_value);
            unsigned last = solutions.GetNumberOfTimeSteps();

            // assert that we are within a [-2,2] in x and [-2,2] in y (on limit cycle)
            TS_ASSERT_DELTA(solutions.rGetSolutions()[last][0], 0, 2);
            TS_ASSERT_DELTA(solutions.rGetSolutions()[last][1], 0, 2);

            delete p_backward_euler;
        }
    }

    void TestRKFehlbergWithExampleFromBook() throw(Exception)
    {
        /*
         * Book is "Numerical Analysis 6th Edition by R.L. Burden and J. D. Faires
         * This example is on P291 Table 5.9
         */
        RkfTestOde ode;

        double max_step_size = 0.25;
        double start_time = 0.0;
        double end_time = 2.0;
        RungeKuttaFehlbergIvpOdeSolver rkf_solver;

        OdeSolution solutions;

        std::vector<double> state_variables = ode.GetInitialConditions();
        double tolerance = 1e-5;
        solutions = rkf_solver.Solve(&ode, state_variables, start_time, end_time, max_step_size, tolerance);

        // Times (from MatLab Code) to check timstepping is being adapted properly
        TS_ASSERT_DELTA(solutions.rGetTimes()[0], 0, 1e-7);
        TS_ASSERT_DELTA(solutions.rGetTimes()[1], 2.500000000000000e-01, 1e-7);
        TS_ASSERT_DELTA(solutions.rGetTimes()[2], 4.868046415733731e-01, 1e-7);
        TS_ASSERT_DELTA(solutions.rGetTimes()[3], 7.298511818781566e-01, 1e-7);
        TS_ASSERT_DELTA(solutions.rGetTimes()[4], 9.798511818781566e-01, 1e-7);
        TS_ASSERT_DELTA(solutions.rGetTimes()[5], 1.229851181878157e+00, 1e-7);
        TS_ASSERT_DELTA(solutions.rGetTimes()[6], 1.479851181878157e+00, 1e-7);
        TS_ASSERT_DELTA(solutions.rGetTimes()[7], 1.729851181878157e+00, 1e-7);
        TS_ASSERT_DELTA(solutions.rGetTimes()[8], 1.979851181878157e+00, 1e-7);
        TS_ASSERT_DELTA(solutions.rGetTimes()[9], 2.000000000000000e+00, 1e-7);

        TS_ASSERT_EQUALS(solutions.GetNumberOfTimeSteps(), 9u);

        // y values (from analytic result)
        for (unsigned i=0; i<solutions.GetNumberOfTimeSteps(); i++)
        {
            double time = solutions.rGetTimes()[i];
            double y = (time+1.0)*(time+1.0) - 0.5*exp(time);

            // Tolerance set to 1e-5, so 2e-5 to pass here
            TS_ASSERT_DELTA(solutions.rGetSolutions()[i][0], y, 2e-5);
        }
    }

    /*
     * === Solving n-dimensional ODEs ===
     *
     * Finally, here's a simple test showing how to solve a 2d ODE using the first method.
     * All that is different is the initial condition has be a vector of length 2, and returned
     * solution is of length 2 at every timestep.
     */
    void TestWith2dOde()
    {
        My2dOde my_2d_ode;
        EulerIvpOdeSolver euler_solver;

        /* Define the initial condition for each state variable. */
        std::vector<double> initial_condition;
        initial_condition.push_back(1.0);
        initial_condition.push_back(0.0);

        /* Solve, and print the solution as [time, y1, y2]. */
        OdeSolution solutions = euler_solver.Solve(&my_2d_ode, initial_condition, 0, 1, 0.01, 0.1);
        for (unsigned i=0; i<solutions.rGetTimes().size(); i++)
        {
            std::cout << solutions.rGetTimes()[i] << " "
                      << solutions.rGetSolutions()[i][0] << " "
                      << solutions.rGetSolutions()[i][1] << "\n";
        }
    }

    void TestRKFehlbergNonlinearEquation() throw(Exception)
    {
        Ode4 ode_system;

        double h_value = 0.1;

        // Euler solver solution worked out
        RungeKuttaFehlbergIvpOdeSolver rkf_solver;
        OdeSolution solutions;

        std::vector<double> state_variables = ode_system.GetInitialConditions();

        solutions = rkf_solver.Solve(&ode_system, state_variables, 0.0, 2.0, h_value, 1e-5);
        int last = solutions.GetNumberOfTimeSteps();

        double numerical_solution;
        numerical_solution = solutions.rGetSolutions()[last][0];

        // The tests
        double analytical_solution = 1.0/(1.0+exp(-12.5));

        TS_ASSERT_DELTA(numerical_solution,analytical_solution,1.0e-4);
    }