MLode15i.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 MLode15i
  :public solvers.MLWrapper {
    public: /* ( setObservable ) */

        double RelTol = 1e-3;
        double AbsTol = 1e-6;
    private:

        fED;

            
    public:

        MLode15i() {
            this = this@solvers.MLWrapper(@ode15i);
            this.Name= " MatLab ode15i implicit solver wrapper ";
            this.registerProps(" RelTol "," AbsTol ");
            this.SolverType= solvers.SolverTypes.Implicit;
            /*  "Disable" MaxStep DPCM warning as implicit solvers are stable */
            this.MaxStep= [];
        }

    
    protected:

        function varargout = solverCall(function_handle odefun,rowvec<double> t,x0,struct opts) {
            opts = odeset(opts, " RelTol ", this.RelTol, " AbsTol ", this.AbsTol);

            if isempty(opts.InitialSlope)
                error(" The ode15i solver must be provided with an initial slope. ");
            end

            /* % Use properties from AJacobianSolver
             * Set Jacobian or Mass matrix */
            if ~isempty(this.JacFun) || ~isempty(this.M)
                 opts = odeset(opts, " Jacobian ", @this.FJAC);
            end

            /*  Process any sparsity patterns */
            JP = [[],[]];
            if ~isempty(this.JPattern)
                JP[1] = this.JPattern;
            end
            /*  Set df/dyp sparsity pattern, derived from mass matrix. Works
             * only (at least can be guaranteed) for constant mass matrices. */
            if ~isempty(this.M)
                JP[2] = this.M.SparsityPattern;
            else
                JP[2] = speye(size(x0,1));
            end
            opts = odeset(opts, " JPattern ", JP);

            /* % Call implicit solver
             * implfun: A handle to the implicit function `f` which describes
             * the ODE via `f(t,x,x') = 0` */
            if ~isempty(this.M)
                if this.M.TimeDependent
                    implfun = @(t,x,xp)this.M.evaluate(t)*xp - odefun(t,x);
                else
                    M = this.M.evaluate(0);
                    implfun = @(t,x,xp)M*xp - odefun(t,x);
                end
            else
                implfun = @(t,x,xp)xp - odefun(t,x);
            end

            /*  Final check for consistent initial condition */
            initialnorm = norm(implfun(0,x0,opts.InitialSlope));
            if initialnorm > 1e-11
                warning(" Initial conditions possibly inconsistent. ||f(0,x0,dx0)|| = %g ",initialnorm);
            end

            /*  Call ode15i solver */
            [varargout[1:nargout]] = this.MLSolver(implfun, t, x0, opts.InitialSlope, opts);
        }
    private:

        function [dfdx , dfdxp ] = FJAC(t,x,xp) {

            /*  Implicit funcion is M*xp - odefun(t,x), so derivatives: */
            dfdx = [];
            if ~isempty(this.JacFun)
                dfdx = -this.JacFun(t, x);
            end
            if ~isempty(this.M)
                dfdxp = this.M.evaluate(t);
            else
                dfdxp = speye(size(x,1));
            end
        }
};
}