BaseCustomSolver.m

Go to the documentation of this file.

namespace solvers{


/* (Autoinserted by mtoc++)
 * This source code has been filtered by the mtoc++ executable,
 * which generates code that can be processed by the doxygen documentation tool.
 *
 * On the other hand, it can neither be interpreted by MATLAB, nor can it be compiled with a C++ compiler.
 * Except for the comments, the function bodies of your M-file functions are untouched.
 * Consequently, the FILTER_SOURCE_FILES doxygen switch (default in our Doxyfile.template) will produce
 * attached source files that are highly readable by humans.
 *
 * Additionally, links in the doxygen generated documentation to the source code of functions and class members refer to
 * the correct locations in the source code browser.
 * However, the line numbers most likely do not correspond to the line numbers in the original MATLAB source files.
 */

class BaseCustomSolver
  :public solvers.BaseSolver {
    public:

        BaseCustomSolver() {
            this = this@solvers.BaseSolver;
            this.Name= " Base Custom Solver ";
        }


        function [rowvect , matrixx ] = solve(function_handle odefun,rowvec t,colvec x0) {

            /*  Get computation times
             *[times, outputtimes] = this.getCompTimes(t); */
            [times, outidx] = this.getCompTimes(t);

            /*  Fire PreSolve event */
            ed = solvers.SolverEventData(times);
            notify(this," PreSolve ",ed);

            x = this.customSolve(odefun, times, x0, outidx);

            /*  Fire PostSolve event */
            ed.States= x;
            notify(this," PostSolve ",ed);

            /*  Extract wanted values */
            t = times(outidx);
            /*  New: x is supposed to be returned in the size of effectively desired timesteps
             * (memory issues for large, long simulations)
             *x = x(:,outputtimes); */
        }
    private:

        function [times , rowvec<integer>outidx ] = getCompTimes(rowvec<double> t) {

            /*  Validity checks */
            if any(abs(sort(t) - t) > 100*eps)
                error(" The time vector must be strictly monotoneously increasing. ");
            end

            /*  Default values
             * outputtimes = true(1,length(t)); */
            outidx = 1:length(t);
            times = t;
            if ~isempty(this.MaxStep)
                if numel(t) == 2
                    times = t(1):this.MaxStep:t(2);
                    /*  outputtimes = true(1,length(times)); */
                    outidx = [1,length(times)];
                else
                    len = length(t);
                    n = max(round((t(2:len)-t(1:len-1))/this.MaxStep),1);  /*  number of steps in interval [t(i-1),t(i)] after refinemant. If t(i-1)-t(i)>= 1.5*MaxStep, interval will be refined */

                    m = ones(1,len);      /*  m(i) is the index of the originally i-th entry in t in times  */

                    m(2:len) = cumsum(n) + m(2:len);/*  - (0:len-2); */

                    times = zeros(1,m(len));
                    /*  outputtimes = times; */
                    for i = 1:len-1
                        times(m(i):m(i+1))=  t(i):(t(i+1)-t(i))/(n(i)):t(i+1); /* linspace(t(i),t(i+1),n(i));
                         * outputtimes(m(i)) = 1; */

                        outidx(i) = m(i);
                    end
                    /* outputtimes(end) = 1; */
                    outidx(end) = length(times);
                end
                /* outputtimes = logical(outputtimes); */

                /*                      % Find refinement indices
                 *                     idx = fliplr(find(t(2:end)-t(1:end-1)-this.MaxStep>100*eps));
                 *                     % If any "gaps" are present, fill up with MaxStep's
                 *
                 *                     if ~isempty(idx)
                 *                         for i=idx
                 *                             add = times(i):this.MaxStep:times(i+1);
                 *                             times = [times(1:i-1) add times(i+2:end)];
                 *                             outputtimes = [outputtimes(1:i) false(1,length(add)-2) outputtimes(i+1:end)];
                 *                         end
                 *                         tout = times(outputtimes);
                 *                         if numel(tout) ~= numel(t) || any(abs(tout - t) > 100*eps)
                 *                             error('Unexpected error: Computed time vector differs from the desired one.');
                 *                             t%#ok
                 *                         end */
            end
        }
    protected: /* ( Abstract ) */

        virtual function colvec<double>x = customSolve(odefun,t,x0,outputtimes) = 0;
    public: /* ( Abstract ) */

        EVENT PreSolve;
        EVENT PostSolve;
};
}