SemiImplicitEuler.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 SemiImplicitEuler
  :public solvers.BaseCustomSolver {
    private:

        model;

        
    public:

        SemiImplicitEuler(models.BaseFullModel model) {
            this.model= model;
            if ~isa(model.System," models.BaseFirstOrderSystem ")
                error(" Solver not implemented for non-first order systems. ");
            end
            this.Name= " Semi-implicit euler method ";
            this.SolverType= solvers.SolverTypes.Implicit;
        }

    
    protected:

        function matrix<double>x = customSolve(unused1,rowvec<double> t,colvec<double> x0,rowvec<integer> outputtimes) {
            s = this.model.System;
            if isempty(s.A)
               error(" This solver requires an (affine) linear system component A. "); 
            end
/*              if isa(s.Model,'models.ReducedModel') && ~isempty(s.Model.ErrorEstimator)...
 *                     && s.Model.ErrorEstimator.Enabled && s.Model.ErrorEstimator.ExtraODEDims > 0
 *                 error('Cannot use this solver with reduced models that have an error estimator enabled with ExtraODEDims > 0 (this solver overrides the odefun handle)');
 *             end
             * Initialize result */
            steps = length(t);

            dt = t(2)-t(1);
            if ~isempty(find((abs(t(2:end)-t(1:end-1) - dt)) / dt > 1e-6,1)) /* any(t(2:end)-t(1:end-1) - dt > 100*eps) */

                error(" non-equidistant dt timesteps. ");
            end

            rtm = this.RealTimeMode;
            /*  Initialize output index counter */
            outidx = 2;
            if rtm
                ed = solvers.SolverEventData;
                x = [];
            else
                effsteps = length(outputtimes);
                /*  Create return matrix in size of effectively desired timesteps */
                x = [x0 zeros(size(x0,1),effsteps-1)];
            end

            oldex = []; newex = []; edim = 0; est = [];
            if isa(this.model," models.ReducedModel ")
                est = this.model.ErrorEstimator;
                if ~isempty(est) && est.Enabled
                    edim = est.ExtraODEDims;
                    oldex = x0(end-edim+1:end);
                end
            end

            /*  Check if a mass matrix is present, otherwise assume identity matrix */
            M = s.getMassMatrix;
            mdep = M.TimeDependent;
            if ~mdep
                M = M.evaluate(0);
            end

            fdep = s.A.TimeDependent;
            if ~fdep
                /*  Evaluation with x=1 "extracts" the matrix A of the (affine) linear system */
                A = s.A.getStateJacobian;
            end

            /*  Precompute lu decomp for iteration steps if all components are not time-dependent */
            if ~mdep && ~fdep
                [l, u] = lu(M - dt * A);
            end

            /*  Solve for each time step */
            oldx = x0(1:end-edim);
            for idx = 2:steps;
                RHS = M*oldx;
                if ~isempty(s.f)
                    RHS = RHS + dt*s.f.evaluate(oldx, t(idx-1));
                end
                if ~isempty(s.B)
                    RHS = RHS + dt*s.B.evaluate(t(idx-1), s.mu)*s.u(t(idx-1));
                end

                /*  Time-independent case: constant */
                if ~fdep && ~mdep
                    newx = u\(l\RHS);
                    /*  Time-dependent case: compute time-dependent components with updated time */
                else
                    if fdep
                        A = s.A.evaluate(1, t(idx));
                    end
                    if mdep
                        M = s.M.evaluate(t(idx));
                    end
                    newx = (M - dt * A)\RHS;
                end

                /*  Implicit error estimation computation */
                if ~isempty(oldex)
                    ut = [];
                    if ~isempty(s.u)
                        ut = s.u(t(idx));
                    end
                    /* al = est.getAlpha(newx, t(idx), s.mu, ut);
                     *bet = est.getBeta(newx, t(idx), s.mu); */

                    /*  Explicit */
                    odepart = est.evalODEPart([newx; oldex], t(idx-1), ut);
                    newex = oldex + dt * odepart;
                    /* newex = oldex + dt*(bet*oldex + al); */

                    /*  Implicit
                     *newex = (dt*al+oldex)/(1-dt*bet); */

/*                      fun = @(y)y-dt*est.evalODEPart([oldx; y], t(idx), s.mu, ut)-oldex/dt;
 *                     opts = optimset('Display','off');
 *                     newex2 = fsolve(fun, oldex, opts); */
                end

                /*  Only produce output at wanted timesteps */
                if outputtimes(outidx) == idx
                    if rtm                        
                        /*  Real time mode: Fire StepPerformed event */
                        ed.Times= t(idx);
                        ed.States= [newx; newex];
                        this.notify(" StepPerformed ",ed);
                        /*  Normal mode: Collect solution in result vector */
                    else
                        x(:,outidx) = [newx; newex];/* #ok */

                    end
                    outidx = outidx+1;
                end
                oldx = newx;
                oldex = newex;
            end
        }
};
}