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

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

//.........这里部分代码省略.........
    mError.resize(number_of_variables);
    mk1.resize(number_of_variables);
    mk2.resize(number_of_variables);
    mk3.resize(number_of_variables);
    mk4.resize(number_of_variables);
    mk5.resize(number_of_variables);
    mk6.resize(number_of_variables);
    myk2.resize(number_of_variables);
    myk3.resize(number_of_variables);
    myk4.resize(number_of_variables);
    myk5.resize(number_of_variables);
    myk6.resize(number_of_variables);

    double current_time = startTime;
    double time_step = maxTimeStep;
    bool got_to_end = false;
    bool accurate_enough = false;
    unsigned number_of_time_steps = 0;

    if (outputSolution)
    {   // Write out ICs
        rSolution.rGetTimes().push_back(current_time);
        rSolution.rGetSolutions().push_back(rYValues);
    }

    // should never get here if this bool has been set to true;
    assert(!mStoppingEventOccurred);
    while (!got_to_end)
    {
        //std::cout << "New timestep\n" << std::flush;
        while (!accurate_enough)
        {
            accurate_enough = true; // assume it is OK until we check and find otherwise

            // Function that calls the appropriate one-step solver
            CalculateNextYValue(pOdeSystem,
                                time_step,
                                current_time,
                                rYValues,
                                rWorkingMemory);

            // Find the maximum error in this vector
            double max_error = -DBL_MAX;
            for (unsigned i=0; i<number_of_variables; i++)
            {
                if (mError[i] > max_error)
                {
                    max_error = mError[i];
                }
            }

            if (max_error > tolerance)
            {// Reject the step-size and do it again.
                accurate_enough = false;
                //std::cout << "Approximation rejected\n" << std::flush;
            }
            else
            {
                // step forward the time since step has now been made
                current_time = current_time + time_step;
                //std::cout << "Approximation accepted with time step = "<< time_step << "\n" << std::flush;
                //std::cout << "max_error = " << max_error << " tolerance = " << tolerance << "\n" << std::flush;
                if (outputSolution)
                {   // Write out ICs
                    //std::cout << "In solver Time = " << current_time << " y = " << rWorkingMemory[0] << "\n" << std::flush;
                    rSolution.rGetTimes().push_back(current_time);
                    rSolution.rGetSolutions().push_back(rWorkingMemory);
                    number_of_time_steps++;
                }
            }

            // Set a new step size based on the accuracy here
            AdjustStepSize(time_step, max_error, tolerance, maxTimeStep, minTimeStep);
        }

        // For the next timestep check the step doesn't go past the end...
        if (current_time + time_step > endTime)
        {   // Allow a smaller timestep for the final step.
            time_step = endTime - current_time;
        }

        if ( pOdeSystem->CalculateStoppingEvent(current_time,
                                                rWorkingMemory) == true )
        {
            mStoppingTime = current_time;
            mStoppingEventOccurred = true;
        }

        if (mStoppingEventOccurred || current_time>=endTime)
        {
            got_to_end = true;
        }

        // Approximation accepted - put it in rYValues
        rYValues.assign(rWorkingMemory.begin(), rWorkingMemory.end());
        accurate_enough = false; // for next loop.
        //std::cout << "Finished Time Step\n" << std::flush;
    }
    rSolution.SetNumberOfTimeSteps(number_of_time_steps);
}