C++ OdeSolution类代码示例

    void TestMirams2010WntOdeSystemSetup()
    {
#ifdef CHASTE_CVODE
        double wnt_level = 0.5;
        boost::shared_ptr<AbstractCellMutationState> p_state(new WildTypeCellMutationState);
        Mirams2010WntOdeSystem wnt_system(wnt_level, p_state);
        // 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;
        //OdeSolution solutions2;

        std::vector<double> initial_conditions = wnt_system.GetInitialConditions();
        std::cout << "Timings for 100 hours\n";
        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;
        // 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);

        // Decent results
        TS_ASSERT_DELTA(solutions.rGetSolutions()[end][0], 67.5011, 1e-4);
        TS_ASSERT_DELTA(solutions.rGetSolutions()[end][1], 67.5011, 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 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 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);
    }

    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);
    }

    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 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;
        }
    }

  // 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);
    }

    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;
        }
    }

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

    void Simulate(const std::string& rOutputDirName,
                  const std::string& rModelName,
                  boost::shared_ptr<AbstractCardiacCellInterface> pCell)
    {
        double end_time = GetAttribute(pCell, "SuggestedCycleLength", 700.0); // ms
        if (pCell->GetSolver() || dynamic_cast<AbstractRushLarsenCardiacCell*>(pCell.get()))
        {
            double dt = GetAttribute(pCell, "SuggestedForwardEulerTimestep", 0.0);
            if (dt > 0.0)
            {
                pCell->SetTimestep(dt);
            }
        }
#ifdef CHASTE_CVODE
        AbstractCvodeSystem* p_cvode_cell = dynamic_cast<AbstractCvodeSystem*>(pCell.get());
        if (p_cvode_cell)
        {
            // Set a larger max internal time steps per sampling interval (CVODE's default is 500)
            p_cvode_cell->SetMaxSteps(1000);
            // Numerical or analytic J for CVODE?
            if (!mUseCvodeJacobian)
            {
                p_cvode_cell->ForceUseOfNumericalJacobian();
            }
        }
#endif
        double sampling_interval = 1.0; // ms; used as max dt for CVODE too
        Timer::Reset();
        OdeSolution solution = pCell->Compute(0.0, end_time, sampling_interval);
        std::stringstream message;
        message << "Model " << rModelName << " writing to " << rOutputDirName << " took";
        Timer::Print(message.str());

        const unsigned output_freq = 10; // Only output every N samples
        solution.WriteToFile(rOutputDirName, rModelName, "ms", output_freq, false);
        // Check an AP was produced
        std::vector<double> voltages = solution.GetVariableAtIndex(pCell->GetVoltageIndex());
        CellProperties props(voltages, solution.rGetTimes());
        props.GetLastActionPotentialDuration(90.0); // Don't catch the exception here if it's thrown
        // Compare against saved results
        CheckResults(rModelName, voltages, solution.rGetTimes(), output_freq);
    }

    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);
    }

    void TestBackwardEulerNonlinearEquation()
    {
        Ode4 ode_system;

        double h_value = 0.01;

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

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

        solutions = backward_euler_solver.Solve(&ode_system, state_variables, 0.0, 2.0, h_value, h_value);
        unsigned 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);
    }

    void TestSimpleSystemWithOrderSwapped()
    {
        // The solution should be the same, but we'll have to construct a new CombinedOdeSystemInformation
        // object, because the order of subsystems has changed.

        SimpleOde1 ode_for_y; // dy/dt = x
        SimpleOde2 ode_for_x; // dx/dt = -y

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

        // 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;

        combined_ode_system.Configure(variable_parameter_map, &ode_for_y, &ode_for_x);

        // ...and vice versa (we can re-use the map in this case)
        combined_ode_system.Configure(variable_parameter_map, &ode_for_x, &ode_for_y);

        // Test solving the combined system.
        // This is dy/dt = x, dx/dt = -y, y(0) = 0, x(0) = 1.
        // The analytic solution is y = sin(t), x = cos(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()[1], sin(2.0), global_error);
        TS_ASSERT_DELTA(solutions.rGetSolutions().back()[0], cos(2.0), global_error);
    }