# Difference between an ODE with manifold projection and a DAE

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.

1 Like

This should explain most of it:

They are related but not the same.

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 = u
du = -u
end
prob = ODEProblem(f,u0,(0.0,100.0))

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

cb = ManifoldProjection(g)

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

@test sol[end]^2 + sol[end]^2 ≈ 2
plot(sol,vars=(1,2),label="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) / d(u) == 0`. Here’s the code that won’t work:
``````u0 = ones(3)
function f(du,u,p,t)
du = u
du = -u
du = u^2 + u^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 = u
du = -u
du = u^2 + u^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
J[3, 1] = 2u
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.

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 = u
du = -u
end
prob = ODEProblem(f,u0,(0.0,100.0))

function g(resid,u,p,t)
resid = u^2 + u^2 - 2
resid = 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^2 + x^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: 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.

1 Like