Difference between an ODE with manifold projection and a DAE

Specific Domains Numerics
1

Hi,

I’d like to know the difference between an ODE set with manifold projection callback and a DAE (formulated as ODE with a singular mass matrix). It seems to me that the former one separates the algebraic equations from the latter and put that into a callback, which is evaluated after integrating each step.

Could anyone explain the fundamental difference between the two (or say, they have no relationship at all)? Thanks in advance.

2

This should explain most of it:

https://scicomp.stackexchange.com/questions/34129/conserving-energy-in-physics-simulation-with-imperfect-numerical-solver/34131#34131

They are related but not the same.

3

Chris,

Thank you for your reply. I’ve got a rough idea about that - both would pretty much work, but the applicability depends on the problem.

I have two follow-up questions:

  1. On the ManifoldProjection documentation, the example code, which has save_everystep=true, plots a different figure.

Code:

u0 = ones(2)
function f(du,u,p,t)
  du[1] = u[2]
  du[2] = -u[1]
end
prob = ODEProblem(f,u0,(0.0,100.0))

function g(resid,u,p,t)
  resid[1] = u[2]^2 + u[1]^2 - 2
  resid[2] = 0
end

cb = ManifoldProjection(g)

sol = solve(prob,Vern7(),save_everystep=false,callback=cb)

@test sol[end][1]^2 + sol[end][2]^2 ≈ 2
plot(sol,vars=(1,2),label="Projection")

manifold-projection

It seems to have plotted the intermediate results before projection. With save_everystep=true, it works as expected. Did I miss something?

  1. I tried to rewrite the above example in the mass-matrix formulation, but it turns out the Jacobian is singular. Namely, d(du[3]) / d(u[3]) == 0. Here’s the code that won’t work:
u0 = ones(3)
function f(du,u,p,t)
  du[1] = u[2]
  du[2] = -u[1]
  du[3] = u[2]^2 + u[1]^2 - 2
end

mm = [1 0 0; 0 1 0; 0 0 0]

func = ODEFunction(f, mass_matrix=mm)
prob = ODEProblem(func, u0,(0.0,100.0))

sol = solve(prob,Rodas5(),save_everystep=false)

It raises SingularException(3).

What I have been doing is to manually provide a Jacobian with a small number at J[3, 3]. It works fine but will likely introduce errors. Here’s the code with a manual Jacobian callback:

u0 = ones(3)
function f(du,u,p,t)
  du[1] = u[2]
  du[2] = -u[1]
  du[3] = u[2]^2 + u[1]^2 - 2
end

mm = [1 0 0; 0 1 0; 0 0 0]

function jac(J, u, p, t)
    J[1, 2] = 1.
    J[2, 1] = -1.
    J[3, 2] = 2u[2]
    J[3, 1] = 2u[1]
    J[3, 3] = 1e-12
end

func = ODEFunction(f, mass_matrix=mm, jac=jac)

prob = ODEProblem(func, u0,(0.0,100.0))

sol = solve(prob,Rodas5())

It seems feasible to formulate it as a fully implicit problem, but then it will lose compatibility with lots of the mass-matrix solvers. Do you know if there are other commonly used techniques or tricks (or an existing option) to workaround such singularity?

Thank you.

4

It only saves the post-projected points by defaults, and you can turn that off.

using OrdinaryDiffEq, DiffEqCallbacks, Test, Plots

u0 = ones(2)
function f(du,u,p,t)
  du[1] = u[2]
  du[2] = -u[1]
end
prob = ODEProblem(f,u0,(0.0,100.0))

function g(resid,u,p,t)
  resid[1] = u[2]^2 + u[1]^2 - 2
  resid[2] = 0
end

cb = ManifoldProjection(g,nlopts=Dict(:ftol=>1e-16,:xtol=>1e-16))
sol = solve(prob,Vern7(),save_everystep=false,callback=cb)
@test all(x->x[1]^2 + x[2]^2 ≈ 2,sol.u) # passes
plot(sol,vars=(1,2),label="Projection")

You can see that all of the saved points satisfy the criteria. You just have to plot it another way:

cool

to see that your mind is deceiving you an everything is on the circle. This is actually a really cool example haha.

For the DAE, you end up with an Index-2 DAE in that form, so you need to do an index reduction to get a 3 Index-1 DAE by hand, which is not trivial but not too hard.