From e1778374fed1ad8923067c350b7b091004f037fb Mon Sep 17 00:00:00 2001 From: Steven Roberts Date: Sat, 20 Apr 2024 22:16:01 -0700 Subject: [PATCH 1/5] Add docs to Problem --- .../+utils/+movie/+recorder/FileRecorder.m | 4 +- .../+utils/+movie/+recorder/MemoryRecorder.m | 8 +- .../+utils/+movie/+recorder/NullRecorder.m | 5 +- toolbox/+otp/+utils/+movie/Movie.m | 18 +- toolbox/+otp/Problem.m | 398 ++++++++++++++++-- toolbox/+otp/RHS.m | 76 ++-- 6 files changed, 428 insertions(+), 81 deletions(-) diff --git a/toolbox/+otp/+utils/+movie/+recorder/FileRecorder.m b/toolbox/+otp/+utils/+movie/+recorder/FileRecorder.m index dfe5ba0f..7e28cf27 100644 --- a/toolbox/+otp/+utils/+movie/+recorder/FileRecorder.m +++ b/toolbox/+otp/+utils/+movie/+recorder/FileRecorder.m @@ -43,8 +43,8 @@ function stop(obj) obj.VideoWriter.close(); end - function h = play(obj) - h = implay(fullfile(obj.VideoWriter.Path, obj.VideoWriter.Filename), obj.VideoWriter.FrameRate); + function play(obj) + implay(fullfile(obj.VideoWriter.Path, obj.VideoWriter.Filename), obj.VideoWriter.FrameRate); end end end diff --git a/toolbox/+otp/+utils/+movie/+recorder/MemoryRecorder.m b/toolbox/+otp/+utils/+movie/+recorder/MemoryRecorder.m index fbd6b56d..ba106d39 100644 --- a/toolbox/+otp/+utils/+movie/+recorder/MemoryRecorder.m +++ b/toolbox/+otp/+utils/+movie/+recorder/MemoryRecorder.m @@ -29,8 +29,12 @@ function stop(~) % Nothing to do end - function h = play(obj) - h = implay(obj.Mov, obj.FrameRate); + function play(obj) + if exist('implay', 'builtin') + implay(obj.Mov, obj.FrameRate); + else + movie(obj.Mov, 1, obj.FrameRate); + end end end end diff --git a/toolbox/+otp/+utils/+movie/+recorder/NullRecorder.m b/toolbox/+otp/+utils/+movie/+recorder/NullRecorder.m index 3a787738..d974c3aa 100644 --- a/toolbox/+otp/+utils/+movie/+recorder/NullRecorder.m +++ b/toolbox/+otp/+utils/+movie/+recorder/NullRecorder.m @@ -24,9 +24,8 @@ function stop(~) % Nothing to do end - function h = play(~) - h = 'Movie not saved for playback'; - error('OTP:movieNotSaved', h); + function play(~) + error('OTP:movieNotSaved', 'Movie not saved for playback'); end end end diff --git a/toolbox/+otp/+utils/+movie/Movie.m b/toolbox/+otp/+utils/+movie/Movie.m index 43c84c7d..4b7679d8 100644 --- a/toolbox/+otp/+utils/+movie/Movie.m +++ b/toolbox/+otp/+utils/+movie/Movie.m @@ -17,7 +17,7 @@ p = inputParser; p.addParameter('Save', false); p.addParameter('FrameRate', otp.utils.movie.Movie.DefaultFramerate); - p.addParameter('TargetDuration', [], @(d) d > 0); + p.addParameter('Duration', [], @(d) d > 0); p.addParameter('Size', [], @(s) length(s) == 2 && all(s > 0)); p.addParameter('Smooth', true, @islogical); p.parse(varargin{:}); @@ -46,9 +46,7 @@ function record(obj, t, y) totalSteps = length(t); [state.numVars, state.totalSteps] = size(y); if length(t) ~= state.totalSteps - error('OTP:invalidSolution', ... - 'Expected y to have %d columns but has %d', ... - length(t), state.totalSteps); + error('OTP:invalidSolution', 'Expected y to have %d columns but has %d', length(t), state.totalSteps); end state.t = t; @@ -56,17 +54,19 @@ function record(obj, t, y) state.step = 0; state.frame = 0; - if isempty(obj.Config.TargetDuration) + if isempty(obj.Config.Duration) state.totalFrames = totalSteps; else - state.totalFrames = round(obj.Config.TargetDuration * obj.FrameRate); + state.totalFrames = round(obj.Config.Duration * obj.FrameRate); end [t0, tEnd] = bounds(t); fig = figure; if ~isempty(obj.Config.Size) - fig.Position = [0; 0; obj.Config.Size(:)]; + pos = get(fig, 'position'); + pos(3:4) = obj.Config.Size; + set(fig, 'position', pos); end obj.init(fig, state); @@ -97,8 +97,8 @@ function record(obj, t, y) obj.Recorder.stop(); end - function h = play(obj) - h = obj.Recorder.play(); + function play(obj) + obj.Recorder.play(); end end diff --git a/toolbox/+otp/Problem.m b/toolbox/+otp/Problem.m index 19e8d6bb..bd4cd3de 100644 --- a/toolbox/+otp/Problem.m +++ b/toolbox/+otp/Problem.m @@ -1,19 +1,52 @@ classdef (Abstract) Problem < handle + % The superclass for all problems in OTP. + % + % A problem is defined by a differential equation + % + % $$M(t, y(t); p) y'(t) = f(t, y(t); p)$$ + % + % a time span $t ∈ [t_0, t_f]$, initial conditions $y(t_0) = y_0$, and any parameters $p$. The unknown solution + % function is $y(t) ∈ ℂ^N$. + % + % This class provides the basic infrastructure for defining, solving, plotting, and animating a problem. + % + % See Also + % -------- + % otp.RHS + % otp.Parameters + % + % ode : https://www.mathworks.com/help/matlab/ref/ode.html properties (SetAccess = private) - % A human-readable representation of the problem + % A human-readable representation of the problem. Name %MATLAB ONLY: (1,:) char end properties (SetAccess = protected) + % The right-hand side and related properties of the differential equation. + % + % This property is an :class:`otp.RHS` which is regenerated when the problem's time span, initial conditions, or + % parameters change. + % + % Example + % ------- + % >>> problem = otp.linear.presets.Canonical; + % >>> problem.RHS.F + % ans = + % + % @(~, y) lambda * y + % + % See Also + % -------- + % otp.RHS RHS %MATLAB ONLY: (1,1) otp.RHS = otp.RHS(@() []) end properties (Access = private) + % The expected number of variables or an empty array for an arbitrary number. ExpectedNumVars %MATLAB ONLY: {mustBeScalarOrEmpty, mustBeInteger, mustBePositive} - % Determines if the Problem properties are still being set. Once set, it - % is ok to call onSettingsChanged. + % Determines if the problem properties are still being set. Once set, it is ok to call onSettingsChanged. Initialized = false % The following internal properties represent the state of the problem @@ -27,22 +60,75 @@ end properties (Dependent) - % The time interval of the problem + % A mutable time interval $[t_0, t_f]$ over which the problem is defined. + % + % Example + % ------- + % >>> problem = otp.protherorobinson.presets.Canonical; + % >>> problem.TimeSpan + % >>> problem.TimeSpan(2) = 100; + % ans = + % + % 0 15 TimeSpan %MATLAB ONLY: (1,2) {otp.utils.validation.mustBeNumerical} - % The initial condition + % A mutable column vector for the initial conditions $y(t_0) = y_0$. + % + % Example + % ------- + % >>> problem = otp.lotkavolterra.presets.Canonical; + % >>> problem.Y0 + % >>> problem.Y0 = [3; 4]; + % ans = + % + % 1 + % 2 Y0 %MATLAB ONLY: (:,1) {otp.utils.validation.mustBeNumerical} - % Additional variables to pass to the F function + % A mutable :class:`otp.Parameters` representing the problem parameters $p$. + % + % Example + % ------- + % >>> problem = otp.nbody.presets.Canonical; + % >>> problem.Parameters.GravitationalConstant + % >>> problem.Parameters.Masses(1) = 2; + % ans = 1 + % + % See Also + % -------- + % otp.Parameters Parameters %MATLAB ONLY: (1,1) otp.Parameters - % The dimension of the ODE + % The number of variables $N$ of the problem. + % + % Example + % ------- + % >>> problem = otp.cusp.presets.Canonical; + % >>> problem.NumVars + % >>> problem.Y0 = ones(30, 1); + % >>> problem.NumVars + % ans = 96 + % ans = 30 NumVars end methods function obj = Problem(name, expectedNumVars, timeSpan, y0, parameters) - % Constructs a problem + % Create a problem object. + % + % Parameters + % ---------- + % name : char(1, :) + % A human-readable representation of the problem. + % expectedNumVars : numeric or [] + % The expected number of variables or an empty array for an arbitrary number. + % timeSpan : numeric(1, 2) + % The start and final time. + % y0 : numeric(:, 1) + % The initial conditions. + % parameters : otp.Parameters + % The parameters. + obj.Name = name; obj.ExpectedNumVars = expectedNumVars; obj.TimeSpan = timeSpan; @@ -65,9 +151,7 @@ function set.Y0(obj, value) numVars = length(value); if ~isempty(obj.ExpectedNumVars) && numVars ~= obj.ExpectedNumVars - error('OTP:invalidY0', ... - 'Expected Y0 to have %d components but has %d', ... - obj.ExpectedNumVars, numVars); + error('OTP:invalidY0', 'Expected Y0 to have %d components but has %d', obj.ExpectedNumVars, numVars); end obj.InternalY0 = value; if obj.Initialized @@ -95,7 +179,70 @@ methods (Sealed) function varargout = plot(obj, sol, varargin) - % Plots all trajectories y versus time + % Plot all solution trajectories with time on the x-axis and $y(t)$ on the y-axis. + % + % Each call to this function creates a new figure and displays all solution variables. It accepts as input + % the two standard output formats of MATLAB's differential equation solvers: a solution structure or a + % ``[t, y]`` pair. + % + % Warning + % ------- + % A plot might appear jagged if there are too few times included in the solution. In MATLAB, + % `deval `_ can be used to evaluate the solution at a + % user-specified set of points. + % + % Parameters + % ---------- + % sol : struct or numeric(1, :) or numeric(:, 1) + % A solution structure with properties ``x`` and ``y``. The length of one dimension of ``y`` must match + % the number of times in ``x``. If ``y`` is square, each column should be a state and each row the + % trajectory of a single variable. Alternatively, this argument can be a vector of times. + % y : numeric(:, :), optional + % If the first argument is a vector of times, this should be the corresponding matrix of solutions. The + % length of one dimension must match the number of times. If square, each row should be a state and each + % column the trajectory of a single variable. + % varargin + % A variable number of name-value pairs. Accepted names include + % + % - ``Title`` – A string for the figure title. + % - ``XLabel`` – A string for the x-axis label. + % - ``YLabel`` – A string for the y-axis label. + % - ``ZLabel`` – A string for the z-axis label. + % - ``View`` – An array specifying the camera line of sight following + % `MATLAB's specification `__. + % - ``ColorOrder`` – The color order palette following + % `MATLAB's specification `__. + % - ``Axis`` – The axis limits and aspect ratios following + % `MATLAB's specification `__. + % - ``LineStyleOrder`` – The line style order following + % `MATLAB's specification `__. + % - ``XScale`` – Either ``'linear'`` or ``'log'``. + % - ``YScale`` – Either ``'linear'`` or ``'log'``. + % - ``ZScale`` – Either ``'linear'`` or ``'log'``. + % - ``FontName`` – A string for the font name. + % - ``FontSize`` – A size in the unit of points for the font. + % - ``Legend`` – A cell array of strings or a function handle with the signature ``label = fun(index)`` + % specifying the label of lines to show in the legend. + % - ``MaxLegendLabels`` – The maximum number of entries to show in the legend. If the number of lines + % exceeds this, only an evenly-spaced subset appears in the legend. + % + % Returns + % ------- + % fig : figure, optional + % The figure created for the plot. + % + % Examples + % -------- + % >>> problem = otp.trigonometricdae.presets.Canonical; + % >>> sol = problem.solve(); + % >>> problem.plot(sol, 'Axis', 'square', 'XLabel', 'Time'); + % >>> x = linspace(p.TimeSpan(1), p.TimeSpan(2), 500); + % >>> fig = problem.plot(x, deval(x, sol), 'Title', 'Another Plot'); + % + % >>> problem = otp.vanderpol.presets.Canonical; + % >>> [t, y] = ode15s(problem.RHS.F, problem.TimeSpan, problem.Y0); + % >>> problem.plot(t, y); + [t, y, params] = obj.parseSolution(sol, varargin{:}); fig = obj.internalPlot(t, y, params{:}); @@ -105,7 +252,47 @@ end function varargout = plotPhaseSpace(obj, sol, varargin) - % Plots a selection of trajectories with respect to each other + % Plot a selection of solution trajectories with respect to each other. + % + % A phase plot is 2D or 3D with one or more lines. For each line, the x-, y-, and z-axis can represent a + % different user-specified component of the solution. Each call to this function creates a new figure. It + % accepts as input the two standard output formats of MATLAB's differential equation solvers: a solution + % structure or a ``[t, y]`` pair. + % + % Warning + % ------- + % A plot might appear jagged if there are too few times included in the solution. In MATLAB, `deval`_ can be + % used to evaluate the solution at a user-specified set of points. + % + % Parameters + % ---------- + % sol : struct or numeric(1, :) or numeric(:, 1) + % A solution structure with properties ``x`` and ``y``. The length of one dimension of ``y`` must match + % the number of times in ``x``. If ``y`` is square, each column should be a state and each row the + % trajectory of a single variable. Alternatively, this argument can be a vector of times. + % y : numeric(:, :), optional + % If the first argument is a vector of times, this should be the corresponding matrix of solutions. The + % length of one dimension must match the number of times. If square, each row should be a state and each + % column the trajectory of a single variable. + % varargin + % A variable number of name-value pairs. Accepted names include those supported by + % :func:`otp.Problem.plot` as well as the following: + % + % - ``Vars`` – a matrix of variable indices where each row corresponds to a line in the plot and columns + % one, two, and optionally three correspond to the x-, y-, and z-axis, respectively. + % + % Returns + % ------- + % fig : figure, optional + % The figure created for the plot. + % + % Examples + % -------- + % >>> problem = otp.lorenz63.presets.Canonical; + % >>> sol = problem.solve(); + % >>> problem.plotPhaseSpace(sol, 'Title', 'Lorenz Butterfly'); + % >>> problem.plotPhaseSpace(sol, 'Vars', [1, 2; 2, 3; 3, 1]); + [t, y, params] = obj.parseSolution(sol, varargin{:}); fig = obj.internalPlotPhaseSpace(t, y, params{:}); @@ -115,7 +302,63 @@ end function varargout = movie(obj, sol, varargin) - % Plots an animation of the trajectories + % Create an animation of solution trajectories. + % + % Each call to this function creates a new figure. It accepts as input the two standard output formats of + % MATLAB's differential equation solvers: a solution structure or a ``[t, y]`` pair. + % + % Warning + % ------- + % Octave has limited movie functionality and support. For most problems, the movie will not work. Saving a + % movie to a file is currently unsupported. + % + % Warning + % ------- + % A movie might appear jerky if there are too few times included in the solution. In MATLAB, `deval`_ can be + % used to evaluate the solution at a user-specified set of points. + % + % Parameters + % ---------- + % sol : struct or numeric(1, :) or numeric(:, 1) + % A solution structure with properties ``x`` and ``y``. The length of one dimension of ``y`` must match + % the number of times in ``x``. If ``y`` is square, each column should be a state and each row the + % trajectory of a single variable. Alternatively, this argument can be a vector of times. + % y : numeric(:, :), optional + % If the first argument is a vector of times, this should be the corresponding matrix of solutions. The + % length of one dimension must match the number of times. If square, each row should be a state and each + % column the trajectory of a single variable. + % varargin + % A variable number of name-value pairs. Accepted names include those supported by + % :func:`otp.Problem.plot` as well as the following: + % + % - ``Save`` – If ``false`` (default) the movie will play but the frames will not be save for playback. + % If ``true`` the frames in memory. If a string or a + % `VideoWriter `__, the movie will be + % written to a file. + % - ``FrameRate`` – The number of frames per second when playing the movie. + % - ``Duration`` – The length in seconds of the movie. Frames may need to be skipped or duplicated to + % accommodate. + % - ``Size`` – An array ``[width, height]`` for the figure size in pixels. + % - ``Smooth`` – If ``true`` (default), solution states are uniformly sampled for frames based on the + % corresponding solution times. If ``false`` solution states are sampled uniformly by time step index. + % This may lead to the inconsistent playback speeds. + % + % Returns + % ------- + % mov : otp.utils.movie.Movie, optional + % A movie object which can be replayed using the ``play()`` method if saved. + % + % Examples + % -------- + % >>> problem = otp.pendulum.presets.Canonical('NumBobs', 8); + % >>> sol = problem.solve(); + % >>> mov = problem.movie(sol, 'Save', true, 'Duration', 5); + % >>> mov.play(); + % + % >>> problem = otp.bouncingball.presets.RandomTerrain; + % >>> sol = problem.solve(); + % >>> mov = problem.movie(sol, 'Save', 'movie.avi', 'FrameRate', 30); + [t, y, params] = obj.parseSolution(sol, varargin{:}); mov = obj.internalMovie(t, y, params{:}); @@ -125,20 +368,100 @@ end function label = index2label(obj, index) - % Gets a human-readable label for a particular component of the ODE + % Get a human-readable label for a particular component of the differential equation + % + % Parameters + % ---------- + % index + % An index $i$ between $1$ and $N$ for component $y_i(t)$. + % + % Returns + % ------- + % label : char(1, :) + % The name of the component. + % + % Example + % ------- + % >>> problem = otp.brusselator.presets.Canonical; + % >>> problem.index2label(2) + % ans = Reactant Y + if ~isscalar(index) || mod(index, 1) || index < 1 || index > obj.NumVars - error('OTP:indexOutOfBounds', ... - 'The index %d is not an integer between 1 and %d', ... - index, obj.NumVars); + error('OTP:indexOutOfBounds', 'The index is not an integer between 1 and %d', obj.NumVars); end label = obj.internalIndex2label(index); end function sol = solve(obj, varargin) + % Numerically compute a solution to the differential equation. + % + % This function will numerically integrate the differential equation until a terminal event is triggered, + % the end of the time span is reached, or the solver throws an error. It is handy for generating a reference + % solution or studying convergence of an integrator with respect to the tolerance. + % + % Note + % ---- + % Different problem subclasses use different default solvers. The solver used is stored in the ``solver`` + % property of the returned solution structure. + % + % Warning + % ------- + % ``ode15s`` in Octave transposes the ``y`` property of the returned solution structure. All Octave solvers + % transpose event data compared to MATLAB. + % + % Parameters + % ---------- + % varargin + % A variable number of name-value pairs. Accepted names include those supported by + % :func:`otp.RHS.odeset` as well as the following: + % + % - ``Solver`` – A function handle to a differential equation solver which can be evaluated as + % ``sol = Solver(odefun, tspan, y0, options)``. + % + % Returns + % ------- + % sol : struct + % A solution structure in the standard MATLAB solver output format. The time and solution steps are + % stored in the ``x`` and ``y`` properties, respectively. + % + % Example + % ------- + % >>> problem = otp.pendulum.presets.Canonical; + % >>> sol1 = problem.solve(); + % >>> sol2 = problem.solve('RelTol', 1e-8, 'Solver', @ode23); + + % No validation is needed, so directly call the protected internal solve function sol = obj.internalSolve(varargin{:}); end function y = solveExactly(obj, t) + % Compute the exact solution to the differential equation if available or throw an error otherwise. + % + % Parameters + % ---------- + % t : numeric or numeric(1, :) or numeric(:, 1), optional + % The times at which to generate the exact solution. If not provided, the end of the time span $t_f$ is + % used. + % + % Returns + % ------- + % y : numeric(:, :) + % A matrix in which element $(i, j)$ is $y_i(t_j)$. + % + % Examples + % -------- + % >>> problem = otp.ascherlineardae.presets.Canonical; + % >>> problem.solveExactly() + % ans = + % + % 1.5772 + % 1.2094 + % + % >>> problem = otp.protherorobinson.presets.Canonical; + % >>> sol = problem.solve(); + % >>> err = sol.y - problem.solveExactly(sol.x); + % >>> problem.plot(sol.x, err, 'YLabel', 'Error'); + if nargin < 2 t = obj.TimeSpan(2); else @@ -149,14 +472,16 @@ end methods (Access = protected) - % This method is called when either TimeSpan, Y0, or parameters are changed. It should - % update F and other properties such as a Jacobian to reflect the changes. function onSettingsChanged(obj) + % This method is called when either TimeSpan, Y0, or parameters are changed. It should update all RHS` to + % reflect the change. + otp.utils.compatibility.abstract(obj); end function fig = internalPlot(obj, t, y, varargin) % Plots all trajectories in y versus time assuming t and y are valid + fig = figure; ax = axes(fig); otp.utils.FancyPlot.plot(ax, t, y.', ... @@ -168,7 +493,7 @@ function onSettingsChanged(obj) end function fig = internalPlotPhaseSpace(obj, ~, y, varargin) - % Plots a selection of trajectories with respect to each other assuming t and y are valid + % Plots a selection of trajectories with respect to each other assuming t and y are valid. p = inputParser; p.KeepUnmatched = true; @@ -178,8 +503,7 @@ function onSettingsChanged(obj) [numLines, dim] = size(vars); if dim < otp.utils.PhysicalConstants.TwoD || dim > otp.utils.PhysicalConstants.ThreeD - error('OTP:invalidPhaseDimension', ... - 'Cannot plot a %dD phase space', dim); + error('OTP:invalidPhaseDimension', 'Cannot plot a %dD phase space', dim); end fig = figure; @@ -215,16 +539,21 @@ function onSettingsChanged(obj) end function mov = internalMovie(obj, t, y, varargin) - mov = otp.utils.movie.TrajectoryMovie('title', obj.Name, varargin{:}); + % Creates an animation assuming t and y are valid. + + mov = otp.utils.movie.TrajectoryMovie('Title', obj.Name, varargin{:}); mov.record(t, y); end function label = internalIndex2label(~, index) % Gets a human-readable label for a particular component of the ODE assuming index is a valid integer + label = sprintf('y_{%d}', index); end function sol = internalSolve(obj, varargin) + % Numerically compute a solution to the differential equation. + p = inputParser; p.KeepUnmatched = true; % Filter name-value pairs not passed to odeset @@ -258,6 +587,9 @@ function onSettingsChanged(obj) end function y = internalSolveExactly(~, ~) + % Throws an error that an exact solution is unavailable. Subclass can override this to provide an exact + % solution when available. + y = 'This problem does not provide an exact solution'; error('OTP:noExactSolution', y); end @@ -265,9 +597,11 @@ function onSettingsChanged(obj) methods (Access = private) function [t, y, params] = parseSolution(obj, sol, varargin) + % Puts the two types of ODE solver output (solution structure or [t, y]) into a consistent format + if isstruct(sol) t = sol.x; - y = sol.y; + y = sol.y.'; params = varargin; else t = sol; @@ -287,21 +621,19 @@ function onSettingsChanged(obj) if m == steps && n == obj.NumVars && m ~= n y = y.'; elseif m ~= obj.NumVars - error('OTP:invalidSolution', ... - 'There are %d timesteps, but %d solution states', steps, m); + error('OTP:invalidSolution', 'There are %d time steps, but %d solution states', steps, m); elseif n ~= steps - error('OTP:invalidSolution', ... - 'Expected solution to have %d variables but has %d', ... - obj.NumVars, n); + error('OTP:invalidSolution', 'Expected solution to have %d variables but has %d', obj.NumVars, n); end end function t = parseTime(~, t) + % Puts a vector of times into a consistent row vector form + if ~isvector(t) || isempty(t) || ~otp.utils.validation.isNumerical(t) - error('OTP:invalidSolution', ... - 'The times must be a nonempty vector of numbers'); + error('OTP:invalidSolution', 'The times must be a nonempty vector of numbers'); end - % Convert to row vector for consistency + t = t(:).'; end end diff --git a/toolbox/+otp/RHS.m b/toolbox/+otp/RHS.m index 42c610b6..6f1a6337 100644 --- a/toolbox/+otp/RHS.m +++ b/toolbox/+otp/RHS.m @@ -4,16 +4,22 @@ % This immutable class contains properties required by time integrators and other numerical methods to describe the % following differential equation: % - % $$M(t, y) y' = f(t, y)$$ + % $$M(t, y(t); p) y'(t) = f(t, y(t); p)$$ % - % The mass matrix $M(t, y)$ is permitted to be singular in which case the problem is a differential-algebraic - % equation. ``RHS`` includes most of the properties available in - % `odeset `_ so that it can easily be used with built-in ODE - % solvers. + % Despite the name, an ``RHS`` includes information about the mass matrix $M(t, y(t); p)$ on the left-hand side. The + % mass matrix is permitted to be singular in which case the problem is a differential-algebraic equation. ``RHS`` + % also includes most of the properties available in MATLAB`s `odeset`_ so that it can easily be used with built-in + % solvers like ``ode15s``. + % + % Note + % ---- + % The parameters $p$ do not appear explicitly in the arguments of functions in this class. Rather, they can be + % provided through the closure of an anonymous function. % % Examples % -------- - % >>> rhs = otp.RHS(@(t, y) [y(2) + t; y(1) * y(2)], ... + % >>> param = 1; + % >>> rhs = otp.RHS(@(t, y) [y(2) + param * t; y(1) * y(2)], ... % ... 'Jacobian', @(~, y) [0, 1; y(2), y(1)], ... % ... 'Mass', [1, 2; 3, 4]); % >>> sol = ode23s(rhs.F, [0, 1], [1, -1], rhs.odeset()); @@ -35,7 +41,7 @@ % odeset : https://www.mathworks.com/help/matlab/ref/odeset.html properties (SetAccess = private) - % The function handle for $f$ in the differential equation $M(t, y) y' = f(t, y)$. + % The function handle for $f$ in the differential equation. % % Parameters % ---------- @@ -55,7 +61,7 @@ % % This follows the same specifications as in `odeset `_. % If set, it is a matrix if it is independent of t and y or a function handle. In either case, it provides a - % square matrix in which element $(i,j)$ is $\frac{\partial f_i(t, y)}{\partial y_j}$. If ``Jacobian`` is a + % square matrix in which element $(i,j)$ is $\frac{\partial f_i(t, y; p)}{\partial y_j}$. If ``Jacobian`` is a % function handle, it has the following signature: % % Parameters @@ -78,12 +84,12 @@ % The sparsity pattern of the Jacobian matrix. % - % This follows the same specifications as in `odeset `_. + % This follows the same specifications as in `odeset`_. JPattern - % The mass matrix $M(t, y)$. + % The mass matrix $M(t, y(t); p)$. % - % This follows the same specifications as in `odeset `_. + % This follows the same specifications as in `odeset`_. % % See Also % -------- @@ -93,27 +99,27 @@ % Whether the mass matrix is singular, i.e., the problem is a differential-algebraic equation. % - % This follows the same specifications as in `odeset `_. + % This follows the same specifications as in `odeset`_. MassSingular % State dependence of the mass matrix. % - % This follows the same specifications as in `odeset `_. + % This follows the same specifications as in `odeset`_. MStateDependence - % The sparsity pattern of $\frac{\partial}{\partial y} M(t, y) v$. + % The sparsity pattern of $\frac{\partial}{\partial y} M(t, y; p) v$. % - % This follows the same specifications as in `odeset `_. + % This follows the same specifications as in `odeset`_. MvPattern % A vector of solution components which should be nonnegative. % - % This follows the same specifications as in `odeset `_. + % This follows the same specifications as in `odeset`_. NonNegative % Whether $f$ can be evaluated at multiple states. % - % This follows the same specifications as in `odeset `_. + % This follows the same specifications as in `odeset`_. % % Warning % ------- @@ -122,10 +128,10 @@ % A function which detects events and determines whether to terminate the integration. % - % This follows the same specifications as in `odeset `_. + % This follows the same specifications as in `odeset`_. Events - % Response to a terminal event occuring. + % Response to a terminal event occurring. % % If set, it is a function handle with the following signature: % @@ -134,7 +140,7 @@ % sol : struct % The solution generated by an ODE solver once an event occurs. % problem : Problem - % The problem for which the event occured. This should not be modified by the function. + % The problem for which the event occurred. This should not be modified by the function. % % Returns % ------- @@ -147,7 +153,7 @@ % Warning % ------- % With Octave solvers, event data in ``sol`` is transposed compared to MATLAB. Therefore, we recommend accessing - % the state at which the event occured with ``sol.y(:, end)`` as opposed to ``sol.ye(:, end)``. + % the state at which the event occurred with ``sol.y(:, end)`` as opposed to ``sol.ye(:, end)``. OnEvent % The action of the Jacobian multiplying a vector. @@ -167,7 +173,7 @@ % Returns % ------- % jvp : numeric(:, 1) - % A vector in which element $i$ is $\sum_j \frac{\partial f_i(t, y)}{\partial y_j} v_j$. + % A vector in which element $i$ is $\sum_j \frac{\partial f_i(t, y; p)}{\partial y_j} v_j$. JacobianVectorProduct % The action of the adjoint (conjugate transpose) of the Jacobian multiplying a vector. @@ -187,7 +193,7 @@ % Returns % ------- % javp : numeric(:, 1) - % A vector in which element $i$ is $\sum_j \frac{\partial \overline{f_j(t, y)}}{\partial y_i} v_j$. + % A vector in which element $i$ is $\sum_j \frac{\partial \overline{f_j(t, y; p)}}{\partial y_i} v_j$. JacobianAdjointVectorProduct % The partial derivative of $f$ with respect to parameters. @@ -222,7 +228,7 @@ % Returns % ------- % pdt : numeric(:, 1) - % A vector in which element $i$ is $\frac{\partial f_i(t, y)}{\partial t}$. + % A vector in which element $i$ is $\frac{\partial f_i(t, y; p)}{\partial t}$. % % Note % ---- @@ -246,7 +252,7 @@ % ------- % hvp : numeric(:, 1) % A vector in which element $i$ is - % $\sum_{j,k} \frac{\partial^2 f_i(t, y)}{\partial y_j \partial y_k} u_j v_k$. + % $\sum_{j,k} \frac{\partial^2 f_i(t, y; p)}{\partial y_j \partial y_k} u_j v_k$. HessianVectorProduct % The action of a vector on the partial derivative of :attr:`JacobianAdjointVectorProduct` with respect to $y$. @@ -266,18 +272,18 @@ % ------- % havp : numeric(:, 1) % A vector in which element $i$ is - % $\sum_{j,k} \frac{\partial^2 \overline{f_j(t, y)}}{\partial y_i \partial y_k} u_j v_k$. + % $\sum_{j,k} \frac{\partial^2 \overline{f_j(t, y; p)}}{\partial y_i \partial y_k} u_j v_k$. HessianAdjointVectorProduct end properties (Dependent) - % A dependent property which retruns :attr:`Jacobian` if it is a matrix and ``[]`` if it is a function handle. + % A dependent property which returns :attr:`Jacobian` if it is a matrix and ``[]`` if it is a function handle. JacobianMatrix % A dependent property which wraps :attr:`Jacobian` in a function handle if necessary. JacobianFunction - % A dependent property which retruns :attr:`Mass` if it is a matrix and ``[]`` if it is a function handle. + % A dependent property which returns :attr:`Mass` if it is a matrix and ``[]`` if it is a function handle. MassMatrix % A dependent property which wraps :attr:`Mass` in a function handle if necessary. @@ -291,7 +297,7 @@ % Parameters % ---------- % F : function_handle - % The function handle for $f$ in the differential equation $M(t, y) y' = f(t, y)$. + % The function handle for $f$ in the differential equation. % varargin % A variable number of name-value pairs. A name can be any property of this class, and the subsequent % value initializes that property. @@ -381,17 +387,23 @@ end function opts = odeset(obj, varargin) - % Convert an ``RHS`` to a struct which can be used by built-in ODE solvers. + % Convert an ``RHS`` to a options structure which can be used by built-in ODE solvers. % % Parameters % ---------- % varargin - % Additional name-value pairs which have precedence over the values in this class. + % Additional name-value pairs which have precedence over the values in this class when constructing the + % options structure. Accept names include those supported by the built-in ``odeset``. % % Returns % ------- % opts : struct - % An options strucutre containing the subset of properties compatible with built-in ODE solvers. + % An options structure containing the subset of properties compatible with built-in ODE solvers. + % + % See Also + % -------- + % + % odeset : https://www.mathworks.com/help/matlab/ref/odeset.html opts = odeset( ... 'Events', obj.Events, ... From 062649ab08120852308ba9fe16eae6dae8d183b3 Mon Sep 17 00:00:00 2001 From: Steven Roberts Date: Tue, 1 Oct 2024 19:31:17 -0700 Subject: [PATCH 2/5] Apply suggestions from code review --- toolbox/+otp/+utils/+movie/+recorder/MemoryRecorder.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/toolbox/+otp/+utils/+movie/+recorder/MemoryRecorder.m b/toolbox/+otp/+utils/+movie/+recorder/MemoryRecorder.m index ba106d39..b09b8d60 100644 --- a/toolbox/+otp/+utils/+movie/+recorder/MemoryRecorder.m +++ b/toolbox/+otp/+utils/+movie/+recorder/MemoryRecorder.m @@ -30,7 +30,7 @@ function stop(~) end function play(obj) - if exist('implay', 'builtin') + if exist('implay', 'file') implay(obj.Mov, obj.FrameRate); else movie(obj.Mov, 1, obj.FrameRate); From 88b4b10ab4238bc0de2e722d9717bbebc845a222 Mon Sep 17 00:00:00 2001 From: Steven Roberts Date: Fri, 4 Oct 2024 19:01:16 -0700 Subject: [PATCH 3/5] Put Problem example output in correct position --- toolbox/+otp/Problem.m | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/toolbox/+otp/Problem.m b/toolbox/+otp/Problem.m index bd4cd3de..d40e022b 100644 --- a/toolbox/+otp/Problem.m +++ b/toolbox/+otp/Problem.m @@ -66,10 +66,10 @@ % ------- % >>> problem = otp.protherorobinson.presets.Canonical; % >>> problem.TimeSpan - % >>> problem.TimeSpan(2) = 100; % ans = % % 0 15 + % >>> problem.TimeSpan(2) = 100; TimeSpan %MATLAB ONLY: (1,2) {otp.utils.validation.mustBeNumerical} % A mutable column vector for the initial conditions $y(t_0) = y_0$. @@ -78,11 +78,11 @@ % ------- % >>> problem = otp.lotkavolterra.presets.Canonical; % >>> problem.Y0 - % >>> problem.Y0 = [3; 4]; % ans = % % 1 % 2 + % >>> problem.Y0 = [3; 4]; Y0 %MATLAB ONLY: (:,1) {otp.utils.validation.mustBeNumerical} % A mutable :class:`otp.Parameters` representing the problem parameters $p$. @@ -91,8 +91,8 @@ % ------- % >>> problem = otp.nbody.presets.Canonical; % >>> problem.Parameters.GravitationalConstant - % >>> problem.Parameters.Masses(1) = 2; % ans = 1 + % >>> problem.Parameters.Masses(1) = 2; % % See Also % -------- @@ -105,9 +105,9 @@ % ------- % >>> problem = otp.cusp.presets.Canonical; % >>> problem.NumVars + % ans = 96 % >>> problem.Y0 = ones(30, 1); % >>> problem.NumVars - % ans = 96 % ans = 30 NumVars end From 4a0916d04e772f33bd095069222b22301df2b771 Mon Sep 17 00:00:00 2001 From: Steven Roberts Date: Fri, 4 Oct 2024 19:39:42 -0700 Subject: [PATCH 4/5] Problem doc fixes --- toolbox/+otp/Problem.m | 57 ++++++++++++++++++++++-------------------- toolbox/+otp/RHS.m | 2 +- 2 files changed, 31 insertions(+), 28 deletions(-) diff --git a/toolbox/+otp/Problem.m b/toolbox/+otp/Problem.m index d40e022b..a5fdf663 100644 --- a/toolbox/+otp/Problem.m +++ b/toolbox/+otp/Problem.m @@ -187,20 +187,21 @@ % % Warning % ------- - % A plot might appear jagged if there are too few times included in the solution. In MATLAB, + % A plot can appear jagged if the solution time points are too coarse. In MATLAB, % `deval `_ can be used to evaluate the solution at a % user-specified set of points. % % Parameters % ---------- % sol : struct or numeric(1, :) or numeric(:, 1) - % A solution structure with properties ``x`` and ``y``. The length of one dimension of ``y`` must match - % the number of times in ``x``. If ``y`` is square, each column should be a state and each row the - % trajectory of a single variable. Alternatively, this argument can be a vector of times. + % A solution structure with properties ``x`` and ``y``. The length of one dimension of the matrix ``y`` + % must match the length (number of time points) of the vector ``x``. If ``y`` is square, each column + % should be a state and each row the trajectory of a single variable. Alternatively, this argument can be + % a vector of time points. % y : numeric(:, :), optional - % If the first argument is a vector of times, this should be the corresponding matrix of solutions. The - % length of one dimension must match the number of times. If square, each row should be a state and each - % column the trajectory of a single variable. + % If the first argument is a vector of time points, this should be the corresponding matrix of solutions. + % The length of one dimension must match the length (number of time points) of ``x``. If square, each row + % should be a state and each column the trajectory of a single variable. % varargin % A variable number of name-value pairs. Accepted names include % @@ -261,19 +262,20 @@ % % Warning % ------- - % A plot might appear jagged if there are too few times included in the solution. In MATLAB, `deval`_ can be - % used to evaluate the solution at a user-specified set of points. + % A plot can appear jagged if the solution time points are too coarse. In MATLAB, `deval`_ can be used to + % evaluate the solution at a user-specified set of points. % % Parameters % ---------- % sol : struct or numeric(1, :) or numeric(:, 1) - % A solution structure with properties ``x`` and ``y``. The length of one dimension of ``y`` must match - % the number of times in ``x``. If ``y`` is square, each column should be a state and each row the - % trajectory of a single variable. Alternatively, this argument can be a vector of times. + % A solution structure with properties ``x`` and ``y``. The length of one dimension of the matrix ``y`` + % must match the length (number of time points) of the vector ``x``. If ``y`` is square, each column + % should be a state and each row the trajectory of a single variable. Alternatively, this argument can be + % a vector of time points. % y : numeric(:, :), optional - % If the first argument is a vector of times, this should be the corresponding matrix of solutions. The - % length of one dimension must match the number of times. If square, each row should be a state and each - % column the trajectory of a single variable. + % If the first argument is a vector of time points, this should be the corresponding matrix of solutions. + % The length of one dimension must match the length (number of time points) of ``x``. If square, each row + % should be a state and each column the trajectory of a single variable. % varargin % A variable number of name-value pairs. Accepted names include those supported by % :func:`otp.Problem.plot` as well as the following: @@ -314,25 +316,26 @@ % % Warning % ------- - % A movie might appear jerky if there are too few times included in the solution. In MATLAB, `deval`_ can be - % used to evaluate the solution at a user-specified set of points. + % A movie can appear jerky if the solution time points are too coarse. In MATLAB, `deval`_ can be used to + % evaluate the solution at a user-specified set of points. % % Parameters % ---------- % sol : struct or numeric(1, :) or numeric(:, 1) - % A solution structure with properties ``x`` and ``y``. The length of one dimension of ``y`` must match - % the number of times in ``x``. If ``y`` is square, each column should be a state and each row the - % trajectory of a single variable. Alternatively, this argument can be a vector of times. + % A solution structure with properties ``x`` and ``y``. The length of one dimension of the matrix ``y`` + % must match the length (number of time points) of the vector ``x``. If ``y`` is square, each column + % should be a state and each row the trajectory of a single variable. Alternatively, this argument can be + % a vector of time points. % y : numeric(:, :), optional - % If the first argument is a vector of times, this should be the corresponding matrix of solutions. The - % length of one dimension must match the number of times. If square, each row should be a state and each - % column the trajectory of a single variable. + % If the first argument is a vector of time points, this should be the corresponding matrix of solutions. + % The length of one dimension must match the length (number of time points) of ``x``. If square, each row + % should be a state and each column the trajectory of a single variable. % varargin % A variable number of name-value pairs. Accepted names include those supported by % :func:`otp.Problem.plot` as well as the following: % % - ``Save`` – If ``false`` (default) the movie will play but the frames will not be save for playback. - % If ``true`` the frames in memory. If a string or a + % If ``true`` the frames will be saved in memory for playback. If a string or a % `VideoWriter `__, the movie will be % written to a file. % - ``FrameRate`` – The number of frames per second when playing the movie. @@ -340,8 +343,8 @@ % accommodate. % - ``Size`` – An array ``[width, height]`` for the figure size in pixels. % - ``Smooth`` – If ``true`` (default), solution states are uniformly sampled for frames based on the - % corresponding solution times. If ``false`` solution states are sampled uniformly by time step index. - % This may lead to the inconsistent playback speeds. + % corresponding solution times. If ``false`` solution states are sampled uniformly by time step index + % which may lead to the inconsistent playback speeds. % % Returns % ------- @@ -628,7 +631,7 @@ function onSettingsChanged(obj) end function t = parseTime(~, t) - % Puts a vector of times into a consistent row vector form + % Puts a vector of time points into a consistent row vector form if ~isvector(t) || isempty(t) || ~otp.utils.validation.isNumerical(t) error('OTP:invalidSolution', 'The times must be a nonempty vector of numbers'); diff --git a/toolbox/+otp/RHS.m b/toolbox/+otp/RHS.m index 6f1a6337..1aa97076 100644 --- a/toolbox/+otp/RHS.m +++ b/toolbox/+otp/RHS.m @@ -48,7 +48,7 @@ % t : numeric % The time at which $f$ is evaluated. % y : numeric(:, 1) or numeric(:, :) - % The state at which $f$ is evaluated. If :attr:`Vectorized` in ``'on'``, it can be a matrix where each + % The state at which $f$ is evaluated. If :attr:`Vectorized` is ``'on'``, it can be a matrix where each % column is a state. % % Returns From 2d4dd34ea9a142c41d71d2ac5a95b696ab268a4b Mon Sep 17 00:00:00 2001 From: Steven Roberts Date: Fri, 4 Oct 2024 21:03:38 -0700 Subject: [PATCH 5/5] Add RHS.odeset example --- toolbox/+otp/RHS.m | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/toolbox/+otp/RHS.m b/toolbox/+otp/RHS.m index 1aa97076..b8855f57 100644 --- a/toolbox/+otp/RHS.m +++ b/toolbox/+otp/RHS.m @@ -400,6 +400,12 @@ % opts : struct % An options structure containing the subset of properties compatible with built-in ODE solvers. % + % Example + % ------- + % >>> problem = otp.oregonator.presets.Canonical; + % >>> opts = problem.RHS.odeset('AbsTol', 1e-12, 'RelTol', 1e-4); + % >>> sol = ode23s(problem.RHS.F, problem.TimeSpan, problem.Y0, opts); + % % See Also % -------- %