Constrained pendulum problem #36
Conversation
For use in half-explicit methods and for data assimilation.
| % Ascher, Uri. "On symmetric schemes and differential-algebraic equations." | ||
| % SIAM journal on scientific and statistical computing 10.5 (1989): 937-949. |
There was a problem hiding this comment.
Is this the right source? I don't see anything on a pendulum.
There was a problem hiding this comment.
Hairer, E., Roche, M., Lubich, C. (1989). Description of differential-algebraic problems. In: The Numerical Solution of Differential-Algebraic Systems by Runge-Kutta Methods. Lecture Notes in Mathematics, vol 1409. Springer, Berlin, Heidelberg. https://doi.org/10.1007/BFb0093948
That is the pendulum problem source.
The Ascher book has a nice formulation as a Hessenberg indeex 2 DAE( linear algebraic variables) we use.
There was a problem hiding this comment.
Yeah the citation is wrong.
| %Gravity defines the magnitude of acceleration towards the ground | ||
| Gravity %MATLAB ONLY: (1,1) {mustBeNonnegative} | ||
| %Mass defines the mass of the pendulum | ||
| Mass %MATLAB ONLY: (1,1) {mustBeNonnegative} |
There was a problem hiding this comment.
Sharper validation: {mustBeReal, mustBeFinite, mustBePositive}
There was a problem hiding this comment.
| @@ -0,0 +1,35 @@ | |||
| classdef ConstrainedPendulumProblem < otp.Problem | |||
There was a problem hiding this comment.
Perhaps just PendulumDAEProblem.
There was a problem hiding this comment.
This class ought to override internalIndex2label
| %CONSTRAINEDPENDULUM PROBLEM This is a Hessenberg Index-2 DAE problem posed in terms of three constraints | ||
| % | ||
|
|
||
| properties |
There was a problem hiding this comment.
These should have SetAccess = private or protected
| initialconstraints = otp.constrainedpendulum.constraints([], obj.Y0, g, m, l, 0); | ||
| E0 = initialconstraints(3); | ||
|
|
||
| obj.RHS = otp.RHS(@(t, y) otp.constrainedpendulum.f(t, y, g, m, l, E0)); |
There was a problem hiding this comment.
This RHS will lead to confusion as a call to solve will not enforce constraints, right? Ideally, RHS should be in a form that a built-in solver can integrate, but we can let that slide for an index-2 DAE. I suggest the default RHS should use a mass matrix and an f that operates on the full state (diff and alg). Then there can be a RHSDifferential and RHSAlgebraic to get the individual pieces.
There was a problem hiding this comment.
Also it looks like 'Vectorized', 'on' can be added to RHS
| @@ -0,0 +1,27 @@ | |||
| function dc = constraintsjacobian(~, state, g, m, ~, ~) | |||
There was a problem hiding this comment.
What Jacobian is this exactly? For the DAE
y' = f(y, z)
0 = g(y)
is it g_y or f_z*g_y
There was a problem hiding this comment.
in this formulation. y' = fhat(y) - ((g_y)^T)*z. T is transpose.
so g_y * f_z = g_y*((g_y)^T)
There was a problem hiding this comment.
this jacobian is actually just g_y
|
Long term, I would like to have a problem class for constrained mechanical systems (Solving ODEs II pg. 464), and the pendulum would simply be a preset. This would also unify the handling of reduced index forms. |
I agree longterm. Right now we need this problem for something. We might also add a constrained shallow water (2D) |
|
|
||
| y0 = [sqrt(2)/2; sqrt(2)/2; 0; 0]; | ||
|
|
||
| obj = [email protected](tspan, y0, params); |
There was a problem hiding this comment.
Will the solve functionality be disabled for this problem?
There was a problem hiding this comment.
That is an interesting question. In theory we should be able to add our own small method to OTP to handle these cases
There was a problem hiding this comment.
Good point. Here's one idea: Add a HighIndexDAE property to https://github.com/ComputationalScienceLaboratory/ODE-Test-Problems/blob/master/src/%2Botp/%2Butils/Solver.m which simply throws an error message that high-index DAEs are not supported by matlab solvers. Eventually when we have a constrained mechanical system problem, we will support RHS reduced to index-1 which can rely on built-in solvers.
| @@ -0,0 +1,27 @@ | |||
| function dc = constraintsjacobian(~, state, g, m, ~, ~) | |||
There was a problem hiding this comment.
Is this to be the standard for otp when getting the Jacobian for Half-Explicit methods? Or is this an internal functionality of the problem?
There was a problem hiding this comment.
if the problem is posed in terms of invariants, this should be standard. I might rename it to invariantJacobian
| @@ -0,0 +1,19 @@ | |||
| function dstate = fdifferential(~, state, g, m, ~, ~) | |||
There was a problem hiding this comment.
The convention for functions is camel case now: https://github.com/ComputationalScienceLaboratory/ODE-Test-Problems/blob/master/docs/CONTRIBUTING.md#functions
|
|
||
| % The right hand size in terms of x, y, x', y', and three control parameters z | ||
| obj.RHS = otp.RHS(@(t, y) otp.pendulumdae.f(t, y, g, m, l, E0), ... | ||
| 'Mass', @(t, y) otp.pendulumdae.mass(t, y, g, m, l, E0)); |
There was a problem hiding this comment.
Since the mass matrix is constant, it can be a matrix instead of a function handle. See
for exampleThere was a problem hiding this comment.
I was just copying the pendulum problem, which had @(t,y). I will change
|
Unless I am missing something, I think I have addressed all the issues. |
| Y0Differential | ||
| Z0Algebraic |
There was a problem hiding this comment.
I'm not opposed to having these, but if we do include this, I think all DAE problems should have these for consistency. Also, I think Y0Algebraic would be better than Z0Algebraic.
There was a problem hiding this comment.
Doesn't really differentiate the two. Algebraic implies there is no direct differential equation thus it is a variable obtained from differentiation of the constraint. I think it should stay.
There was a problem hiding this comment.
I think this naming is pretty clear. It's mostly motivated by established naming conventions. For partitioned problems, we use <property><partition name>. E.g.
RHSStiff and RHSNonstiff
jacobianDifferential and jacobianAlgebraic
Y0Position and Y0Velocity
This is nice for autocompletion as well.
There was a problem hiding this comment.
My bad I misinterpreted what was being proposed. I thought you wanted Y0Differential -> Y0Algebraic not Z0Algebraic -> Y0Algebraic. I agree on the convention you proposed.
There was a problem hiding this comment.
I was thinking of changing it to that already, glad to see we all think alike
| % Generate the constituent RHS for the differential part | ||
| obj.RHSDifferential = otp.RHS(@(t, y) otp.pendulumdae.fDifferential(t, y, g, m, l, E0), ... | ||
| 'JacobianVectorProduct', ... | ||
| @(t, y, v) otp.pendulumdae.jacobianDifferentialVectorProduct(t, y, g, m, l, E0, v), ... |
There was a problem hiding this comment.
jacobianVectorProductDifferential is the naming convention https://github.com/ComputationalScienceLaboratory/ODE-Test-Problems/blob/master/docs/CONTRIBUTING.md#functions-for-rhs
| 'Vectorized', 'on'); | ||
|
|
||
| % Generate the constituent RHS for the algebraic part | ||
| obj.RHSAlgebraic = otp.RHS(@(t, y) otp.pendulumdae.invariants(t, y, g, m, l, E0), ... |
There was a problem hiding this comment.
The terms "algebraic" and "invariant" are used inconsistently. I think "algebraic" is the more standard DAE term, so I say we should name the functions
fAlgebraic
jacobianAlgebraic
...
Added some final pedantic comments |
src/+otp/+utils/+solvers/dae34.m
Outdated
| reltol = odeget(options, 'RelTol', 1e-3); | ||
| abstol = odeget(options, 'AbsTol', 1e-6); | ||
| J = odeget(options, 'Jacobian', @(t, y) jacapprox(f, t, y)); | ||
| M = odeget(options, 'Mass', speye(numel(y0))); |
There was a problem hiding this comment.
We should determine if the mass matrix needs to be sparse from the Jacobian. Ideally, we could call eye(n, 'like', J) but this is not supported in Octave.
| h = tend - tc; | ||
| end | ||
|
|
||
| stages = zeros(numel(yc), stagenum); |
There was a problem hiding this comment.
No need to reallocate every step
src/+otp/+utils/+solvers/dae34.m
Outdated
|
|
||
| np = inf; | ||
|
|
||
| ntol = 1e-6; |
There was a problem hiding this comment.
This should be based on AbsTol and RelTol
src/+otp/+utils/+solvers/dae34.m
Outdated
| nmaxits = 5; | ||
| its = 0; | ||
| while norm(np) > ntol || its < nmaxits | ||
| Mc = M(staget, yc + stagedy + gh*newtonk0); |
There was a problem hiding this comment.
If the mass matrix is state dependent, a Newton iteration would need dM/dy. We should parse the MStateDependence property to handle this better
src/+otp/+utils/+solvers/dae34.m
Outdated
| Mc = M(staget, yc + stagedy + gh*newtonk0); | ||
| g = Mc*newtonk0 - f(staget, yc + stagedy + gh*newtonk0); | ||
| H = Mc - gh*J(staget, yc + stagedy + gh*newtonk0); | ||
| [np, ~] = gmres(H, g); |
There was a problem hiding this comment.
A direct solver would be more standard. Unfortunately, Octave does not support decomposition, so lu is probably best.
There was a problem hiding this comment.
I found that the matrix is often very poorly scaled, so a direct solve was throwing tons of warnings.
src/+otp/+utils/+solvers/dae34.m
Outdated
| while norm(np) > ntol || its < nmaxits | ||
| Mc = M(staget, yc + stagedy + gh*newtonk0); | ||
| g = Mc*newtonk0 - f(staget, yc + stagedy + gh*newtonk0); | ||
| H = Mc - gh*J(staget, yc + stagedy + gh*newtonk0); |
There was a problem hiding this comment.
We can upgrade Newton iterations to algorithm on pg 119 of Solving ODEs II. This will evaluate the Jacobian only as needed.
src/+otp/+utils/+solvers/dae34.m
Outdated
|
|
||
| Mc = M(tc + h, yc); | ||
|
|
||
| err = rms((Mc*(yc-yhat))./sc); |
There was a problem hiding this comment.
rms not supported in Octave
There was a problem hiding this comment.
Not sure why error is multiplied by Mc
There was a problem hiding this comment.
I can do manual rms. Error is multiplied by Mc, as one way to do error is to simply look at the differential variables and ignore the algebraic.
|
|
||
| end | ||
|
|
||
| function J = jacapprox(f, t, y) |
There was a problem hiding this comment.
We'll need something a bit more robust or not support finite differences at all.
…lgorithm. The starting algorithm for FIRK can actually just use the previous stages. This speeds up computation and just works. For SDIRK this does not work, just for some Lobatto methods
src/+otp/+utils/+solvers/dae43.m
Outdated
| reltol = odeget(options, 'RelTol', 1e-3); | ||
| abstol = odeget(options, 'AbsTol', 1e-6); | ||
| J = odeget(options, 'Jacobian', []); | ||
| M = odeget(options, 'Mass', speye(numel(y0))); |
There was a problem hiding this comment.
Default value should be []
src/+otp/+utils/+solvers/dae43.m
Outdated
|
|
||
| % Compute initial step size | ||
| if isempty(h) | ||
| sc = abstol + reltol*abs(y0); |
There was a problem hiding this comment.
RelTol could be a vector. Should use .*
src/+otp/+utils/+solvers/dae43.m
Outdated
| step = 1; | ||
|
|
||
| Mc = M(tspan(1), y0); | ||
| Mfull = zeros(n*stagenum, n*stagenum, 'like', Mc); |
There was a problem hiding this comment.
Unfortunately, 'like' option is not supported in Octave
src/+otp/+utils/+solvers/dae43.m
Outdated
|
|
||
| end | ||
|
|
||
| function J = jacapprox(f, t, y) |
There was a problem hiding this comment.
This seems like a useful enough to warrant it's own publicly-accessible function. We can also look into more sophisticated and robust approaches.
| M = @(t, y) eye(numel(y0)); | ||
| end | ||
| elseif ~isa(M, 'function_handle') | ||
| M = @(t, y) M; |
There was a problem hiding this comment.
We also need to support mass matrices that are only a function of time, e.g., M(t)
|
All resolved. |
|
Should I remove the other DAE solvers (SDIRK and ESDIRK)? |
For use in half-explicit methods and for data assimilation.