From 829d7c7ef609aebb52ac7055170edb6767d65d96 Mon Sep 17 00:00:00 2001 From: rjg18 Date: Fri, 4 Jun 2021 22:57:25 -0400 Subject: [PATCH 01/10] MATLODEV2 ERK and Infrastructure Added --- CHANGE_LOG.md | 7 ++ README.md | 27 +++++ src/+matlode/+rk/+erk/ClassicRK4.m | 31 ++++++ src/+matlode/+rk/+erk/DormandPrince.m | 34 +++++++ src/+matlode/+rk/+erk/ERK.m | 99 +++++++++++++++++++ src/+matlode/+rk/+erk/ForwardEuler.m | 27 +++++ src/+matlode/+rk/RungeKutta.m | 45 +++++++++ src/+matlode/+stepsizecontroller/Fixed.m | 34 +++++++ .../+stepsizecontroller/StepSizeController.m | 15 +++ src/+matlode/+util/+errnorm/standardNorm.m | 9 ++ src/+matlode/Integrator.m | 53 ++++++++++ src/+matlode/OneStepIntegrator.m | 23 +++++ 12 files changed, 404 insertions(+) create mode 100644 CHANGE_LOG.md create mode 100644 README.md create mode 100644 src/+matlode/+rk/+erk/ClassicRK4.m create mode 100644 src/+matlode/+rk/+erk/DormandPrince.m create mode 100644 src/+matlode/+rk/+erk/ERK.m create mode 100644 src/+matlode/+rk/+erk/ForwardEuler.m create mode 100644 src/+matlode/+rk/RungeKutta.m create mode 100644 src/+matlode/+stepsizecontroller/Fixed.m create mode 100644 src/+matlode/+stepsizecontroller/StepSizeController.m create mode 100644 src/+matlode/+util/+errnorm/standardNorm.m create mode 100644 src/+matlode/Integrator.m create mode 100644 src/+matlode/OneStepIntegrator.m diff --git a/CHANGE_LOG.md b/CHANGE_LOG.md new file mode 100644 index 0000000..d6002b0 --- /dev/null +++ b/CHANGE_LOG.md @@ -0,0 +1,7 @@ + + +### 2.0.0beta1: ++ Change log ++ ERK, RungeKutta, DormanPrince, ClassicRK4, ForwardEuler ++ Fixed +- Removed MATLODE v1 \ No newline at end of file diff --git a/README.md b/README.md new file mode 100644 index 0000000..ded399f --- /dev/null +++ b/README.md @@ -0,0 +1,27 @@ + +# MATLODE V2.0.0beta1 + +##Example Code +*Using ODE Test Problems + +```matlab +problem = otp.hires.presets.Canonical; +integrator = matlode.rk.erk.ClassicRK4; +options = integrator.matlodeSets('StepsN', 100000); +sol = integrator.integrate(problem.Rhs.F, problem.TimeSpan, problem.Y0); +``` +OR +```matlab +problem = otp.hires.presets.Canonical; +integrator = matlode.rk.erk.ClassicRK4; +sol = integrator.integrate(problem.Rhs.F, problem.TimeSpan, problem.Y0, 'StepsN', 100000); +``` +OR +```matlab +problem = otp.hires.presets.Canonical; +integrator = matlode.rk.erk.ClassicRK4; +options = integrator.matlodeSets('StepsN', 100000); +sol = integrator.integrate(problem.Rhs.F, problem.TimeSpan, problem.Y0, options, 'AbsTol', 1e-8); +``` + +[insert other useful stuff] \ No newline at end of file diff --git a/src/+matlode/+rk/+erk/ClassicRK4.m b/src/+matlode/+rk/+erk/ClassicRK4.m new file mode 100644 index 0000000..245b82f --- /dev/null +++ b/src/+matlode/+rk/+erk/ClassicRK4.m @@ -0,0 +1,31 @@ +classdef ClassicRK4 < matlode.rk.erk.ERK + methods + function obj = ClassicRK4() + + %TODO + %include reference + + a = [[0, 0, 0, 0]; + [1/2, 0, 0, 0]; + [0, 1/2, 0, 0]; + [0, 0, 0, 1]]; + + b = [1/6, 1/3, 1/3, 1/6]; + + bHat = []; + + c = [0, 1/2, 1/2, 1]; + + order = 4; + + embbededOrder = 0; + + obj = obj@matlode.rk.erk.ERK(a, b, bHat, c, order, embbededOrder); + end + + function sol = integrate(obj, f, tspan, y0, varargin) + sol = integrate@matlode.rk.erk.ERK(obj, f, tspan, y0, 'StepSizeController', matlode.stepsizecontroller.Fixed(1000), varargin{:}); + end + end +end + diff --git a/src/+matlode/+rk/+erk/DormandPrince.m b/src/+matlode/+rk/+erk/DormandPrince.m new file mode 100644 index 0000000..eb7215c --- /dev/null +++ b/src/+matlode/+rk/+erk/DormandPrince.m @@ -0,0 +1,34 @@ +classdef DormandPrince < matlode.rk.erk.ERK + methods + function obj = DormandPrince() + + %TODO + %include reference + + a = [[0, 0, 0, 0, 0, 0, 0]; + [1/5, 0, 0, 0, 0, 0, 0]; + [3/40, 9/40, 0, 0, 0, 0, 0]; + [44/45, -56/15, 32/9, 0, 0, 0, 0]; + [19372/6561, -25360/2187, 64448/6561, -212/729, 0, 0, 0]; + [9017/3168, -355/33, 46732/5247, 49/176, -5103/18656, 0, 0]; + [35/384, 0, 500/1113, 125/192, -2187/6784, 11/84, 0]]; + + b = [35/384, 0, 500/1113, 125/192, -2187/6784, 11/84, 0]; + + bHat = [5179/57600, 0, 7571/16695, 393/640, -92097/339200, 187/2100, 1/40]; + + c = [0, 1/5, 3/10, 4/5, 8/9, 1, 1]; + + order = 5; + + embbededOrder = 4; + + obj = obj@matlode.rk.erk.ERK(a, b, bHat, c, order, embbededOrder); + end + + function sol = integrate(obj, f, tspan, y0, varargin) + sol = integrate@matlode.rk.erk.ERK(obj, f, tspan, y0, 'StepSizeController', matlode.stepsizecontroller.Fixed(1000), varargin{:}); + end + end +end + diff --git a/src/+matlode/+rk/+erk/ERK.m b/src/+matlode/+rk/+erk/ERK.m new file mode 100644 index 0000000..745491d --- /dev/null +++ b/src/+matlode/+rk/+erk/ERK.m @@ -0,0 +1,99 @@ +classdef ERK < matlode.rk.RungeKutta + %Explicit Runge Kutta class + + properties (SetAccess = immutable) + a + b + bHat + c + stage + order + embeddedOrder + adaptable + FSAL + end + + methods + + + + function obj = ERK(a, b, bHat, c, order, embeddedOrder) + + obj.a = a; + obj.b = b; + obj.c = c; + obj.bHat = bHat; + obj.embeddedOrder = embeddedOrder; + obj.order = order; + obj.stage = length(b); + obj.FSAL = all(a(end, :) == b) && all(a(1, :) == 0); + obj.adaptable = ~isempty(bHat); + end + + function sol = integrate(obj, f, tspan, y0, varargin) + + opts = obj.matlodeSets(varargin{:}); + + numVars = length(y0); + t = tspan(1); + k = zeros(numVars, obj.stage); + sol.y = y0; + + sol.stats.nsteps = 0; + sol.stats.nfailed = 0; + sol.stats.nfevals = 0; + + fsalStart = uint32(obj.FSAL) + 1; + fevalIterCounts = obj.stage - fsalStart + 1; + + + [h, k(:,1)] = opts.StepSizeController.startingStep(f, tspan, y0, obj.order, opts); + + + while t < tspan(end) + + for i = fsalStart:obj.stage + k(:, i) = f(t + h * obj.c(i), sol.y + h * (k(:, 1:i-1) * obj.a(i, 1:i-1)')); + end + + sol.stats.nfevals = sol.stats.nfevals + fevalIterCounts; + + yCurrent = sol.y + h * (k * obj.b'); + err = 0; + + + if opts.StepSizeController.adaptive + yEmbbeded = sol.y + h * (k * obj.bHat); + err = opts.ErrNorm(yCurrent, yEmbbeded, opts); + end + + [accept, h, t] = opts.StepSizeController.newStepSize(t, tspan, h, err, obj.embeddedOrder, sol.stats, opts); + + if accept + sol.y = yCurrent; + sol.stats.nsteps = sol.stats.nsteps + 1; + + + if obj.FSAL + k(:, 1) = k(:, end); + end + else + + sol.stats.nfailed = sol.stats.nfailed + 1; + end + end + + + end + + function opts = matlodeSets(obj, varargin) + p = inputParser; + + %ERK specific options + + opts = obj.matlodeSetsHelper(p, varargin{:}); + end + end + +end + diff --git a/src/+matlode/+rk/+erk/ForwardEuler.m b/src/+matlode/+rk/+erk/ForwardEuler.m new file mode 100644 index 0000000..1218a2d --- /dev/null +++ b/src/+matlode/+rk/+erk/ForwardEuler.m @@ -0,0 +1,27 @@ +classdef ForwardEuler < matlode.rk.erk.ERK + methods + function obj = ForwardEuler() + %TODO + %include reference + + a = [[0]]; + + b = [1]; + + bHat = []; + + c = [0]; + + order = 1; + + embbededOrder = 0; + + obj = obj@matlode.rk.erk.ERK(a, b, bHat, c, order, embbededOrder); + end + + function sol = integrate(obj, f, tspan, y0, varargin) + sol = integrate@matlode.rk.erk.ERK(obj, f, tspan, y0, 'StepSizeController', matlode.stepsizecontroller.Fixed(1000), varargin{:}); + end + end +end + diff --git a/src/+matlode/+rk/RungeKutta.m b/src/+matlode/+rk/RungeKutta.m new file mode 100644 index 0000000..a08eee4 --- /dev/null +++ b/src/+matlode/+rk/RungeKutta.m @@ -0,0 +1,45 @@ +classdef (Abstract) RungeKutta < matlode.OneStepIntegrator + %RungeKutta integrator + + properties (Abstract, SetAccess = immutable) + bHat + embeddedOrder + order + FSAL + end + + methods + function obj = fromCoeffs(a, b, bHat, c, order, embededOrder) + + + if istril(a) + if any(diag(a)) + %TODO + % set as DIRK + return + end + + obj = ERK(a, b, bHat, c, order, embededOrder); + return + end + + %TODO + + %Set as FIRK + + end + end + + methods (Access = protected) + function opts = matlodeSetsHelper(obj, p, varargin) + + %RungeKutta sepcific options + + opts = matlodeSetsHelper@matlode.OneStepIntegrator(obj, p, varargin{:}); + + end + end + + +end + diff --git a/src/+matlode/+stepsizecontroller/Fixed.m b/src/+matlode/+stepsizecontroller/Fixed.m new file mode 100644 index 0000000..e65b775 --- /dev/null +++ b/src/+matlode/+stepsizecontroller/Fixed.m @@ -0,0 +1,34 @@ +classdef Fixed < matlode.stepsizecontroller.StepSizeController + %Fixed stepsize controller + + properties (Constant) + adaptive = false; + end + + properties (SetAccess = immutable, GetAccess = private) + steps + end + + methods + + function obj = Fixed(steps) + obj.steps = steps; + end + + function [h0, f0] = startingStep(obj, f, tspan, y0, ~, options) + if options.InitialStep ~= 0 + warning('Ignoring initialStep option, Controller is Fixed'); + end + + f0 = f(tspan(1), y0); + h0 = (tspan(2) - tspan(1)) / obj.steps; + end + + function [accept, hCur, tNew] = newStepSize(~, ~, tspan, hCur, ~, ~, stepdlx, ~) + accept = true; + tNew = (stepdlx.nsteps + 1) * hCur + tspan(1); + end + + end +end + diff --git a/src/+matlode/+stepsizecontroller/StepSizeController.m b/src/+matlode/+stepsizecontroller/StepSizeController.m new file mode 100644 index 0000000..b6f45d0 --- /dev/null +++ b/src/+matlode/+stepsizecontroller/StepSizeController.m @@ -0,0 +1,15 @@ +classdef (Abstract) StepSizeController < handle + %Template for a step size controller + + properties (Abstract, Constant) + adaptive + end + + methods (Abstract) + + [h0, f0] = startingStep(obj, f, tspan, y0, order, options); + [accept, hNew, tNew] = newStepSize(obj, tCur, tspan, hCur, err, embeddedOrder, stepdlx, options); + end + +end + diff --git a/src/+matlode/+util/+errnorm/standardNorm.m b/src/+matlode/+util/+errnorm/standardNorm.m new file mode 100644 index 0000000..e711a8d --- /dev/null +++ b/src/+matlode/+util/+errnorm/standardNorm.m @@ -0,0 +1,9 @@ +function err = standardNorm(y, yHat, options) + %procedure described in solving ODEs I Harrier + %cannot assume that user has signal processing toolbox + + sc = (options.AbsTol + max(abs(y), abs(yHat)) .* options.RelTol); + + err = sqrt(1/length(y)) * norm((y - yHat) ./ sc); +end + diff --git a/src/+matlode/Integrator.m b/src/+matlode/Integrator.m new file mode 100644 index 0000000..5d341c6 --- /dev/null +++ b/src/+matlode/Integrator.m @@ -0,0 +1,53 @@ +classdef (Abstract) Integrator < handle + %INTEGRATOR Summary of this class goes here + % Detailed explanation goes here + + properties (Abstract, SetAccess = immutable) + c + adaptable + end + + methods (Abstract) + + %integrate + sol = integrate(obj, func, tspan, y0, varargin) + opts = matlodeSets(obj, varargin) + + end + + methods (Access = protected) + + + function opts = matlodeSetsHelper(obj, p, varargin) + + p.addParameter('AbsTol', 1e-6); + p.addParameter('RelTol', 1e-6); + p.addParameter('InitialStep', 0); + p.addParameter('StepSizeController', matlode.stepsizecontroller.Fixed(1000)); + p.addParameter('StepsN', 0); + p.addParameter('ErrNorm', @(y, yHat, options) matlode.util.errnorm.standardNorm(y, yHat, options)); + p.addParameter('NewtonTol', 1e-7); + p.addParameter('NewtonMaxItr', 25); + p.addParameter('Jacobian', []); + + p.parse(varargin{:}); + + opts = p.Results; + + %shorthanded fixed steps + if opts.StepsN > 0 + opts.StepSizeController = matlode.stepsizecontroller.Fixed(opts.StepsN); + end + opts = rmfield(opts, 'StepsN'); + + %make sure a method can use a adaptive controller + if ~obj.adaptable && opts.StepSizeController.adaptive + warning('Current method does not have the ability to use a adaptive step size controller'); + %overrider user choice??? + end + + end + end + +end + diff --git a/src/+matlode/OneStepIntegrator.m b/src/+matlode/OneStepIntegrator.m new file mode 100644 index 0000000..f6653e6 --- /dev/null +++ b/src/+matlode/OneStepIntegrator.m @@ -0,0 +1,23 @@ +classdef (Abstract) OneStepIntegrator < matlode.Integrator + %One Step Integrator template + + properties (Abstract, SetAccess = immutable) + a + b + stage + end + + methods (Access = protected) + function opts = matlodeSetsHelper(obj, p, varargin) + + + + %OneStepIntegrator Specific options + + opts = matlodeSetsHelper@matlode.Integrator(obj, p, varargin{:}); + + end + + end +end + From b16d3f4f366ea84bfd0a14f1f79f8fb33606d722 Mon Sep 17 00:00:00 2001 From: Steven Roberts Date: Wed, 23 Jun 2021 10:47:29 -0400 Subject: [PATCH 02/10] Add linear solver framework --- .gitignore | 7 +++ .../+linearsolver/DecompositionLinearSolver.m | 22 +++++++++ src/+matlode/+linearsolver/LinearSolver.m | 8 +++ .../+linearsolver/MatrixFreeLinearSolver.m | 34 +++++++++++++ .../+linearsolver/MatrixLinearSolver.m | 49 +++++++++++++++++++ 5 files changed, 120 insertions(+) create mode 100644 .gitignore create mode 100644 src/+matlode/+linearsolver/DecompositionLinearSolver.m create mode 100644 src/+matlode/+linearsolver/LinearSolver.m create mode 100644 src/+matlode/+linearsolver/MatrixFreeLinearSolver.m create mode 100644 src/+matlode/+linearsolver/MatrixLinearSolver.m diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..b352972 --- /dev/null +++ b/.gitignore @@ -0,0 +1,7 @@ +*.asv +*.m~ +*.mex* +*.mlappinstall +*.mltbx +octave-workspace + diff --git a/src/+matlode/+linearsolver/DecompositionLinearSolver.m b/src/+matlode/+linearsolver/DecompositionLinearSolver.m new file mode 100644 index 0000000..f0d58bb --- /dev/null +++ b/src/+matlode/+linearsolver/DecompositionLinearSolver.m @@ -0,0 +1,22 @@ +classdef DecompositionLinearSolver < matlode.linearsolver.MatrixLinearSolver + properties (SetAccess = immutable, GetAccess = private) + Args + end + + methods + function obj = DecompositionLinearSolver(varargin) + obj = obj@matlode.linearsolver.MatrixLinearSolver(); + obj.Args = varargin; + end + + function system = preprocess(obj, updateState, updateTimestep, system, t, y, f, m1, mass, m2, jac) + system = preprocess@matlode.linearsolver.MatrixLinearSolver( ... + obj, updateState, updateTimestep, system, t, y, f, m1, mass, m2, jac); + + if updateState || updateTimestep + system.matrix = decomposition(system.matrix, obj.Args{:}); + end + end + end +end + diff --git a/src/+matlode/+linearsolver/LinearSolver.m b/src/+matlode/+linearsolver/LinearSolver.m new file mode 100644 index 0000000..e5aee67 --- /dev/null +++ b/src/+matlode/+linearsolver/LinearSolver.m @@ -0,0 +1,8 @@ +classdef (Abstract) LinearSolver + methods (Abstract) + system = preprocess(obj, updateState, updateTimestep, oldSystem, t, y, f, m1, mass, m2, jac); + + sol = solve(obj, x, system, t, y, f, m1, mass, m2, jac); + end +end + diff --git a/src/+matlode/+linearsolver/MatrixFreeLinearSolver.m b/src/+matlode/+linearsolver/MatrixFreeLinearSolver.m new file mode 100644 index 0000000..311f535 --- /dev/null +++ b/src/+matlode/+linearsolver/MatrixFreeLinearSolver.m @@ -0,0 +1,34 @@ +classdef MatrixFreeLinearSolver < matlode.linearsolver.LinearSolver + properties (SetAccess = immutable, GetAccess = private) + Solver + end + + methods + function obj = MatrixFreeLinearSolver(solver) + if nargin < 1 + solver = @gmres; + end + + obj.Solver = solver; + end + + function system = preprocess(~, ~, ~, system, ~, ~, ~, ~, ~, ~, ~) + end + + function sol = solve(obj, x, ~, t, y, ~, m1, mass, m2, jac) + if length(m1) == 1 + if isempty(jac) + jvp = @(v) m1 * mass(t, y, v); + elseif isempty(mass) + jvp = @(v) m1 * v - m2 * jac(t, y, v); + else + jvp = @(v) m1 * mass(t, y, v) - m2 * jac(t, y, v); + end + else + error('Matrix-free block systems are not supported yet'); + end + sol = obj.Solver(jvp, x); + end + end +end + diff --git a/src/+matlode/+linearsolver/MatrixLinearSolver.m b/src/+matlode/+linearsolver/MatrixLinearSolver.m new file mode 100644 index 0000000..f036bc8 --- /dev/null +++ b/src/+matlode/+linearsolver/MatrixLinearSolver.m @@ -0,0 +1,49 @@ +classdef MatrixLinearSolver < matlode.linearsolver.LinearSolver + properties (SetAccess = immutable, GetAccess = private) + Solver + end + + methods + function obj = MatrixLinearSolver(solver) + if nargin < 1 + solver = @mldivide; + end + + obj.Solver = solver; + end + + function system = preprocess(~, updateState, updateTimestep, system, t, y, ~, m1, mass, m2, jac) + init = isempty(system); + + if updateState + if ~isnumeric(jac) + system.jac = jac(t, y); + elseif init + system.jac = jac; + end + + if ~isnumeric(mass) + system.mass = mass(t, y); + elseif init + if isempty(mass) + system.mass = eye(size(system.jac), 'like', system.jac); + else + system.mass = mass; + end + end + end + + if updateState || updateTimestep + system.matrix = kron(m1, system.mass); + if m2 ~= 0 + system.matrix = system.matrix - kron(m2, jac); + end + end + end + + function sol = solve(obj, x, system, ~, ~, ~, ~, ~, ~, ~) + sol = obj.Solver(system.matrix, x); + end + end +end + From d35bbf5a82dd9ae77bbe0e3b0a555cb1021688af Mon Sep 17 00:00:00 2001 From: rjg18 Date: Tue, 29 Jun 2021 14:26:40 -0400 Subject: [PATCH 03/10] Added StepSizeControllers, ErrNorm, DenseOutput, and Optimizations --- CHANGE_LOG.md | 23 +- README.md | 22 +- src/+matlode/+denseoutput/DenseOutput.m | 6 + src/+matlode/+denseoutput/Direct.m | 24 ++ src/+matlode/+denseoutput/Hermite.m | 18 + src/+matlode/+denseoutput/Linear.m | 20 ++ src/+matlode/+denseoutput/RKDense.m | 22 ++ src/+matlode/+errnorm/HistNorm.m | 11 + src/+matlode/+errnorm/InfNorm.m | 11 + src/+matlode/+errnorm/NormErr.m | 6 + src/+matlode/+errnorm/StandardNorm.m | 11 + src/+matlode/+rk/+erk/ClassicRK4.m | 38 ++- src/+matlode/+rk/+erk/DormandPrince.m | 58 +++- src/+matlode/+rk/+erk/ERK.m | 315 +++++++++++++++--- src/+matlode/+rk/+erk/ForwardEuler.m | 24 +- src/+matlode/+rk/RungeKutta.m | 19 +- .../+stepsizecontroller/+soberpresets/H321.m | 17 + src/+matlode/+stepsizecontroller/Fixed.m | 21 +- src/+matlode/+stepsizecontroller/Gustafson.m | 70 ++++ .../+stepsizecontroller/SoberlandController.m | 78 +++++ .../+stepsizecontroller/StandardController.m | 49 +++ .../+stepsizecontroller/StartingController.m | 41 +++ .../+stepsizecontroller/StepSizeController.m | 13 +- src/+matlode/+util/+errnorm/standardNorm.m | 9 - src/+matlode/+util/CoefficentTransformers.m | 15 + src/+matlode/Integrator.m | 93 ++++-- src/+matlode/OneStepIntegrator.m | 11 +- 27 files changed, 881 insertions(+), 164 deletions(-) create mode 100644 src/+matlode/+denseoutput/DenseOutput.m create mode 100644 src/+matlode/+denseoutput/Direct.m create mode 100644 src/+matlode/+denseoutput/Hermite.m create mode 100644 src/+matlode/+denseoutput/Linear.m create mode 100644 src/+matlode/+denseoutput/RKDense.m create mode 100644 src/+matlode/+errnorm/HistNorm.m create mode 100644 src/+matlode/+errnorm/InfNorm.m create mode 100644 src/+matlode/+errnorm/NormErr.m create mode 100644 src/+matlode/+errnorm/StandardNorm.m create mode 100644 src/+matlode/+stepsizecontroller/+soberpresets/H321.m create mode 100644 src/+matlode/+stepsizecontroller/Gustafson.m create mode 100644 src/+matlode/+stepsizecontroller/SoberlandController.m create mode 100644 src/+matlode/+stepsizecontroller/StandardController.m create mode 100644 src/+matlode/+stepsizecontroller/StartingController.m delete mode 100644 src/+matlode/+util/+errnorm/standardNorm.m create mode 100644 src/+matlode/+util/CoefficentTransformers.m diff --git a/CHANGE_LOG.md b/CHANGE_LOG.md index d6002b0..699fe6d 100644 --- a/CHANGE_LOG.md +++ b/CHANGE_LOG.md @@ -1,7 +1,20 @@ -### 2.0.0beta1: -+ Change log -+ ERK, RungeKutta, DormanPrince, ClassicRK4, ForwardEuler -+ Fixed -- Removed MATLODE v1 \ No newline at end of file +### 2.0.0beta2: ### +\+ StepSizeControllers: Standard, Gustafson, Soberland, H321 +\+ DenseOutput: DirectFunctions, RKDense, Linear, Hermite +\+ ErrNorm: Standard, History, Inf +\+ CoefficentCasting +\+ IntegratTo +\~ Updated Options +\~ Improved Time Loop Structure +\~ Optimizations +\~ DenseOuput Coefficents +\- Removed MatlodeSets and some option changes + + +### 2.0.0beta1: ### +\+ Change log +\+ ERK, RungeKutta, DormanPrince, ClassicRK4, ForwardEuler +\+ Fixed +\- Removed MATLODE v1 \ No newline at end of file diff --git a/README.md b/README.md index ded399f..e6383ea 100644 --- a/README.md +++ b/README.md @@ -6,22 +6,12 @@ ```matlab problem = otp.hires.presets.Canonical; -integrator = matlode.rk.erk.ClassicRK4; -options = integrator.matlodeSets('StepsN', 100000); -sol = integrator.integrate(problem.Rhs.F, problem.TimeSpan, problem.Y0); -``` -OR -```matlab -problem = otp.hires.presets.Canonical; -integrator = matlode.rk.erk.ClassicRK4; -sol = integrator.integrate(problem.Rhs.F, problem.TimeSpan, problem.Y0, 'StepsN', 100000); -``` -OR -```matlab -problem = otp.hires.presets.Canonical; -integrator = matlode.rk.erk.ClassicRK4; -options = integrator.matlodeSets('StepsN', 100000); -sol = integrator.integrate(problem.Rhs.F, problem.TimeSpan, problem.Y0, options, 'AbsTol', 1e-8); +integrator = matlode.rk.erk.DormanPrince; + +options.ErrNorm = matlode.errnorm.HistNorm.errEstimate(1e-6, 1e-6); +options.StepSizeController = matlode.stepsizecontroller.Gustafson; + +sol = integrator.integrate(problem.Rhs.F, problem.TimeSpan, problem.Y0, options); ``` [insert other useful stuff] \ No newline at end of file diff --git a/src/+matlode/+denseoutput/DenseOutput.m b/src/+matlode/+denseoutput/DenseOutput.m new file mode 100644 index 0000000..f904736 --- /dev/null +++ b/src/+matlode/+denseoutput/DenseOutput.m @@ -0,0 +1,6 @@ +classdef (Abstract) DenseOutput < handle + methods (Abstract) + [denseY, fEvals] = denseOut(obj, f, t, tNeeded, yC, yN, k, h) + end +end + diff --git a/src/+matlode/+denseoutput/Direct.m b/src/+matlode/+denseoutput/Direct.m new file mode 100644 index 0000000..c879b2e --- /dev/null +++ b/src/+matlode/+denseoutput/Direct.m @@ -0,0 +1,24 @@ +classdef Direct < matlode.denseoutput.DenseOutput + %This class is for custom denseoutputs that methods could require + %The custum dense output cannot require f evaluations + %Function must be in row format + + properties (SetAccess = immutable) + kFuncs + end + + methods + function obj = Direct(kfuncs) + obj.kFuncs = kfuncs; + end + + function [denseY, fEvals] = denseOut(obj, ~, t, tNeeded, yC, ~, k, h) + + theta = (tNeeded - t) / h; + + denseY = yC + (k * h) * obj.kFuncs(theta); + fEvals = 0; + end + end +end + diff --git a/src/+matlode/+denseoutput/Hermite.m b/src/+matlode/+denseoutput/Hermite.m new file mode 100644 index 0000000..90a2281 --- /dev/null +++ b/src/+matlode/+denseoutput/Hermite.m @@ -0,0 +1,18 @@ +classdef Hermite < matlode.denseoutput.DenseOutput + %HERMITE Formula obtained from Solving ODEs I Chapter II.6 + + methods + function obj = Hermite() + end + + function [denseY, fEvals] = denseOut(~, f, t, tNeeded, yC, yN, ~, h) + + theta = (tNeeded - t) / h; + f0 = f(t, yC); + f1 = f(t + h, yN); + denseY = (1 - theta) * yC + theta * yN + theta*(theta - 1) * ((1-2*theta) * (yN - yC) + (theta - 1) * (h * f0 + theta * h * f1)); + fEvals = 2; + end + end +end + diff --git a/src/+matlode/+denseoutput/Linear.m b/src/+matlode/+denseoutput/Linear.m new file mode 100644 index 0000000..acfed0f --- /dev/null +++ b/src/+matlode/+denseoutput/Linear.m @@ -0,0 +1,20 @@ +classdef Linear < matlode.denseoutput.DenseOutput + + properties (SetAccess = immutable) + b + end + + methods + function obj = Linear(b) + obj.b = b; + end + + function [denseY, fEvals] = denseOut(obj, ~, t, tNeeded, yC, ~, k, h) + + theta = (tNeeded - t) / h; + denseY = yC + (k * h) * (obj.b' * theta) ; + fEvals = 0; + end + end +end + diff --git a/src/+matlode/+denseoutput/RKDense.m b/src/+matlode/+denseoutput/RKDense.m new file mode 100644 index 0000000..d28301a --- /dev/null +++ b/src/+matlode/+denseoutput/RKDense.m @@ -0,0 +1,22 @@ +classdef RKDense < matlode.denseoutput.DenseOutput + %RKDENSE Used for dense output computed from RK trees + + properties (SetAccess = immutable) + bTheta + end + + methods + function obj = RKDense(b) + obj.bTheta = @(theta) b * theta.^(1:size(b,2))'; + end + + function [denseY, fEvals] = denseOut(obj, ~, t, tNeeded, yC, ~, k, h) + + theta = (tNeeded - t) / h; + + denseY = yC + (k * h) * obj.bTheta(theta); + fEvals = 0; + end + end +end + diff --git a/src/+matlode/+errnorm/HistNorm.m b/src/+matlode/+errnorm/HistNorm.m new file mode 100644 index 0000000..6c44923 --- /dev/null +++ b/src/+matlode/+errnorm/HistNorm.m @@ -0,0 +1,11 @@ +classdef HistNorm < matlode.errnorm.NormErr + methods (Static) + function errFunc = errEstimate(AbsTol, RelTol) + + scFunc = @(y) (AbsTol + max(abs(y(:, 1)), abs(y(:, 2))) * RelTol); + errFunc = @(y, yE) sqrt(1/size(y,1)) * norm(yE ./ scFunc(y)); + + end + end +end + diff --git a/src/+matlode/+errnorm/InfNorm.m b/src/+matlode/+errnorm/InfNorm.m new file mode 100644 index 0000000..e83702a --- /dev/null +++ b/src/+matlode/+errnorm/InfNorm.m @@ -0,0 +1,11 @@ +classdef InfNorm < matlode.errnorm.NormErr + methods (Static) + function errFunc = errEstimate(AbsTol, RelTol) + + scFunc = @(y) (AbsTol + abs(y(:, 1)) * RelTol); + errFunc = @(y, yE) norm(yE ./ scFunc(y), 'Inf'); + + end + end +end + diff --git a/src/+matlode/+errnorm/NormErr.m b/src/+matlode/+errnorm/NormErr.m new file mode 100644 index 0000000..5f5a08a --- /dev/null +++ b/src/+matlode/+errnorm/NormErr.m @@ -0,0 +1,6 @@ +classdef (Abstract) NormErr + methods (Abstract,Static) + errFunc = errEstimate(AbsTol, RelTol) + end +end + diff --git a/src/+matlode/+errnorm/StandardNorm.m b/src/+matlode/+errnorm/StandardNorm.m new file mode 100644 index 0000000..c6eb454 --- /dev/null +++ b/src/+matlode/+errnorm/StandardNorm.m @@ -0,0 +1,11 @@ +classdef StandardNorm < matlode.errnorm.NormErr + methods (Static) + function errFunc = errEstimate(AbsTol, RelTol) + + scFunc = @(y) (AbsTol + abs(y(:, 1)) * RelTol); + errFunc = @(y, yE) sqrt(1/size(y,1)) * norm(yE ./ scFunc(y)); + + end + end +end + diff --git a/src/+matlode/+rk/+erk/ClassicRK4.m b/src/+matlode/+rk/+erk/ClassicRK4.m index 245b82f..ae72c45 100644 --- a/src/+matlode/+rk/+erk/ClassicRK4.m +++ b/src/+matlode/+rk/+erk/ClassicRK4.m @@ -1,26 +1,48 @@ classdef ClassicRK4 < matlode.rk.erk.ERK + properties (SetAccess = immutable) + DenseOut + end + methods - function obj = ClassicRK4() + function obj = ClassicRK4(datatype) %TODO %include reference - a = [[0, 0, 0, 0]; - [1/2, 0, 0, 0]; - [0, 1/2, 0, 0]; - [0, 0, 0, 1]]; + % p = 4 & hatp = 0 & s = 4 + + if nargin == 0 + datatype = 'double'; + end + + caster = @(x) matlode.util.CoefficentTransformers.transform(x,datatype); + + a = caster(join(['[[0, 0, 0, 0];',... + '[1/2, 0, 0, 0];',... + '[0, 1/2, 0, 0];',... + '[0, 0, 0, 1]]'])); - b = [1/6, 1/3, 1/3, 1/6]; + b = caster('[1/6, 1/3, 1/3, 1/6]'); bHat = []; - c = [0, 1/2, 1/2, 1]; + e = []; + + c = caster('[0, 1/2, 1/2, 1]'); order = 4; embbededOrder = 0; - obj = obj@matlode.rk.erk.ERK(a, b, bHat, c, order, embbededOrder); + obj = obj@matlode.rk.erk.ERK(a, b, bHat, c, e, order, embbededOrder); + + %p* = 3 & s* = 3 + denseMatrix = caster(join(['[1, -3/2, 2/3];',... + '[0, 1, -2/3]',... + '[0, 1, -2/3]',... + '[0, 1/2, 2/3]'])); + + obj.DenseOut = matlode.denseoutput.RKDense(denseMatrix); end function sol = integrate(obj, f, tspan, y0, varargin) diff --git a/src/+matlode/+rk/+erk/DormandPrince.m b/src/+matlode/+rk/+erk/DormandPrince.m index eb7215c..7a313a4 100644 --- a/src/+matlode/+rk/+erk/DormandPrince.m +++ b/src/+matlode/+rk/+erk/DormandPrince.m @@ -1,33 +1,61 @@ classdef DormandPrince < matlode.rk.erk.ERK + properties (SetAccess = immutable) + DenseOut + end + methods - function obj = DormandPrince() + function obj = DormandPrince(datatype) %TODO - %include reference - - a = [[0, 0, 0, 0, 0, 0, 0]; - [1/5, 0, 0, 0, 0, 0, 0]; - [3/40, 9/40, 0, 0, 0, 0, 0]; - [44/45, -56/15, 32/9, 0, 0, 0, 0]; - [19372/6561, -25360/2187, 64448/6561, -212/729, 0, 0, 0]; - [9017/3168, -355/33, 46732/5247, 49/176, -5103/18656, 0, 0]; - [35/384, 0, 500/1113, 125/192, -2187/6784, 11/84, 0]]; + %include reference for method + + if nargin == 0 + datatype = 'double'; + end + + caster = @(x) matlode.util.CoefficentTransformers.transform(x,datatype); + + + a = caster(join(['[[0, 0, 0, 0, 0, 0, 0];',... + '[1/5, 0, 0, 0, 0, 0, 0];',... + '[3/40, 9/40, 0, 0, 0, 0, 0];',... + '[44/45, -56/15, 32/9, 0, 0, 0, 0];',... + '[19372/6561, -25360/2187, 64448/6561, -212/729, 0, 0, 0];',... + '[9017/3168, -355/33, 46732/5247, 49/176, -5103/18656, 0, 0];',... + '[35/384, 0, 500/1113, 125/192, -2187/6784, 11/84, 0]]'])); - b = [35/384, 0, 500/1113, 125/192, -2187/6784, 11/84, 0]; + b = caster('[35/384, 0, 500/1113, 125/192, -2187/6784, 11/84, 0]'); + + bHat = caster('[5179/57600, 0, 7571/16695, 393/640, -92097/339200, 187/2100, 1/40]'); - bHat = [5179/57600, 0, 7571/16695, 393/640, -92097/339200, 187/2100, 1/40]; + e = caster('[71/57600, 0, -71/16695, 71/1920, -17253/339200, 22/525, -1/40]'); - c = [0, 1/5, 3/10, 4/5, 8/9, 1, 1]; + c = caster('[0, 1/5, 3/10, 4/5, 8/9, 1, 1]'); order = 5; embbededOrder = 4; - obj = obj@matlode.rk.erk.ERK(a, b, bHat, c, order, embbededOrder); + obj = obj@matlode.rk.erk.ERK(a, b, bHat, c, e, order, embbededOrder); + + %Shampien formula for dense output given in the Solving ODE + %book on pg192 (6.12) + %the Dense output has been transformed from a function form to + %a matrix form + %p* == 4 s* == 5 + denseMatrix = caster(join(['[[1, -(4034104133/1410260304), 105330401/33982176, -(13107642775/11282082432), 6542295/470086768];',... + '[ 0, 0, 0, 0, 0];',... + '[ 0, 132343189600/32700410799, -(833316000/131326951), 91412856700/32700410799, -(523383600/10900136933)];',... + '[ 0, -(115792950/29380423), 185270875/16991088, -(12653452475/1880347072), 98134425/235043384];',... + '[ 0, 70805911779/24914598704, -(4531260609/600351776), 988140236175/199316789632, -(14307999165/24914598704)];',... + '[ 0, -(331320693/205662961), 31361737/7433601, -(2426908385/822651844), 97305120/205662961];',... + '[ 0, 44764047/29380423, -(1532549/353981), 90730570/29380423, -(8293050/29380423)]]'])); + obj.DenseOut = matlode.denseoutput.RKDense(denseMatrix); + end function sol = integrate(obj, f, tspan, y0, varargin) - sol = integrate@matlode.rk.erk.ERK(obj, f, tspan, y0, 'StepSizeController', matlode.stepsizecontroller.Fixed(1000), varargin{:}); + sol = integrate@matlode.rk.erk.ERK(obj, f, tspan, y0, 'StepSizeController', matlode.stepsizecontroller.StandardController, 'Dense', obj.DenseOut, varargin{:}); end end end diff --git a/src/+matlode/+rk/+erk/ERK.m b/src/+matlode/+rk/+erk/ERK.m index 745491d..0e3007f 100644 --- a/src/+matlode/+rk/+erk/ERK.m +++ b/src/+matlode/+rk/+erk/ERK.m @@ -1,99 +1,304 @@ classdef ERK < matlode.rk.RungeKutta %Explicit Runge Kutta class + properties (Constant) + PartitionMethod = false + end + properties (SetAccess = immutable) - a - b - bHat - c - stage - order - embeddedOrder - adaptable + A + B + BHat + C + E + Stage + Order + EmbeddedOrder + Adaptive + Datatype FSAL end methods - - - function obj = ERK(a, b, bHat, c, order, embeddedOrder) - - obj.a = a; - obj.b = b; - obj.c = c; - obj.bHat = bHat; - obj.embeddedOrder = embeddedOrder; - obj.order = order; - obj.stage = length(b); + function obj = ERK(a, b, bHat, c, e, order, embeddedOrder) + + obj.A = a; + obj.B = b; + obj.BHat = bHat; + obj.C = c; + obj.E = e; + obj.EmbeddedOrder = embeddedOrder; + obj.Order = order; + obj.Stage = length(b); obj.FSAL = all(a(end, :) == b) && all(a(1, :) == 0); - obj.adaptable = ~isempty(bHat); + obj.Adaptive = ~isempty(bHat); + obj.Datatype = class(a); end - function sol = integrate(obj, f, tspan, y0, varargin) + + end + + methods (Access = protected) + function opts = matlodeSets(obj, p, varargin) - opts = obj.matlodeSets(varargin{:}); + %ERK specific options + + opts = matlodeSets@matlode.rk.RungeKutta(obj, p, varargin{:}); + end + + function [t, y, stats] = timeLoop(obj, f, tspan, y0, opts, varargin) numVars = length(y0); - t = tspan(1); - k = zeros(numVars, obj.stage); - sol.y = y0; + multiTspan = length(tspan) > 2; + isDense = ~isempty(opts.Dense); + overstep = false; + startupErr = true; + + %Structure access optimizations + errNormFunc = opts.ErrNorm; + stepController = opts.StepSizeController; + newStepFunc = @(prevAccept, tCur, tspan, hHist, err, q, nSteps, nFailed) opts.StepSizeController.newStepSize(prevAccept, tCur, tspan, hHist, err, q, nSteps, nFailed); + + + k = zeros(numVars, obj.Stage); + y = []; + t = []; - sol.stats.nsteps = 0; - sol.stats.nfailed = 0; - sol.stats.nfevals = 0; + if multiTspan + y = zeros(numVars, length(tspan)); + t = zeros(1, length(tspan)); + else + y = zeros(numVars, opts.ChunkSize); + t = zeros(1, opts.ChunkSize); + end + + %inital values + t(1,1) = tspan(1); + tCur = tspan(1); + tNext = tCur; + ti = 2; + hmax = opts.MaxStep; + + y(:, 1) = y0; + yNext = y0; + + %temporary stats + nSteps = 1; + nFailed = 0; + nFevals = 1; + nSmallSteps = 0; + + tdir = sign(tspan(end) - tspan(1)); + + q = min(obj.Order, obj.EmbeddedOrder); + + %%%% + % Start of method specific intializing code + %%%% fsalStart = uint32(obj.FSAL) + 1; - fevalIterCounts = obj.stage - fsalStart + 1; + fevalIterCounts = double(obj.Stage - fsalStart + 1); + %%%% + % End of method specific intializing code + %%%% - [h, k(:,1)] = opts.StepSizeController.startingStep(f, tspan, y0, obj.order, opts); + %First step + %allocates memory for first step + [hHist, k(:,end)] = stepController.startingStep(f, tspan, y0, obj.Order, errNormFunc, opts.InitialStep); + hN = hHist(1); + err = zeros(1, length(hHist)); + accept = true; - while t < tspan(end) + %Time Loop + while ti <= length(tspan) + %Cycle history + err = circshift(err, 1, 2); + hHist = circshift(hHist, 1, 2); - for i = fsalStart:obj.stage - k(:, i) = f(t + h * obj.c(i), sol.y + h * (k(:, 1:i-1) * obj.a(i, 1:i-1)')); - end + yCur = yNext; + tCur = tNext; + hC = hN; + hHist(1) = hN; - sol.stats.nfevals = sol.stats.nfevals + fevalIterCounts; + %Subject to change + hmin = 64 * eps(tCur); - yCurrent = sol.y + h * (k * obj.b'); - err = 0; + %%%% + % Start of method specific before time loop code + %%%% - - if opts.StepSizeController.adaptive - yEmbbeded = sol.y + h * (k * obj.bHat); - err = opts.ErrNorm(yCurrent, yEmbbeded, opts); + if obj.FSAL + k(:, 1) = k(:, end); end - [accept, h, t] = opts.StepSizeController.newStepSize(t, tspan, h, err, obj.embeddedOrder, sol.stats, opts); + %%%% + % End of method specific before time loop code + %%%% - if accept - sol.y = yCurrent; - sol.stats.nsteps = sol.stats.nsteps + 1; + %Accept Loop + %Will keep looping until accepted step + while true + + + %%%% + % Start of method specific time loop code + %%%% + + for i = fsalStart:obj.Stage + k(:, i) = f(tCur + hC * obj.C(i), yCur + k(:, 1:i-1) * (hC * obj.A(i, 1:i-1)')); + end + + yNext = yCur + k * (hC * obj.B'); + + nFevals = nFevals + fevalIterCounts; + prevAccept = accept; + + %Perform error estimates + %Could maybe leave out of method specific zone? + if stepController.Adaptive + yEmbbeded = k * (hC * obj.E'); + err(:, 1) = errNormFunc([yNext, yCur], yEmbbeded); + + if startupErr + err(:, 2:end) = err(:, 1); + end + startupErr = false; + end + + %%%% + % End of method specific time loop code + %%%% + + %Find next Step + [accept, hN, tNext] = newStepFunc(prevAccept, tCur, tspan, hHist, err, q, nSteps, nFailed); + + %set new step to be in range + if stepController.Adaptive + hN = min(hmax, hN); + end + + %Check if step is really small + if hN < hmin + if nSmallSteps == 0 + warning('The step the integrator is taking is extremely small, results may not be good') + end + %accept step since the step cannot get any smaller + nSmallSteps = nSmallSteps + 1; + hN = hmin; + tNext = tCur + hmin; + break; + end - if obj.FSAL - k(:, 1) = k(:, end); + %check step acception + if accept + break; end - else - sol.stats.nfailed = sol.stats.nfailed + 1; + nFailed = nFailed + 1; + hC = hN; + hHist(1) = hN; + + end + + nSteps = nSteps + 1; + + %check if controller failed to overstep for dense + if isDense && overstep && tNext < tspan(ti) + overstep = false; + end + + %Check for all memory allocation + if ~multiTspan + + %Allocate more memory if non-dense + if nSteps > length(t) + y = [y, zeros(numVars, opts.ChunkSize)]; + t = [t, zeros(1, opts.ChunkSize)]; + end + + t(:, nSteps) = tNext; + y(:, nSteps) = yNext; + end + + %check if condition holds + %used for end checking,integrate to, and dense output + if ~(tspan(ti) * tdir > (tNext + hN) * tdir) + %Dense Output + if isDense && multiTspan && ti ~= length(tspan) + + if overstep + while tNext > tspan(ti) && ti ~= length(tspan) + %Call Dense output function and record + + %y0 = yCur y1 = yNext + t(:, ti) = tspan(ti); + [y(:, ti), tempFevals] = opts.Dense.denseOut(f, tCur, tspan(ti), yCur, yNext, k, hC); + nFevals = nFevals + tempFevals; + + ti = ti + 1; + end + + if tspan(ti) - tNext < hmin + t(:, ti) = tspan(ti); + y(:, ti) = yNext; + ti = ti + 1; + end + + overstep = false; + else + overstep = true; + + %check for tspan(end) to hit exactly + if stepController.Adaptive && ~(tspan(end) * tdir > (tNext + hN) * tdir) + hN = tspan(end) - tNext; + end + end + else + %integrate to/ End point + %check if close enough with hmin + if tspan(ti) - tNext < hmin + + if multiTspan + t(:, ti) = tspan(ti); + y(:, ti) = yNext; + end + ti = ti + 1; + else + if stepController.Adaptive + hN = tspan(ti) - tNext; + end + + end + end end end + %End of Timeloop work + %Truncate extra memory + if ~multiTspan + t = t(:, 1:nSteps); + y = y(:, 1:nSteps); + else + t(:, end) = tspan(end); + y(:, end) = yNext; + end - end - - function opts = matlodeSets(obj, varargin) - p = inputParser; + %Create Stats + stats.nSteps = nSteps - 1; + stats.nFailed = nFailed; + stats.nFevals = nFevals; + stats.nSmallSteps = nSmallSteps; - %ERK specific options - opts = obj.matlodeSetsHelper(p, varargin{:}); end + + end + end diff --git a/src/+matlode/+rk/+erk/ForwardEuler.m b/src/+matlode/+rk/+erk/ForwardEuler.m index 1218a2d..5038c6d 100644 --- a/src/+matlode/+rk/+erk/ForwardEuler.m +++ b/src/+matlode/+rk/+erk/ForwardEuler.m @@ -1,26 +1,38 @@ classdef ForwardEuler < matlode.rk.erk.ERK methods - function obj = ForwardEuler() + function obj = ForwardEuler(datatype) %TODO %include reference - a = [[0]]; + %p = 1 & hatp = 0 & s = 1 + + if nargin == 0 + datatype = 'double'; + end + + caster = @(x) matlode.util.CoefficentTransformers.transform(x,datatype); + + a = caster('[[0]]'); - b = [1]; + b = caster('[1]'); bHat = []; - c = [0]; + e = []; + + c = caster('[0]'); order = 1; embbededOrder = 0; - obj = obj@matlode.rk.erk.ERK(a, b, bHat, c, order, embbededOrder); + obj = obj@matlode.rk.erk.ERK(a, b, bHat, c, e, order, embbededOrder); + + %p* = 1 & s* = 1 end function sol = integrate(obj, f, tspan, y0, varargin) - sol = integrate@matlode.rk.erk.ERK(obj, f, tspan, y0, 'StepSizeController', matlode.stepsizecontroller.Fixed(1000), varargin{:}); + sol = integrate@matlode.rk.erk.ERK(obj, f, tspan, y0,'Dense', matlode.denseoutput.Linear(obj.B), 'StepSizeController', matlode.stepsizecontroller.Fixed(1000), varargin{:}); end end end diff --git a/src/+matlode/+rk/RungeKutta.m b/src/+matlode/+rk/RungeKutta.m index a08eee4..0fa4a8e 100644 --- a/src/+matlode/+rk/RungeKutta.m +++ b/src/+matlode/+rk/RungeKutta.m @@ -2,14 +2,19 @@ %RungeKutta integrator properties (Abstract, SetAccess = immutable) - bHat - embeddedOrder - order + A + B + C + BHat + E + Stage + EmbeddedOrder + Order FSAL end methods - function obj = fromCoeffs(a, b, bHat, c, order, embededOrder) + function obj = fromCoeffs(a, b, bHat, c, e, bTheta, order, embededOrder) if istril(a) @@ -19,7 +24,7 @@ return end - obj = ERK(a, b, bHat, c, order, embededOrder); + obj = ERK(a, b, bHat, c, e, bTheta, order, embededOrder); return end @@ -31,11 +36,11 @@ end methods (Access = protected) - function opts = matlodeSetsHelper(obj, p, varargin) + function opts = matlodeSets(obj, p, varargin) %RungeKutta sepcific options - opts = matlodeSetsHelper@matlode.OneStepIntegrator(obj, p, varargin{:}); + opts = matlodeSets@matlode.OneStepIntegrator(obj, p, varargin{:}); end end diff --git a/src/+matlode/+stepsizecontroller/+soberpresets/H321.m b/src/+matlode/+stepsizecontroller/+soberpresets/H321.m new file mode 100644 index 0000000..15a743f --- /dev/null +++ b/src/+matlode/+stepsizecontroller/+soberpresets/H321.m @@ -0,0 +1,17 @@ +classdef H321 < matlode.stepsizecontroller.SoberlandController + %TODO Reference NASA Paper and Soberland paper + %This controller works and provides a smooth solution + + methods + function obj = H321() + a = [-1/3, -1/18, 5/18]; + b = [5/6, 1/6]; + + AQfunc = @(q, A) (A / q); + + %This is stupid, it ignores the first parameter for some reason + obj = obj@matlode.stepsizecontroller.SoberlandController(0, 'A', a, 'B', b, 'AQfunc', AQfunc); + end + end +end + diff --git a/src/+matlode/+stepsizecontroller/Fixed.m b/src/+matlode/+stepsizecontroller/Fixed.m index e65b775..aecb339 100644 --- a/src/+matlode/+stepsizecontroller/Fixed.m +++ b/src/+matlode/+stepsizecontroller/Fixed.m @@ -2,31 +2,36 @@ %Fixed stepsize controller properties (Constant) - adaptive = false; + Adaptive = false; end properties (SetAccess = immutable, GetAccess = private) - steps + Steps + end + + properties (SetAccess = immutable) + History end methods function obj = Fixed(steps) - obj.steps = steps; + obj.Steps = steps; + obj.History = 1; end - function [h0, f0] = startingStep(obj, f, tspan, y0, ~, options) - if options.InitialStep ~= 0 + function [h0, f0] = startingStep(obj, f, tspan, y0, ~, ~, intialStep) + if intialStep ~= 0 warning('Ignoring initialStep option, Controller is Fixed'); end f0 = f(tspan(1), y0); - h0 = (tspan(2) - tspan(1)) / obj.steps; + h0 = (tspan(2) - tspan(1)) / obj.Steps; end - function [accept, hCur, tNew] = newStepSize(~, ~, tspan, hCur, ~, ~, stepdlx, ~) + function [accept, h, tNew] = newStepSize(~, ~, ~, tspan, h, ~, ~, nSteps, ~, ~) accept = true; - tNew = (stepdlx.nsteps + 1) * hCur + tspan(1); + tNew = (nSteps) * h + tspan(1); end end diff --git a/src/+matlode/+stepsizecontroller/Gustafson.m b/src/+matlode/+stepsizecontroller/Gustafson.m new file mode 100644 index 0000000..389f2fe --- /dev/null +++ b/src/+matlode/+stepsizecontroller/Gustafson.m @@ -0,0 +1,70 @@ +classdef Gustafson < matlode.stepsizecontroller.StartingController + %TODO Reference Gustaffson paper(1991) + %Slight derivation of the Gustafson paper on the coefficents + + properties (Constant) + Adaptive = true; + end + + properties (SetAccess = immutable) + Fac + FacMin + FacMax + Ki + Kp + A + History + end + + methods + function obj = Gustafson(varargin) + + obj.History = 2; + + p = inputParser; + + p.addParameter('Fac', 0.95); + p.addParameter('FacMin', 0.1); + p.addParameter('FacMax', 2); + p.addParameter('Ki', 0.3); + p.addParameter('Kp', -0.4); + p.addParameter('A', 1); + + p.parse(varargin{:}); + + opts = p.Results; + + obj.Fac = opts.Fac; + obj.FacMin = opts.FacMin; + obj.FacMax = opts.FacMax; + obj.Ki = opts.Ki; + obj.Kp = opts.Kp; + obj.A = opts.A; + + end + + function [accept, hNew, tNew] = newStepSize(obj, prevAccept, t, ~, h, err, q, ~, ~) + accept = err(1) <= 1; + + scalh = 1; + + if accept + tNew = t + h(1); + if ~prevAccept + %H_0220 Controller + scalh = (h(1) / h(2))^obj.A; + end + %PI Controller + scalerr = err(1)^(obj.Ki/q) * err(2)^(obj.Kp/q); + else + tNew = t; + %Standard Controller(I controller) + scalerr = err(1)^(-1/(q + 1)); + end + + hNew = h(1) * min(obj.FacMax, max(obj.FacMin, obj.Fac * scalerr * scalh)); + + end + end +end + diff --git a/src/+matlode/+stepsizecontroller/SoberlandController.m b/src/+matlode/+stepsizecontroller/SoberlandController.m new file mode 100644 index 0000000..24c1136 --- /dev/null +++ b/src/+matlode/+stepsizecontroller/SoberlandController.m @@ -0,0 +1,78 @@ +classdef SoberlandController < matlode.stepsizecontroller.StartingController + %SOBERLANDCONTROLLER + %A general controller for all history based non-conditional controllers + %can be expanded up to any amount of history gien proper coefficents + + properties (Constant) + Adaptive = true; + end + + properties (SetAccess = immutable) + Fac + FacMax + FacMin + A + B + AQfunc % Geneall @(q, A) (A / (1+q)) + History + end + + methods + function obj = SoberlandController(varargin) + + p = inputParser; + + p.addParameter('Fac', 0.95); + p.addParameter('FacMin', 0.1); + p.addParameter('FacMax', 2); + p.addParameter('B', [0]); + p.addParameter('A', [1]); + p.addRequired('AQfunc'); + + p.parse(varargin{:}); + + opts = p.Results; + + obj.Fac = opts.Fac; + obj.FacMin = opts.FacMin; + obj.FacMax = opts.FacMax; + + if isempty(opts.A) + error('A cannot be empty') + end + + obj.A = opts.A; + obj.B = opts.B; + + %This is need because input parser struggles with + %function_handles + if ~isa(opts.AQfunc, 'function_handle') + error('AQfunc must be a function handle') + end + + obj.AQfunc = opts.AQfunc; + obj.History = length(obj.A); + end + + function [accept, hNew, tNew] = newStepSize(obj, ~, t, ~, h, err, q, ~, ~) + + accept = mean(err) <= 1; + if accept + tNew = t + h(1); + + else + tNew = t; + end + + + facScal = prod(err(1:obj.History).^obj.AQfunc(q,obj.A)); + hScal = prod((h(1:obj.History - 1) ./ h(2:obj.History)).^obj.B); + prediction = obj.Fac * facScal * hScal; + + hNew = h(1) * min(obj.FacMax, max(obj.FacMin, prediction)); + + + end + end +end + diff --git a/src/+matlode/+stepsizecontroller/StandardController.m b/src/+matlode/+stepsizecontroller/StandardController.m new file mode 100644 index 0000000..9f1058b --- /dev/null +++ b/src/+matlode/+stepsizecontroller/StandardController.m @@ -0,0 +1,49 @@ +classdef StandardController < matlode.stepsizecontroller.StartingController + %Standard Error Controller as described in Solving ODES I book + %Seperate from Soberland to help with performance(TODO test if true) + + properties (Constant) + Adaptive = true; + end + + properties (SetAccess = immutable) + Fac + FacMin + FacMax + History + end + + methods + function obj = StandardController(varargin) + + obj.History = 1; + p = inputParser; + + p.addParameter('Fac', 0.95); + p.addParameter('FacMin', 0.1); + p.addParameter('FacMax', 2); + + p.parse(varargin{:}); + + opts = p.Results; + + obj.Fac = opts.Fac; + obj.FacMin = opts.FacMin; + obj.FacMax = opts.FacMax; + end + + function [accept, hNew, tNew] = newStepSize(obj, ~, t, ~, h, err, q, ~, ~) + accept = err <= 1; + + if accept + tNew = t + h; + else + tNew = t; + end + + hNew = h * min(obj.FacMax, max(obj.FacMin, obj.Fac * err^(-1 / (q + 1)))); + + end + end +end + diff --git a/src/+matlode/+stepsizecontroller/StartingController.m b/src/+matlode/+stepsizecontroller/StartingController.m new file mode 100644 index 0000000..502898e --- /dev/null +++ b/src/+matlode/+stepsizecontroller/StartingController.m @@ -0,0 +1,41 @@ +classdef (Abstract) StartingController < matlode.stepsizecontroller.StepSizeController + %Adaptive controller + % This class is mostly used to contian the starting procedure of all + % adaptive controllers and some properties + + methods + function [h0, f0] = startingStep(obj, f, tspan, y0, order, errFunc, intialStep) + if intialStep ~= 0 + f0 = f(tspan(1), y0); + h0 = intialStep; + else + f0 = f(tspan(1), y0); + d0 = errFunc([y0, y0], 0); + d1 = errFunc([f0, f0], 0); + if any([d0, d1] < 1e-5) + h0 = 1e-6; + else + h0 = 0.01 * d0 / d1; + end + + y1 = y0 + h0 * f0; + f1 = f(tspan(1) + h0, y1); + err0 = errFunc([f1 - f0, f1 - f0], 0); + d2 = err0 / h0; + dm = max(d1, d2); + + if dm <= 1e-15 + h1 = max(1e-6, h0 * 1e-3); + else + h1 = (0.01 / dm)^(1 / (order + 1)); + end + + h0 = min(100 * h0, h1); + end + + %allocate for length of history + h0 = ones(1, obj.History) * h0; + end + end +end + diff --git a/src/+matlode/+stepsizecontroller/StepSizeController.m b/src/+matlode/+stepsizecontroller/StepSizeController.m index b6f45d0..e29c0ca 100644 --- a/src/+matlode/+stepsizecontroller/StepSizeController.m +++ b/src/+matlode/+stepsizecontroller/StepSizeController.m @@ -2,13 +2,18 @@ %Template for a step size controller properties (Abstract, Constant) - adaptive + Adaptive + end + + properties (Abstract, SetAccess = immutable) + History end methods (Abstract) - - [h0, f0] = startingStep(obj, f, tspan, y0, order, options); - [accept, hNew, tNew] = newStepSize(obj, tCur, tspan, hCur, err, embeddedOrder, stepdlx, options); + %accept0, h0, and err0, respect the amount of history needed for a + %controller + [h0, f0] = startingStep(obj, f, tspan, y0, order, errFunc, intialStep); + [accept, hNew, tNew] = newStepSize(obj, prevAccept, t, tspan, h, err, q, nSteps, nFailed); end end diff --git a/src/+matlode/+util/+errnorm/standardNorm.m b/src/+matlode/+util/+errnorm/standardNorm.m deleted file mode 100644 index e711a8d..0000000 --- a/src/+matlode/+util/+errnorm/standardNorm.m +++ /dev/null @@ -1,9 +0,0 @@ -function err = standardNorm(y, yHat, options) - %procedure described in solving ODEs I Harrier - %cannot assume that user has signal processing toolbox - - sc = (options.AbsTol + max(abs(y), abs(yHat)) .* options.RelTol); - - err = sqrt(1/length(y)) * norm((y - yHat) ./ sc); -end - diff --git a/src/+matlode/+util/CoefficentTransformers.m b/src/+matlode/+util/CoefficentTransformers.m new file mode 100644 index 0000000..9b0846e --- /dev/null +++ b/src/+matlode/+util/CoefficentTransformers.m @@ -0,0 +1,15 @@ +classdef CoefficentTransformers + methods (Static) + + function valuetype = transform(str, datatype) + + if isa(datatype, 'sym') + valuetype = str2sym(str); + else + %Better than str2double as it can deal with arrays and sqrt + valuetype = cast(str2num(str), datatype); + end + end + end +end + diff --git a/src/+matlode/Integrator.m b/src/+matlode/Integrator.m index 5d341c6..142308b 100644 --- a/src/+matlode/Integrator.m +++ b/src/+matlode/Integrator.m @@ -2,51 +2,100 @@ %INTEGRATOR Summary of this class goes here % Detailed explanation goes here + properties (Abstract, Constant) + PartitionMethod + end + properties (Abstract, SetAccess = immutable) - c - adaptable + Adaptive + Datatype end - methods (Abstract) - - %integrate - sol = integrate(obj, func, tspan, y0, varargin) - opts = matlodeSets(obj, varargin) - + methods (Abstract, Access = protected) + [t, y, stats] = timeLoop(obj, f, tspan, y0, opts, varargin); + end + + methods + function sol = integrate(obj, f, tspan, y0, varargin) + + if(length(tspan) < 2) + error('Unsupported tspan entry'); + end + + if ~issorted(tspan) && ~issorted(tspan, 'descend') + error('tspan must be sorted in order') + end + + %TODO + %optimize options for when only options struct is given + + %Create options + p = inputParser; + opts = obj.matlodeSets(p, varargin{:}); + + if isempty(opts.Dense) && length(tspan) > 2 && isa(opts.StepSizeController, 'matlode.stepsizecontroller.Fixed') + error('IntegrateTo is not supported for fixed step size, please use DenseOutput instead') + end + + t = []; + y = []; + stats = []; + + %Compute + if isa(f, 'function_handle') && ~obj.PartitionMethod + + [t, y, stats] = obj.timeLoop(f, tspan, y0, opts); + + elseif isa(f, 'cell') && cellfun(@(x) isa(x, 'function_handle'), f) + + %determine size of partition + if length(f) == 1 && ~obj.PartitionMethod + [t, y, stats] = obj.timeLoop(f{1}, tspan, y0, opts); + elseif length(f) > 1 + %TODO + %paritioning setup + + else + error('Need at least one function provided') + end + + else + error('f must be a ''function_handle'' or a ''cell'''); + end + + sol.t = t; + sol.y = y; + sol.stats = stats; + end end methods (Access = protected) - function opts = matlodeSetsHelper(obj, p, varargin) + function opts = matlodeSets(obj, p, varargin) - p.addParameter('AbsTol', 1e-6); - p.addParameter('RelTol', 1e-6); + %Assign to controllers p.addParameter('InitialStep', 0); p.addParameter('StepSizeController', matlode.stepsizecontroller.Fixed(1000)); - p.addParameter('StepsN', 0); - p.addParameter('ErrNorm', @(y, yHat, options) matlode.util.errnorm.standardNorm(y, yHat, options)); - p.addParameter('NewtonTol', 1e-7); - p.addParameter('NewtonMaxItr', 25); + p.addParameter('ErrNorm', matlode.errnorm.HistNorm.errEstimate(sqrt(eps), sqrt(eps))); p.addParameter('Jacobian', []); + p.addParameter('ChunkSize', 1000); + p.addParameter('Dense', [] ); + p.addParameter('MaxStep', inf); p.parse(varargin{:}); opts = p.Results; - %shorthanded fixed steps - if opts.StepsN > 0 - opts.StepSizeController = matlode.stepsizecontroller.Fixed(opts.StepsN); - end - opts = rmfield(opts, 'StepsN'); - %make sure a method can use a adaptive controller - if ~obj.adaptable && opts.StepSizeController.adaptive + if ~obj.Adaptive && opts.StepSizeController.Adaptive warning('Current method does not have the ability to use a adaptive step size controller'); %overrider user choice??? end + end + end end diff --git a/src/+matlode/OneStepIntegrator.m b/src/+matlode/OneStepIntegrator.m index f6653e6..c45ce68 100644 --- a/src/+matlode/OneStepIntegrator.m +++ b/src/+matlode/OneStepIntegrator.m @@ -1,20 +1,13 @@ classdef (Abstract) OneStepIntegrator < matlode.Integrator %One Step Integrator template - properties (Abstract, SetAccess = immutable) - a - b - stage - end methods (Access = protected) - function opts = matlodeSetsHelper(obj, p, varargin) - - + function opts = matlodeSets(obj, p, varargin) %OneStepIntegrator Specific options - opts = matlodeSetsHelper@matlode.Integrator(obj, p, varargin{:}); + opts = matlodeSets@matlode.Integrator(obj, p, varargin{:}); end From d0a294f65dfef611cde80aa69aea59d7aa341b8b Mon Sep 17 00:00:00 2001 From: rjg18 Date: Tue, 13 Jul 2021 10:57:18 -0400 Subject: [PATCH 04/10] Added Generic Starting Step, Updates to Backwards, Removed Abstract Immutable Properties --- src/+matlode/+errnorm/HistNorm.m | 11 --- src/+matlode/+errnorm/InfNorm.m | 2 +- src/+matlode/+errnorm/StandardNorm.m | 2 +- src/+matlode/+rk/+erk/ERK.m | 82 ++++++++----------- src/+matlode/+rk/RungeKutta.m | 18 +++- src/+matlode/+startingstep/BookStarting.m | 43 ++++++++++ src/+matlode/+startingstep/SetStep.m | 31 +++++++ src/+matlode/+startingstep/StartingStep.m | 19 +++++ src/+matlode/+startingstep/WattsStarting.m | 70 ++++++++++++++++ .../+stepsizecontroller/+soberpresets/H321.m | 17 ---- .../+stepsizecontroller/+soderpresets/H321.m | 16 ++++ src/+matlode/+stepsizecontroller/Fixed.m | 6 +- src/+matlode/+stepsizecontroller/Gustafson.m | 23 ++---- ...landController.m => SoderlandController.m} | 46 +++++------ .../+stepsizecontroller/StandardController.m | 26 +----- .../+stepsizecontroller/StartingController.m | 66 ++++++++------- .../+stepsizecontroller/StepSizeController.m | 17 +++- src/+matlode/Integrator.m | 44 +++++----- src/+matlode/OneStepIntegrator.m | 4 + 19 files changed, 337 insertions(+), 206 deletions(-) delete mode 100644 src/+matlode/+errnorm/HistNorm.m create mode 100644 src/+matlode/+startingstep/BookStarting.m create mode 100644 src/+matlode/+startingstep/SetStep.m create mode 100644 src/+matlode/+startingstep/StartingStep.m create mode 100644 src/+matlode/+startingstep/WattsStarting.m delete mode 100644 src/+matlode/+stepsizecontroller/+soberpresets/H321.m create mode 100644 src/+matlode/+stepsizecontroller/+soderpresets/H321.m rename src/+matlode/+stepsizecontroller/{SoberlandController.m => SoderlandController.m} (62%) diff --git a/src/+matlode/+errnorm/HistNorm.m b/src/+matlode/+errnorm/HistNorm.m deleted file mode 100644 index 6c44923..0000000 --- a/src/+matlode/+errnorm/HistNorm.m +++ /dev/null @@ -1,11 +0,0 @@ -classdef HistNorm < matlode.errnorm.NormErr - methods (Static) - function errFunc = errEstimate(AbsTol, RelTol) - - scFunc = @(y) (AbsTol + max(abs(y(:, 1)), abs(y(:, 2))) * RelTol); - errFunc = @(y, yE) sqrt(1/size(y,1)) * norm(yE ./ scFunc(y)); - - end - end -end - diff --git a/src/+matlode/+errnorm/InfNorm.m b/src/+matlode/+errnorm/InfNorm.m index e83702a..a065cb2 100644 --- a/src/+matlode/+errnorm/InfNorm.m +++ b/src/+matlode/+errnorm/InfNorm.m @@ -2,7 +2,7 @@ methods (Static) function errFunc = errEstimate(AbsTol, RelTol) - scFunc = @(y) (AbsTol + abs(y(:, 1)) * RelTol); + scFunc = @(y) (AbsTol + max(abs(y(:, 1)), abs(y(:, 2))) * RelTol); errFunc = @(y, yE) norm(yE ./ scFunc(y), 'Inf'); end diff --git a/src/+matlode/+errnorm/StandardNorm.m b/src/+matlode/+errnorm/StandardNorm.m index c6eb454..b96befe 100644 --- a/src/+matlode/+errnorm/StandardNorm.m +++ b/src/+matlode/+errnorm/StandardNorm.m @@ -2,7 +2,7 @@ methods (Static) function errFunc = errEstimate(AbsTol, RelTol) - scFunc = @(y) (AbsTol + abs(y(:, 1)) * RelTol); + scFunc = @(y) (AbsTol + max(abs(y(:, 1)), abs(y(:, 2))) * RelTol); errFunc = @(y, yE) sqrt(1/size(y,1)) * norm(yE ./ scFunc(y)); end diff --git a/src/+matlode/+rk/+erk/ERK.m b/src/+matlode/+rk/+erk/ERK.m index 0e3007f..3409397 100644 --- a/src/+matlode/+rk/+erk/ERK.m +++ b/src/+matlode/+rk/+erk/ERK.m @@ -5,38 +5,10 @@ PartitionMethod = false end - properties (SetAccess = immutable) - A - B - BHat - C - E - Stage - Order - EmbeddedOrder - Adaptive - Datatype - FSAL - end - methods - function obj = ERK(a, b, bHat, c, e, order, embeddedOrder) - - obj.A = a; - obj.B = b; - obj.BHat = bHat; - obj.C = c; - obj.E = e; - obj.EmbeddedOrder = embeddedOrder; - obj.Order = order; - obj.Stage = length(b); - obj.FSAL = all(a(end, :) == b) && all(a(1, :) == 0); - obj.Adaptive = ~isempty(bHat); - obj.Datatype = class(a); + obj = obj@matlode.rk.RungeKutta(a, b, bHat, c, e, order, embeddedOrder); end - - end methods (Access = protected) @@ -78,7 +50,10 @@ tCur = tspan(1); tNext = tCur; ti = 2; - hmax = opts.MaxStep; + + tdir = sign(tspan(end) - tspan(1)); + hmax = opts.MaxStep * tdir; + hmin = 64 * eps(tspan(1)) * tdir; y(:, 1) = y0; yNext = y0; @@ -86,13 +61,20 @@ %temporary stats nSteps = 1; nFailed = 0; - nFevals = 1; nSmallSteps = 0; - tdir = sign(tspan(end) - tspan(1)); - q = min(obj.Order, obj.EmbeddedOrder); + %First step + %allocates memory for first step + [h0, f0, nFevals] = stepController.startingStep(f, tspan, y0, obj.Order, errNormFunc, hmin, hmax); + + hHist = ones(1, stepController.History) * h0; + hN = hHist(1); + err = zeros(1, length(hHist)); + + accept = true; + %%%% % Start of method specific intializing code %%%% @@ -100,18 +82,19 @@ fsalStart = uint32(obj.FSAL) + 1; fevalIterCounts = double(obj.Stage - fsalStart + 1); + if obj.FSAL + if nFevals == 0 + k(:, end) = f(tspan(1), y0); + nFevals = nFevals + 1; + else + k(:, end) = f0; + end + end + %%%% % End of method specific intializing code %%%% - %First step - %allocates memory for first step - [hHist, k(:,end)] = stepController.startingStep(f, tspan, y0, obj.Order, errNormFunc, opts.InitialStep); - hN = hHist(1); - err = zeros(1, length(hHist)); - - accept = true; - %Time Loop while ti <= length(tspan) %Cycle history @@ -124,7 +107,7 @@ hHist(1) = hN; %Subject to change - hmin = 64 * eps(tCur); + hmin = 64 * eps(tCur) * tdir; %%%% % Start of method specific before time loop code @@ -178,13 +161,13 @@ %set new step to be in range if stepController.Adaptive - hN = min(hmax, hN); + hN = min(hmax * tdir, hN * tdir) * tdir; end %Check if step is really small - if hN < hmin + if hN * tdir < hmin * tdir if nSmallSteps == 0 - warning('The step the integrator is taking is extremely small, results may not be good') + warning('The step the integrator is taking extremely small, results may not be optimum') end %accept step since the step cannot get any smaller nSmallSteps = nSmallSteps + 1; @@ -207,7 +190,7 @@ nSteps = nSteps + 1; %check if controller failed to overstep for dense - if isDense && overstep && tNext < tspan(ti) + if isDense && overstep && tNext * tdir < tspan(ti) * tdir overstep = false; end @@ -231,10 +214,9 @@ if isDense && multiTspan && ti ~= length(tspan) if overstep - while tNext > tspan(ti) && ti ~= length(tspan) + while tNext * tdir > tspan(ti) * tdir && ti ~= length(tspan) %Call Dense output function and record - %y0 = yCur y1 = yNext t(:, ti) = tspan(ti); [y(:, ti), tempFevals] = opts.Dense.denseOut(f, tCur, tspan(ti), yCur, yNext, k, hC); nFevals = nFevals + tempFevals; @@ -242,7 +224,7 @@ ti = ti + 1; end - if tspan(ti) - tNext < hmin + if (tspan(ti) - tNext) * tdir < hmin * tdir t(:, ti) = tspan(ti); y(:, ti) = yNext; ti = ti + 1; @@ -260,7 +242,7 @@ else %integrate to/ End point %check if close enough with hmin - if tspan(ti) - tNext < hmin + if (tspan(ti) - tNext) * tdir < hmin * tdir if multiTspan t(:, ti) = tspan(ti); diff --git a/src/+matlode/+rk/RungeKutta.m b/src/+matlode/+rk/RungeKutta.m index 0fa4a8e..3eb9cd5 100644 --- a/src/+matlode/+rk/RungeKutta.m +++ b/src/+matlode/+rk/RungeKutta.m @@ -1,7 +1,7 @@ classdef (Abstract) RungeKutta < matlode.OneStepIntegrator %RungeKutta integrator - properties (Abstract, SetAccess = immutable) + properties (SetAccess = immutable) A B C @@ -14,6 +14,22 @@ end methods + function obj = RungeKutta(a, b, bHat, c, e, order, embeddedOrder) + + obj = obj@matlode.OneStepIntegrator(~isempty(bHat), class(a)); + + obj.A = a; + obj.B = b; + obj.BHat = bHat; + obj.C = c; + obj.E = e; + obj.EmbeddedOrder = embeddedOrder; + obj.Order = order; + obj.Stage = size(b, 2); + obj.FSAL = all(a(end, :) == b) && all(a(1, :) == 0); + + end + function obj = fromCoeffs(a, b, bHat, c, e, bTheta, order, embededOrder) diff --git a/src/+matlode/+startingstep/BookStarting.m b/src/+matlode/+startingstep/BookStarting.m new file mode 100644 index 0000000..0d630af --- /dev/null +++ b/src/+matlode/+startingstep/BookStarting.m @@ -0,0 +1,43 @@ +classdef BookStarting < matlode.startingstep.StartingStep + %This starting step procedure is derived from the Solving ODEs I book + + + methods + function obj = BookStarting(varargin) + obj = obj@matlode.startingstep.StartingStep(varargin{:}); + end + + function [h0, f0, fevals] = startingStep(obj, f, tspan, y0, order, errFunc, minStep, maxStep) + fevals = 2; + tdir = sign(tspan(end) - tspan(1)); + + f0 = f(tspan(1), y0); + + d0 = errFunc([y0, y0], 0); + d1 = errFunc([f0, f0], 0); + + if any([d0, d1] < 1e-5) + h0 = 1e-6; + else + h0 = 0.01 * d0 / d1; + end + + y1 = y0 + h0 * f0 * tdir; + f1 = f(tspan(1) + h0 * tdir, y1); + err0 = errFunc([f1 - f0, f1 - f0], 0); + d2 = err0 / h0; + dm = max(d1, d2); + + if dm <= eps + h1 = max(1e-6, h0 * 1e-3); + else + h1 = (0.01 / dm)^(1 / (order + 1)); + end + + + h0 = min(100 * h0, h1); + h0 = max([min([h0, maxStep * tdir]), minStep * tdir]) * tdir; + end + end +end + diff --git a/src/+matlode/+startingstep/SetStep.m b/src/+matlode/+startingstep/SetStep.m new file mode 100644 index 0000000..73f1f92 --- /dev/null +++ b/src/+matlode/+startingstep/SetStep.m @@ -0,0 +1,31 @@ +classdef SetStep < matlode.startingstep.StartingStep + %SETSTEP Summary of this class goes here + % Detailed explanation goes here + + properties + StartingStep + end + + methods + function obj = SetStep(ssp, varargin) + obj = obj@matlode.startingstep.StartingStep(varargin{:}); + obj.StartingStep = ssp; + end + + function [h0, f0, fevals] = startingStep(obj, ~, tspan, ~, ~, ~, minStep, maxStep) + fevals = 0; + f0 = []; + if obj.StartingStep > maxStep + error('The starting step is greater then the max step') + end + if obj.StartingStep < minStep + error('The starting step cannot be so small') + end + if sign(obj.StartingStep) ~= sign(tspan(end) - tspan(1)) + error('The starting step cannot go in the opposite direction') + end + h0 = max([min([obj.StartingStep, maxStep]), minStep]); + end + end +end + diff --git a/src/+matlode/+startingstep/StartingStep.m b/src/+matlode/+startingstep/StartingStep.m new file mode 100644 index 0000000..81779c2 --- /dev/null +++ b/src/+matlode/+startingstep/StartingStep.m @@ -0,0 +1,19 @@ +classdef StartingStep + %This is an abstract setup for dealing with different kinds of starting + %step procedures + + methods + function obj = StartingStep(varargin) + p = inputParser; + p.KeepUnmatched = false; + + %error checking using the parser + p.parse(varargin{:}); + end + end + + methods (Abstract) + [h0, f0, fevals] = startingStep(obj, f, tspan, y0, order, errFunc, minStep, maxStep); + end +end + diff --git a/src/+matlode/+startingstep/WattsStarting.m b/src/+matlode/+startingstep/WattsStarting.m new file mode 100644 index 0000000..353ac4e --- /dev/null +++ b/src/+matlode/+startingstep/WattsStarting.m @@ -0,0 +1,70 @@ +classdef WattsStarting < matlode.startingstep.StartingStep + %Based of Watts starting procedure + %Watts, H. A. (1983). Starting step size for an ODE solver. Journal of Computational and Applied Mathematics + + properties (SetAccess = immutable) + InitTol + InitEpsMax + InitEpsMin + Fac + end + + methods + function obj = WattsStarting(varargin) + + p = inputParser; + p.KeepUnmatched = true; + + p.addParameter('InitTol', eps^(1/2)); + p.addParameter('InitEpsMax', eps^(1/2)); + p.addParameter('InitEpsMin', 100 * eps); + p.addParameter('Fac', 0.8); + + p.parse(varargin{:}); + + opts = p.Results; + varargout = p.Unmatched; + + obj = obj@matlode.startingstep.StartingStep(varargout); + + obj.InitTol = opts.InitTol; + obj.InitEpsMax = opts.InitEpsMax; + obj.InitEpsMin = opts.InitEpsMin; + obj.Fac = opts.Fac; + + end + + function [h0, f0, fevals] = startingStep(obj, f, tspan, y0, order, errFunc, minStep, maxStep) + fevals = 2; + tdir = sign(tspan(end) - tspan(1)); + + f0 = f(tspan(1), y0); + fnorm = errFunc([f0, f0], 0); + + hstep = abs(tspan(1) - tspan(end)); + + h0 = max([min([obj.InitEpsMax * abs(tspan(1)), hstep]), obj.InitEpsMin * abs(tspan(1))]); + + if h0 < eps + h0 = eps^obj.InitEpsMax * hstep; + end + + y1 = y0 + h0 * tdir * f0; + + ydd = (1/(h0)) * (f(tspan(1) + h0 * tdir, y1) - f0) ; + + fddnorm = errFunc([ydd, ydd], 0); + + if fddnorm > eps + h0 = obj.Fac * sqrt(2) * obj.InitTol^(1/(order + 1)) / sqrt(fddnorm); + elseif fnorm > eps + h0 = obj.Fac * obj.InitTol^(1/(order + 1)) / fnorm; + else + h0 = obj.Fac * obj.InitTol^(1/(order + 1)) * abs(tspan(2) - tspan(1)); + end + + h0 = max([min([h0, maxStep * tdir]), minStep * tdir]) * tdir; + end + end +end + diff --git a/src/+matlode/+stepsizecontroller/+soberpresets/H321.m b/src/+matlode/+stepsizecontroller/+soberpresets/H321.m deleted file mode 100644 index 15a743f..0000000 --- a/src/+matlode/+stepsizecontroller/+soberpresets/H321.m +++ /dev/null @@ -1,17 +0,0 @@ -classdef H321 < matlode.stepsizecontroller.SoberlandController - %TODO Reference NASA Paper and Soberland paper - %This controller works and provides a smooth solution - - methods - function obj = H321() - a = [-1/3, -1/18, 5/18]; - b = [5/6, 1/6]; - - AQfunc = @(q, A) (A / q); - - %This is stupid, it ignores the first parameter for some reason - obj = obj@matlode.stepsizecontroller.SoberlandController(0, 'A', a, 'B', b, 'AQfunc', AQfunc); - end - end -end - diff --git a/src/+matlode/+stepsizecontroller/+soderpresets/H321.m b/src/+matlode/+stepsizecontroller/+soderpresets/H321.m new file mode 100644 index 0000000..e5b3e16 --- /dev/null +++ b/src/+matlode/+stepsizecontroller/+soderpresets/H321.m @@ -0,0 +1,16 @@ +classdef H321 < matlode.stepsizecontroller.SoderlandController + %Kennedy, C. A., & Carpenter, M. H. (2016). + %Diagonally Implicit Runge-Kutta methods for ordinary differential equations, a review + + methods + function obj = H321(varargin) + a = [-1/3, -1/18, 5/18]; + b = [5/6, 1/6]; + + AQfunc = @(q, A) (A / q); + + obj = obj@matlode.stepsizecontroller.SoderlandController('A', a, 'B', b, 'AQfunc', AQfunc, varargin{:}); + end + end +end + diff --git a/src/+matlode/+stepsizecontroller/Fixed.m b/src/+matlode/+stepsizecontroller/Fixed.m index aecb339..35aa4da 100644 --- a/src/+matlode/+stepsizecontroller/Fixed.m +++ b/src/+matlode/+stepsizecontroller/Fixed.m @@ -9,15 +9,11 @@ Steps end - properties (SetAccess = immutable) - History - end - methods function obj = Fixed(steps) + obj = obj@matlode.stepsizecontroller.StepSizeController(1); obj.Steps = steps; - obj.History = 1; end function [h0, f0] = startingStep(obj, f, tspan, y0, ~, ~, intialStep) diff --git a/src/+matlode/+stepsizecontroller/Gustafson.m b/src/+matlode/+stepsizecontroller/Gustafson.m index 389f2fe..c0052ea 100644 --- a/src/+matlode/+stepsizecontroller/Gustafson.m +++ b/src/+matlode/+stepsizecontroller/Gustafson.m @@ -1,31 +1,20 @@ classdef Gustafson < matlode.stepsizecontroller.StartingController - %TODO Reference Gustaffson paper(1991) - %Slight derivation of the Gustafson paper on the coefficents + %Gustafsson, K. (1991). Control theoretic techniques for stepsize + % selection in explicit Runge-Kutta methods. ACM Transactions on Mathematical Software - properties (Constant) - Adaptive = true; - end properties (SetAccess = immutable) - Fac - FacMin - FacMax Ki Kp A - History end methods function obj = Gustafson(varargin) - obj.History = 2; - p = inputParser; + p.KeepUnmatched = true; - p.addParameter('Fac', 0.95); - p.addParameter('FacMin', 0.1); - p.addParameter('FacMax', 2); p.addParameter('Ki', 0.3); p.addParameter('Kp', -0.4); p.addParameter('A', 1); @@ -33,10 +22,10 @@ p.parse(varargin{:}); opts = p.Results; + varargout = p.Unmatched; + + obj = obj@matlode.stepsizecontroller.StartingController(2, varargout); - obj.Fac = opts.Fac; - obj.FacMin = opts.FacMin; - obj.FacMax = opts.FacMax; obj.Ki = opts.Ki; obj.Kp = opts.Kp; obj.A = opts.A; diff --git a/src/+matlode/+stepsizecontroller/SoberlandController.m b/src/+matlode/+stepsizecontroller/SoderlandController.m similarity index 62% rename from src/+matlode/+stepsizecontroller/SoberlandController.m rename to src/+matlode/+stepsizecontroller/SoderlandController.m index 24c1136..543692f 100644 --- a/src/+matlode/+stepsizecontroller/SoberlandController.m +++ b/src/+matlode/+stepsizecontroller/SoderlandController.m @@ -1,57 +1,49 @@ -classdef SoberlandController < matlode.stepsizecontroller.StartingController +classdef SoderlandController < matlode.stepsizecontroller.StartingController %SOBERLANDCONTROLLER %A general controller for all history based non-conditional controllers %can be expanded up to any amount of history gien proper coefficents - properties (Constant) - Adaptive = true; - end + %Söderlind, G. (2003). Digital filters in adaptive time-stepping. + % ACM Transactions on Mathematical Software properties (SetAccess = immutable) - Fac - FacMax - FacMin A B - AQfunc % Geneall @(q, A) (A / (1+q)) - History + AQfunc % Geneall @(q, A) (A / q) end methods - function obj = SoberlandController(varargin) + function obj = SoderlandController(varargin) - p = inputParser; + p = inputParser; + p.KeepUnmatched = true; - p.addParameter('Fac', 0.95); - p.addParameter('FacMin', 0.1); - p.addParameter('FacMax', 2); - p.addParameter('B', [0]); p.addParameter('A', [1]); - p.addRequired('AQfunc'); + p.addParameter('B', [0]); + p.addParameter('AQfunc', @(q, A) (A/q)); p.parse(varargin{:}); opts = p.Results; - - obj.Fac = opts.Fac; - obj.FacMin = opts.FacMin; - obj.FacMax = opts.FacMax; + varargout = p.Unmatched; if isempty(opts.A) - error('A cannot be empty') + error('A cannot be empty'); end - obj.A = opts.A; - obj.B = opts.B; + if isempty(opts.B) + error('B cannots be empty'); + end - %This is need because input parser struggles with - %function_handles if ~isa(opts.AQfunc, 'function_handle') error('AQfunc must be a function handle') end + obj = obj@matlode.stepsizecontroller.StartingController(max(length(opts.A), length(opts.B)), varargout); + + obj.A = opts.A; + obj.B = opts.B; obj.AQfunc = opts.AQfunc; - obj.History = length(obj.A); end function [accept, hNew, tNew] = newStepSize(obj, ~, t, ~, h, err, q, ~, ~) @@ -66,7 +58,7 @@ facScal = prod(err(1:obj.History).^obj.AQfunc(q,obj.A)); - hScal = prod((h(1:obj.History - 1) ./ h(2:obj.History)).^obj.B); + hScal = prod(abs(h(1:obj.History - 1) ./ h(2:obj.History)).^obj.B); prediction = obj.Fac * facScal * hScal; hNew = h(1) * min(obj.FacMax, max(obj.FacMin, prediction)); diff --git a/src/+matlode/+stepsizecontroller/StandardController.m b/src/+matlode/+stepsizecontroller/StandardController.m index 9f1058b..517a220 100644 --- a/src/+matlode/+stepsizecontroller/StandardController.m +++ b/src/+matlode/+stepsizecontroller/StandardController.m @@ -1,35 +1,13 @@ classdef StandardController < matlode.stepsizecontroller.StartingController %Standard Error Controller as described in Solving ODES I book - %Seperate from Soberland to help with performance(TODO test if true) + %Seperate from Soberland to help with performance - properties (Constant) - Adaptive = true; - end - - properties (SetAccess = immutable) - Fac - FacMin - FacMax - History - end methods function obj = StandardController(varargin) - obj.History = 1; - p = inputParser; - - p.addParameter('Fac', 0.95); - p.addParameter('FacMin', 0.1); - p.addParameter('FacMax', 2); - - p.parse(varargin{:}); - - opts = p.Results; + obj = obj@matlode.stepsizecontroller.StartingController(1, varargin{:}); - obj.Fac = opts.Fac; - obj.FacMin = opts.FacMin; - obj.FacMax = opts.FacMax; end function [accept, hNew, tNew] = newStepSize(obj, ~, t, ~, h, err, q, ~, ~) diff --git a/src/+matlode/+stepsizecontroller/StartingController.m b/src/+matlode/+stepsizecontroller/StartingController.m index 502898e..a44b612 100644 --- a/src/+matlode/+stepsizecontroller/StartingController.m +++ b/src/+matlode/+stepsizecontroller/StartingController.m @@ -3,38 +3,44 @@ % This class is mostly used to contian the starting procedure of all % adaptive controllers and some properties + properties (Constant) + Adaptive = true; + end + + properties (SetAccess = immutable) + Fac + FacMax + FacMin + InitalStep + end + + methods - function [h0, f0] = startingStep(obj, f, tspan, y0, order, errFunc, intialStep) - if intialStep ~= 0 - f0 = f(tspan(1), y0); - h0 = intialStep; - else - f0 = f(tspan(1), y0); - d0 = errFunc([y0, y0], 0); - d1 = errFunc([f0, f0], 0); - if any([d0, d1] < 1e-5) - h0 = 1e-6; - else - h0 = 0.01 * d0 / d1; - end - - y1 = y0 + h0 * f0; - f1 = f(tspan(1) + h0, y1); - err0 = errFunc([f1 - f0, f1 - f0], 0); - d2 = err0 / h0; - dm = max(d1, d2); - - if dm <= 1e-15 - h1 = max(1e-6, h0 * 1e-3); - else - h1 = (0.01 / dm)^(1 / (order + 1)); - end - - h0 = min(100 * h0, h1); - end + function obj = StartingController(hist, varargin) + p = inputParser; + p.KeepUnmatched = true; + + p.addParameter('Fac', 0.8); + p.addParameter('FacMin', 0.1); + p.addParameter('FacMax', 2); + p.addParameter('InitalStep', matlode.startingstep.BookStarting()); + + p.parse(varargin{:}); - %allocate for length of history - h0 = ones(1, obj.History) * h0; + opts = p.Results; + varargout = p.Unmatched; + + obj = obj@matlode.stepsizecontroller.StepSizeController(hist, varargout); + + obj.Fac = opts.Fac; + obj.FacMin = opts.FacMin; + obj.FacMax = opts.FacMax; + obj.InitalStep = opts.InitalStep; + + end + + function [h0, f0, fevals] = startingStep(obj, f, tspan, y0, order, errFunc, minStep, maxStep) + [h0, f0, fevals] = obj.InitalStep.startingStep(f, tspan, y0, order, errFunc, minStep, maxStep); end end end diff --git a/src/+matlode/+stepsizecontroller/StepSizeController.m b/src/+matlode/+stepsizecontroller/StepSizeController.m index e29c0ca..0677b9b 100644 --- a/src/+matlode/+stepsizecontroller/StepSizeController.m +++ b/src/+matlode/+stepsizecontroller/StepSizeController.m @@ -5,14 +5,27 @@ Adaptive end - properties (Abstract, SetAccess = immutable) + properties (SetAccess = immutable) History end + methods (Access = protected) + function obj = StepSizeController(hist, varargin) + + p = inputParser; + p.KeepUnmatched = false; + + %Gives better error checking + p.parse(varargin{:}); + + obj.History = hist; + end + end + methods (Abstract) %accept0, h0, and err0, respect the amount of history needed for a %controller - [h0, f0] = startingStep(obj, f, tspan, y0, order, errFunc, intialStep); + [h0, f0, fevals] = startingStep(obj, f, tspan, y0, order, errFunc, minStep, maxStep); [accept, hNew, tNew] = newStepSize(obj, prevAccept, t, tspan, h, err, q, nSteps, nFailed); end diff --git a/src/+matlode/Integrator.m b/src/+matlode/Integrator.m index 142308b..2a4c36a 100644 --- a/src/+matlode/Integrator.m +++ b/src/+matlode/Integrator.m @@ -6,7 +6,7 @@ PartitionMethod end - properties (Abstract, SetAccess = immutable) + properties (SetAccess = immutable) Adaptive Datatype end @@ -16,25 +16,32 @@ end methods + function obj = Integrator(adap, datatype) + + obj.Adaptive = adap; + obj.Datatype = datatype; + end + function sol = integrate(obj, f, tspan, y0, varargin) - if(length(tspan) < 2) - error('Unsupported tspan entry'); - end + [n,m] = size(tspan); - if ~issorted(tspan) && ~issorted(tspan, 'descend') - error('tspan must be sorted in order') + if ~ismatrix(tspan) + error('tspan cannot have this many dimensions'); + elseif n ~= 1 && m ~= 1 + error('tspan must be a vector') + elseif (length(tspan) < 2) + error('tspan must have a initial and final entry'); + elseif ~issorted(tspan) && ~issorted(tspan, 'descend') + error('tspan must be sorted in either descending or ascending order') end - %TODO - %optimize options for when only options struct is given - %Create options p = inputParser; opts = obj.matlodeSets(p, varargin{:}); - if isempty(opts.Dense) && length(tspan) > 2 && isa(opts.StepSizeController, 'matlode.stepsizecontroller.Fixed') - error('IntegrateTo is not supported for fixed step size, please use DenseOutput instead') + if isempty(opts.Dense) && length(tspan) > 2 && ~opts.StepSizeController.Adaptive + error('IntegrateTo is not supported for a non-adaptable controller, please use DenseOutput instead') end t = []; @@ -51,12 +58,11 @@ %determine size of partition if length(f) == 1 && ~obj.PartitionMethod [t, y, stats] = obj.timeLoop(f{1}, tspan, y0, opts); - elseif length(f) > 1 + elseif isempty(f) + error('Partitioned system solves, require atleast one partition to be provided') + else %TODO %paritioning setup - - else - error('Need at least one function provided') end else @@ -75,12 +81,11 @@ function opts = matlodeSets(obj, p, varargin) %Assign to controllers - p.addParameter('InitialStep', 0); p.addParameter('StepSizeController', matlode.stepsizecontroller.Fixed(1000)); - p.addParameter('ErrNorm', matlode.errnorm.HistNorm.errEstimate(sqrt(eps), sqrt(eps))); + p.addParameter('ErrNorm', matlode.errnorm.InfNorm.errEstimate(sqrt(eps), sqrt(eps))); p.addParameter('Jacobian', []); p.addParameter('ChunkSize', 1000); - p.addParameter('Dense', [] ); + p.addParameter('Dense', []); p.addParameter('MaxStep', inf); p.parse(varargin{:}); @@ -89,8 +94,7 @@ %make sure a method can use a adaptive controller if ~obj.Adaptive && opts.StepSizeController.Adaptive - warning('Current method does not have the ability to use a adaptive step size controller'); - %overrider user choice??? + error('Current method does not have the ability to use a adaptive step size controller'); end diff --git a/src/+matlode/OneStepIntegrator.m b/src/+matlode/OneStepIntegrator.m index c45ce68..65ded7e 100644 --- a/src/+matlode/OneStepIntegrator.m +++ b/src/+matlode/OneStepIntegrator.m @@ -3,6 +3,10 @@ methods (Access = protected) + function obj = OneStepIntegrator(varargin) + obj = obj@matlode.Integrator(varargin{:}); + end + function opts = matlodeSets(obj, p, varargin) %OneStepIntegrator Specific options From 5aaa799933ea642cc5ab1de329f8fbe1b180ef66 Mon Sep 17 00:00:00 2001 From: rjg18 Date: Wed, 14 Jul 2021 13:06:06 -0400 Subject: [PATCH 05/10] Simplifications to the time loop, Validations, small changes --- src/+matlode/+rk/+erk/ERK.m | 115 +++++++----------- src/+matlode/+startingstep/BookStarting.m | 4 +- src/+matlode/+startingstep/SetStep.m | 12 +- src/+matlode/+startingstep/WattsStarting.m | 10 +- src/+matlode/+stepsizecontroller/Gustafson.m | 6 +- .../+stepsizecontroller/StartingController.m | 14 +-- src/+matlode/+util/scalarValidationFunc.m | 5 + src/+matlode/Integrator.m | 11 +- 8 files changed, 83 insertions(+), 94 deletions(-) create mode 100644 src/+matlode/+util/scalarValidationFunc.m diff --git a/src/+matlode/+rk/+erk/ERK.m b/src/+matlode/+rk/+erk/ERK.m index 3409397..69d5fc9 100644 --- a/src/+matlode/+rk/+erk/ERK.m +++ b/src/+matlode/+rk/+erk/ERK.m @@ -24,8 +24,6 @@ numVars = length(y0); multiTspan = length(tspan) > 2; isDense = ~isempty(opts.Dense); - overstep = false; - startupErr = true; %Structure access optimizations errNormFunc = opts.ErrNorm; @@ -33,6 +31,7 @@ newStepFunc = @(prevAccept, tCur, tspan, hHist, err, q, nSteps, nFailed) opts.StepSizeController.newStepSize(prevAccept, tCur, tspan, hHist, err, q, nSteps, nFailed); + k = zeros(numVars, obj.Stage); y = []; t = []; @@ -50,9 +49,15 @@ tCur = tspan(1); tNext = tCur; ti = 2; + tl = 2; + if isDense + ti = length(tspan); + else + tl = length(tspan); + end tdir = sign(tspan(end) - tspan(1)); - hmax = opts.MaxStep * tdir; + hmax = min([abs(opts.MaxStep), abs(tspan(end) - tspan(1))]) * tdir; hmin = 64 * eps(tspan(1)) * tdir; y(:, 1) = y0; @@ -71,7 +76,7 @@ hHist = ones(1, stepController.History) * h0; hN = hHist(1); - err = zeros(1, length(hHist)); + err = ones(1, length(hHist)); accept = true; @@ -125,7 +130,6 @@ %Will keep looping until accepted step while true - %%%% % Start of method specific time loop code %%%% @@ -145,11 +149,6 @@ if stepController.Adaptive yEmbbeded = k * (hC * obj.E'); err(:, 1) = errNormFunc([yNext, yCur], yEmbbeded); - - if startupErr - err(:, 2:end) = err(:, 1); - end - startupErr = false; end %%%% @@ -159,15 +158,10 @@ %Find next Step [accept, hN, tNext] = newStepFunc(prevAccept, tCur, tspan, hHist, err, q, nSteps, nFailed); - %set new step to be in range - if stepController.Adaptive - hN = min(hmax * tdir, hN * tdir) * tdir; - end - %Check if step is really small - if hN * tdir < hmin * tdir + if abs(hN) < abs(hmin) if nSmallSteps == 0 - warning('The step the integrator is taking extremely small, results may not be optimum') + warning('The step the integrator is taking extremely small, results may not be optimal') end %accept step since the step cannot get any smaller nSmallSteps = nSmallSteps + 1; @@ -187,13 +181,13 @@ end - nSteps = nSteps + 1; - - %check if controller failed to overstep for dense - if isDense && overstep && tNext * tdir < tspan(ti) * tdir - overstep = false; + %set new step to be in range + if stepController.Adaptive + hN = min(abs(hmax), abs(hN)) * tdir; end + nSteps = nSteps + 1; + %Check for all memory allocation if ~multiTspan @@ -209,53 +203,41 @@ %check if condition holds %used for end checking,integrate to, and dense output - if ~(tspan(ti) * tdir > (tNext + hN) * tdir) - %Dense Output - if isDense && multiTspan && ti ~= length(tspan) - - if overstep - while tNext * tdir > tspan(ti) * tdir && ti ~= length(tspan) - %Call Dense output function and record - - t(:, ti) = tspan(ti); - [y(:, ti), tempFevals] = opts.Dense.denseOut(f, tCur, tspan(ti), yCur, yNext, k, hC); - nFevals = nFevals + tempFevals; - - ti = ti + 1; - end - - if (tspan(ti) - tNext) * tdir < hmin * tdir - t(:, ti) = tspan(ti); - y(:, ti) = yNext; - ti = ti + 1; - end + if ~(tspan(ti) * tdir > (tNext + hN) * tdir) && tl == length(tspan) + + %integrate to/ End point + %check if close enough with hmin + if (tspan(ti) - tNext) * tdir < hmin * tdir - overstep = false; - else - overstep = true; - - %check for tspan(end) to hit exactly - if stepController.Adaptive && ~(tspan(end) * tdir > (tNext + hN) * tdir) - hN = tspan(end) - tNext; - end + if multiTspan + t(:, ti) = tspan(ti); + y(:, ti) = yNext; end + ti = ti + 1; else - %integrate to/ End point - %check if close enough with hmin - if (tspan(ti) - tNext) * tdir < hmin * tdir - - if multiTspan - t(:, ti) = tspan(ti); - y(:, ti) = yNext; - end - ti = ti + 1; - else - if stepController.Adaptive - hN = tspan(ti) - tNext; - end - + + if stepController.Adaptive + hN = tspan(ti) - tNext; end end + + elseif isDense && multiTspan && tl ~= length(tspan) && ~(tspan(tl) * tdir > (tNext + hN) * tdir) + + while tNext * tdir > tspan(tl) * tdir && tl ~= length(tspan) + %Call Dense output function and record + + t(:, tl) = tspan(tl); + [y(:, tl), tempFevals] = opts.Dense.denseOut(f, tCur, tspan(tl), yCur, yNext, k, hC); + nFevals = nFevals + tempFevals; + + tl = tl + 1; + end + + if tl == length(tspan) && (tspan(ti) - tNext) * tdir < hmin * tdir + t(:, ti) = tspan(ti); + y(:, ti) = yNext; + ti = ti + 1; + end end end @@ -275,12 +257,7 @@ stats.nFevals = nFevals; stats.nSmallSteps = nSmallSteps; - end - - end - - end diff --git a/src/+matlode/+startingstep/BookStarting.m b/src/+matlode/+startingstep/BookStarting.m index 0d630af..4afd5bc 100644 --- a/src/+matlode/+startingstep/BookStarting.m +++ b/src/+matlode/+startingstep/BookStarting.m @@ -7,7 +7,7 @@ obj = obj@matlode.startingstep.StartingStep(varargin{:}); end - function [h0, f0, fevals] = startingStep(obj, f, tspan, y0, order, errFunc, minStep, maxStep) + function [h0, f0, fevals] = startingStep(~, f, tspan, y0, order, errFunc, minStep, maxStep) fevals = 2; tdir = sign(tspan(end) - tspan(1)); @@ -21,6 +21,8 @@ else h0 = 0.01 * d0 / d1; end + + h0 = max([min([h0, maxStep * tdir]), minStep * tdir]); y1 = y0 + h0 * f0 * tdir; f1 = f(tspan(1) + h0 * tdir, y1); diff --git a/src/+matlode/+startingstep/SetStep.m b/src/+matlode/+startingstep/SetStep.m index 73f1f92..2373fd9 100644 --- a/src/+matlode/+startingstep/SetStep.m +++ b/src/+matlode/+startingstep/SetStep.m @@ -15,15 +15,15 @@ function [h0, f0, fevals] = startingStep(obj, ~, tspan, ~, ~, ~, minStep, maxStep) fevals = 0; f0 = []; - if obj.StartingStep > maxStep - error('The starting step is greater then the max step') - end - if obj.StartingStep < minStep - error('The starting step cannot be so small') - end if sign(obj.StartingStep) ~= sign(tspan(end) - tspan(1)) error('The starting step cannot go in the opposite direction') end + if obj.StartingStep > maxStep + error('The starting step is greater then the maximum step %f', maxStep) + elseif obj.StartingStep < minStep + error('The starting step cannot be smaller than the minimum step %f', minStep) + end + h0 = max([min([obj.StartingStep, maxStep]), minStep]); end end diff --git a/src/+matlode/+startingstep/WattsStarting.m b/src/+matlode/+startingstep/WattsStarting.m index 353ac4e..39f117d 100644 --- a/src/+matlode/+startingstep/WattsStarting.m +++ b/src/+matlode/+startingstep/WattsStarting.m @@ -1,5 +1,5 @@ classdef WattsStarting < matlode.startingstep.StartingStep - %Based of Watts starting procedure + %Based of Watt's starting procedure %Watts, H. A. (1983). Starting step size for an ODE solver. Journal of Computational and Applied Mathematics properties (SetAccess = immutable) @@ -15,10 +15,10 @@ p = inputParser; p.KeepUnmatched = true; - p.addParameter('InitTol', eps^(1/2)); - p.addParameter('InitEpsMax', eps^(1/2)); - p.addParameter('InitEpsMin', 100 * eps); - p.addParameter('Fac', 0.8); + p.addParameter('InitTol', eps^(1/2), matlode.util.scalarValidationFunc); + p.addParameter('InitEpsMax', eps^(1/2), matlode.util.scalarValidationFunc); + p.addParameter('InitEpsMin', 100 * eps, matlode.util.scalarValidationFunc); + p.addParameter('Fac', 0.8, matlode.util.scalarValidationFunc); p.parse(varargin{:}); diff --git a/src/+matlode/+stepsizecontroller/Gustafson.m b/src/+matlode/+stepsizecontroller/Gustafson.m index c0052ea..f13d01e 100644 --- a/src/+matlode/+stepsizecontroller/Gustafson.m +++ b/src/+matlode/+stepsizecontroller/Gustafson.m @@ -15,9 +15,9 @@ p = inputParser; p.KeepUnmatched = true; - p.addParameter('Ki', 0.3); - p.addParameter('Kp', -0.4); - p.addParameter('A', 1); + p.addParameter('Ki', 0.3, matlode.util.scalarValidationFunc); + p.addParameter('Kp', -0.4, matlode.util.scalarValidationFunc); + p.addParameter('A', 1, matlode.util.scalarValidationFunc); p.parse(varargin{:}); diff --git a/src/+matlode/+stepsizecontroller/StartingController.m b/src/+matlode/+stepsizecontroller/StartingController.m index a44b612..87f6093 100644 --- a/src/+matlode/+stepsizecontroller/StartingController.m +++ b/src/+matlode/+stepsizecontroller/StartingController.m @@ -1,7 +1,7 @@ classdef (Abstract) StartingController < matlode.stepsizecontroller.StepSizeController - %Adaptive controller - % This class is mostly used to contian the starting procedure of all - % adaptive controllers and some properties + % This class is used to contain common adaptive stepsizecontroler + % properties and to allow the user to change the starting step + % procedure properties (Constant) Adaptive = true; @@ -20,10 +20,10 @@ p = inputParser; p.KeepUnmatched = true; - p.addParameter('Fac', 0.8); - p.addParameter('FacMin', 0.1); - p.addParameter('FacMax', 2); - p.addParameter('InitalStep', matlode.startingstep.BookStarting()); + p.addParameter('Fac', 0.8, matlode.util.scalarValidationFunc); + p.addParameter('FacMin', 0.1, matlode.util.scalarValidationFunc); + p.addParameter('FacMax', 2, matlode.util.scalarValidationFunc); + p.addParameter('InitalStep', matlode.startingstep.WattsStarting()); p.parse(varargin{:}); diff --git a/src/+matlode/+util/scalarValidationFunc.m b/src/+matlode/+util/scalarValidationFunc.m new file mode 100644 index 0000000..8ef24b8 --- /dev/null +++ b/src/+matlode/+util/scalarValidationFunc.m @@ -0,0 +1,5 @@ +function validFunc = scalarValidationFunc() +msg = 'Value must be a scalar'; +validFunc = @(x) assert(isscalar(x), msg); +end + diff --git a/src/+matlode/Integrator.m b/src/+matlode/Integrator.m index 2a4c36a..e4ae86e 100644 --- a/src/+matlode/Integrator.m +++ b/src/+matlode/Integrator.m @@ -32,8 +32,8 @@ error('tspan must be a vector') elseif (length(tspan) < 2) error('tspan must have a initial and final entry'); - elseif ~issorted(tspan) && ~issorted(tspan, 'descend') - error('tspan must be sorted in either descending or ascending order') + elseif ~issorted(tspan, 'strictascend') && ~issorted(tspan, 'strictdescend') + error('tspan must have unique values and be sorted in either descending or ascending order') end %Create options @@ -44,6 +44,7 @@ error('IntegrateTo is not supported for a non-adaptable controller, please use DenseOutput instead') end + t = []; y = []; stats = []; @@ -61,8 +62,12 @@ elseif isempty(f) error('Partitioned system solves, require atleast one partition to be provided') else - %TODO %paritioning setup + if obj.PartitionMethod + [t, y, stats] = obj.timeLoop(f, tspan, y0, opts); + else + error('This method does not work with a partitioned f'); + end end else From c7f70cab9902c6d328319e8867741356db7d78bd Mon Sep 17 00:00:00 2001 From: rjg18 Date: Thu, 15 Jul 2021 13:06:55 -0400 Subject: [PATCH 06/10] Quick Fixes, and Time Loop Performance --- src/+matlode/+rk/+erk/ERK.m | 73 +++++++++++-------- src/+matlode/+startingstep/SetStep.m | 2 +- src/+matlode/+stepsizecontroller/Fixed.m | 9 +-- .../+stepsizecontroller/SoderlandController.m | 2 +- 4 files changed, 46 insertions(+), 40 deletions(-) diff --git a/src/+matlode/+rk/+erk/ERK.m b/src/+matlode/+rk/+erk/ERK.m index 69d5fc9..48e24ea 100644 --- a/src/+matlode/+rk/+erk/ERK.m +++ b/src/+matlode/+rk/+erk/ERK.m @@ -36,14 +36,19 @@ y = []; t = []; + tlen = 0; + if multiTspan y = zeros(numVars, length(tspan)); t = zeros(1, length(tspan)); else y = zeros(numVars, opts.ChunkSize); t = zeros(1, opts.ChunkSize); + tlen = length(t); end + + %inital values t(1,1) = tspan(1); tCur = tspan(1); @@ -56,6 +61,8 @@ tl = length(tspan); end + tspanlen = length(tspan); + tdir = sign(tspan(end) - tspan(1)); hmax = min([abs(opts.MaxStep), abs(tspan(end) - tspan(1))]) * tdir; hmin = 64 * eps(tspan(1)) * tdir; @@ -101,7 +108,7 @@ %%%% %Time Loop - while ti <= length(tspan) + while ti <= tspanlen %Cycle history err = circshift(err, 1, 2); hHist = circshift(hHist, 1, 2); @@ -137,17 +144,14 @@ for i = fsalStart:obj.Stage k(:, i) = f(tCur + hC * obj.C(i), yCur + k(:, 1:i-1) * (hC * obj.A(i, 1:i-1)')); end - - yNext = yCur + k * (hC * obj.B'); - - nFevals = nFevals + fevalIterCounts; - prevAccept = accept; + khc = k * hC; + yNext = yCur + khc * obj.B'; %Perform error estimates %Could maybe leave out of method specific zone? if stepController.Adaptive - yEmbbeded = k * (hC * obj.E'); + yEmbbeded = khc * obj.E'; err(:, 1) = errNormFunc([yNext, yCur], yEmbbeded); end @@ -155,6 +159,10 @@ % End of method specific time loop code %%%% + nFevals = nFevals + fevalIterCounts; + + prevAccept = accept; + %Find next Step [accept, hN, tNext] = newStepFunc(prevAccept, tCur, tspan, hHist, err, q, nSteps, nFailed); @@ -192,53 +200,54 @@ if ~multiTspan %Allocate more memory if non-dense - if nSteps > length(t) + if nSteps > tlen y = [y, zeros(numVars, opts.ChunkSize)]; t = [t, zeros(1, opts.ChunkSize)]; + tlen = length(t); end t(:, nSteps) = tNext; y(:, nSteps) = yNext; + elseif isDense + %Dense output + while tl ~= tspanlen && tNext * tdir > tspan(tl) * tdir + + %incase lucky and land on point, can be cheaper + %than dense output if dense requires fevals + if abs(tspan(tl) - tNext) < abs(hmin) + t(:, tl) = tspan(tl); + y(:, tl) = yNext; + else + %Call Dense output function and record + t(:, tl) = tspan(tl); + [y(:, tl), tempFevals] = opts.Dense.denseOut(f, tCur, tspan(tl), yCur, yNext, k, hC); + nFevals = nFevals + tempFevals; + end + tl = tl + 1; + end end %check if condition holds - %used for end checking,integrate to, and dense output - if ~(tspan(ti) * tdir > (tNext + hN) * tdir) && tl == length(tspan) + %used for end checking and integrate to + if tspan(ti) * tdir <= (tNext + hN) * tdir %integrate to/ End point %check if close enough with hmin - if (tspan(ti) - tNext) * tdir < hmin * tdir + if abs(tspan(ti) - tNext) < abs(hmin) if multiTspan t(:, ti) = tspan(ti); y(:, ti) = yNext; end ti = ti + 1; - else + elseif stepController.Adaptive - if stepController.Adaptive - hN = tspan(ti) - tNext; - end + hN = tspan(ti) - tNext; end - elseif isDense && multiTspan && tl ~= length(tspan) && ~(tspan(tl) * tdir > (tNext + hN) * tdir) - - while tNext * tdir > tspan(tl) * tdir && tl ~= length(tspan) - %Call Dense output function and record - - t(:, tl) = tspan(tl); - [y(:, tl), tempFevals] = opts.Dense.denseOut(f, tCur, tspan(tl), yCur, yNext, k, hC); - nFevals = nFevals + tempFevals; - - tl = tl + 1; - end - - if tl == length(tspan) && (tspan(ti) - tNext) * tdir < hmin * tdir - t(:, ti) = tspan(ti); - y(:, ti) = yNext; - ti = ti + 1; - end end + + end %End of Timeloop work diff --git a/src/+matlode/+startingstep/SetStep.m b/src/+matlode/+startingstep/SetStep.m index 2373fd9..352fb8b 100644 --- a/src/+matlode/+startingstep/SetStep.m +++ b/src/+matlode/+startingstep/SetStep.m @@ -16,7 +16,7 @@ fevals = 0; f0 = []; if sign(obj.StartingStep) ~= sign(tspan(end) - tspan(1)) - error('The starting step cannot go in the opposite direction') + error('The starting step cannot go in the opposite direction, %f', obj.StartingStep) end if obj.StartingStep > maxStep error('The starting step is greater then the maximum step %f', maxStep) diff --git a/src/+matlode/+stepsizecontroller/Fixed.m b/src/+matlode/+stepsizecontroller/Fixed.m index 35aa4da..df8dd8c 100644 --- a/src/+matlode/+stepsizecontroller/Fixed.m +++ b/src/+matlode/+stepsizecontroller/Fixed.m @@ -16,13 +16,10 @@ obj.Steps = steps; end - function [h0, f0] = startingStep(obj, f, tspan, y0, ~, ~, intialStep) - if intialStep ~= 0 - warning('Ignoring initialStep option, Controller is Fixed'); - end - - f0 = f(tspan(1), y0); + function [h0, f0, fevals] = startingStep(obj, ~, tspan, ~, ~, ~, ~, ~) h0 = (tspan(2) - tspan(1)) / obj.Steps; + f0 = []; + fevals = 0; end function [accept, h, tNew] = newStepSize(~, ~, ~, tspan, h, ~, ~, nSteps, ~, ~) diff --git a/src/+matlode/+stepsizecontroller/SoderlandController.m b/src/+matlode/+stepsizecontroller/SoderlandController.m index 543692f..9b7a27a 100644 --- a/src/+matlode/+stepsizecontroller/SoderlandController.m +++ b/src/+matlode/+stepsizecontroller/SoderlandController.m @@ -48,7 +48,7 @@ function [accept, hNew, tNew] = newStepSize(obj, ~, t, ~, h, err, q, ~, ~) - accept = mean(err) <= 1; + accept = err(1) <= 1; if accept tNew = t + h(1); From c26725ddeed03d73bb6696daa3c01d7c12b22b60 Mon Sep 17 00:00:00 2001 From: amitns Date: Thu, 23 Sep 2021 19:33:14 -0400 Subject: [PATCH 07/10] Added RK65SSP33. --- src/+matlode/+rk/+erk/RK65SSP33.m | 65 +++++++++++++++++++++++++++++++ 1 file changed, 65 insertions(+) create mode 100644 src/+matlode/+rk/+erk/RK65SSP33.m diff --git a/src/+matlode/+rk/+erk/RK65SSP33.m b/src/+matlode/+rk/+erk/RK65SSP33.m new file mode 100644 index 0000000..760d971 --- /dev/null +++ b/src/+matlode/+rk/+erk/RK65SSP33.m @@ -0,0 +1,65 @@ +classdef RK65SSP33 < matlode.rk.erk.ERK + properties (SetAccess = immutable) + DenseOut + end + + methods + function obj = RK65SSP33(datatype) + + % A strong stability preserving RK method + % Reference : Page 5.6, Table 4.12 a) https://core.ac.uk/download/pdf/56372752.pdf + % Macdonald, Colin Barr. "Constructing high-order Runge-Kutta methods with embedded strong-stability-preserving pairs." + % PhD diss., Theses (Dept. of Mathematics)/Simon Fraser University, 2003. + + if nargin == 0 + datatype = 'double'; + end + + caster = @(x) matlode.util.CoefficentTransformers.transform(x,datatype); + + + a = caster(join(['[[0, 0, 0, 0, 0, 0];',... + '[1, 0, 0, 0, 0, 0];',... + '[1/4, 1/4, 0, 0, 0, 0];',... + '[2046/15625, -454/15625, 1533/15625, 0, 0, 0];',... + '[-739/5625, 511/5625, -566/16875, 20/27, 0, 0];',... + '[11822/21875, -6928/21875, -4269/21875, -4/7, 54/35, 0]]'])); + + b = caster('[1/6, 1/6, 2/3, 0, 0, 0]'); + + bHat = caster('[1/24, 0, 0, 125/336, 27/56, 5/48]'); + + e = caster('[1/8, 1/6, 2/3, -125/336, -27/56, -5/48]'); + + c = caster('[0, 1, 1/2, 1/5, 2/3, 1]'); + + order = 3; + + embbededOrder = 5; + + obj = obj@matlode.rk.erk.ERK(a, b, bHat, c, e, order, embbededOrder); + + % Fix the Dense thing. + +% %Shampien formula for dense output given in the Solving ODE +% %book on pg192 (6.12) +% %the Dense output has been transformed from a function form to +% %a matrix form +% %p* == 4 s* == 5 +% denseMatrix = caster(join(['[[1, -(4034104133/1410260304), 105330401/33982176, -(13107642775/11282082432), 6542295/470086768];',... +% '[ 0, 0, 0, 0, 0];',... +% '[ 0, 132343189600/32700410799, -(833316000/131326951), 91412856700/32700410799, -(523383600/10900136933)];',... +% '[ 0, -(115792950/29380423), 185270875/16991088, -(12653452475/1880347072), 98134425/235043384];',... +% '[ 0, 70805911779/24914598704, -(4531260609/600351776), 988140236175/199316789632, -(14307999165/24914598704)];',... +% '[ 0, -(331320693/205662961), 31361737/7433601, -(2426908385/822651844), 97305120/205662961];',... +% '[ 0, 44764047/29380423, -(1532549/353981), 90730570/29380423, -(8293050/29380423)]]'])); +% obj.DenseOut = matlode.denseoutput.RKDense(denseMatrix); + + end + + function sol = integrate(obj, f, tspan, y0, varargin) + sol = integrate@matlode.rk.erk.ERK(obj, f, tspan, y0, 'StepSizeController', matlode.stepsizecontroller.StandardController, 'Dense', obj.DenseOut, varargin{:}); + end + end +end + From 80df554a36e762e9562c8ebe3e4b15d33eca7d82 Mon Sep 17 00:00:00 2001 From: amitns Date: Thu, 23 Sep 2021 19:47:27 -0400 Subject: [PATCH 08/10] Updated RK65SSP33 --- src/+matlode/+rk/+erk/RK65SSP33.m | 21 +-------------------- 1 file changed, 1 insertion(+), 20 deletions(-) diff --git a/src/+matlode/+rk/+erk/RK65SSP33.m b/src/+matlode/+rk/+erk/RK65SSP33.m index 760d971..372b876 100644 --- a/src/+matlode/+rk/+erk/RK65SSP33.m +++ b/src/+matlode/+rk/+erk/RK65SSP33.m @@ -1,7 +1,4 @@ classdef RK65SSP33 < matlode.rk.erk.ERK - properties (SetAccess = immutable) - DenseOut - end methods function obj = RK65SSP33(datatype) @@ -38,27 +35,11 @@ embbededOrder = 5; obj = obj@matlode.rk.erk.ERK(a, b, bHat, c, e, order, embbededOrder); - - % Fix the Dense thing. - -% %Shampien formula for dense output given in the Solving ODE -% %book on pg192 (6.12) -% %the Dense output has been transformed from a function form to -% %a matrix form -% %p* == 4 s* == 5 -% denseMatrix = caster(join(['[[1, -(4034104133/1410260304), 105330401/33982176, -(13107642775/11282082432), 6542295/470086768];',... -% '[ 0, 0, 0, 0, 0];',... -% '[ 0, 132343189600/32700410799, -(833316000/131326951), 91412856700/32700410799, -(523383600/10900136933)];',... -% '[ 0, -(115792950/29380423), 185270875/16991088, -(12653452475/1880347072), 98134425/235043384];',... -% '[ 0, 70805911779/24914598704, -(4531260609/600351776), 988140236175/199316789632, -(14307999165/24914598704)];',... -% '[ 0, -(331320693/205662961), 31361737/7433601, -(2426908385/822651844), 97305120/205662961];',... -% '[ 0, 44764047/29380423, -(1532549/353981), 90730570/29380423, -(8293050/29380423)]]'])); -% obj.DenseOut = matlode.denseoutput.RKDense(denseMatrix); end function sol = integrate(obj, f, tspan, y0, varargin) - sol = integrate@matlode.rk.erk.ERK(obj, f, tspan, y0, 'StepSizeController', matlode.stepsizecontroller.StandardController, 'Dense', obj.DenseOut, varargin{:}); + sol = integrate@matlode.rk.erk.ERK(obj, f, tspan, y0, 'StepSizeController', matlode.stepsizecontroller.StandardController, 'Dense', matlode.denseoutput.Linear(obj.B), varargin{:}); end end end From f4866c85abc8a96206da58a49054c5c4a0a88cc6 Mon Sep 17 00:00:00 2001 From: Andrey A Popov Date: Fri, 17 Dec 2021 19:16:41 -0500 Subject: [PATCH 09/10] minor fix to spacing --- src/+matlode/+stepsizecontroller/StartingController.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/+matlode/+stepsizecontroller/StartingController.m b/src/+matlode/+stepsizecontroller/StartingController.m index 87f6093..4300320 100644 --- a/src/+matlode/+stepsizecontroller/StartingController.m +++ b/src/+matlode/+stepsizecontroller/StartingController.m @@ -30,7 +30,7 @@ opts = p.Results; varargout = p.Unmatched; - obj = obj@matlode.stepsizecontroller.StepSizeController(hist, varargout); + obj = obj@matlode.stepsizecontroller.StepSizeController(hist, varargout); obj.Fac = opts.Fac; obj.FacMin = opts.FacMin; From 78ca8a6d597a5dc047ce13e3fb7f2b94004da2f0 Mon Sep 17 00:00:00 2001 From: Andrey A Popov Date: Thu, 30 Dec 2021 19:04:07 -0500 Subject: [PATCH 10/10] updated the warning to contain an errID --- src/+matlode/+rk/+erk/ERK.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/+matlode/+rk/+erk/ERK.m b/src/+matlode/+rk/+erk/ERK.m index 48e24ea..8b8ffe1 100644 --- a/src/+matlode/+rk/+erk/ERK.m +++ b/src/+matlode/+rk/+erk/ERK.m @@ -169,7 +169,7 @@ %Check if step is really small if abs(hN) < abs(hmin) if nSmallSteps == 0 - warning('The step the integrator is taking extremely small, results may not be optimal') + warning('MATLODE:smallStepSize', 'The step the integrator is taking extremely small, results may not be optimal') end %accept step since the step cannot get any smaller nSmallSteps = nSmallSteps + 1;