# Difference between an ODE with manifold projection and a DAE

**URL:** <https://discourse.julialang.org/t/difference-between-an-ode-with-manifold-projection-and-a-dae/45136>\
**Category:** Numerics\
**Created:** [August 18, 2020, 1:20am UTC](https://discourse.julialang.org/t/difference-between-an-ode-with-manifold-projection-and-a-dae/45136 "2020-08-18T01:20:35Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![cuihantao](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cuihantao/32/8783_2.png) [@cuihantao](https://discourse.julialang.org/u/cuihantao)\
**Post date:** [August 18, 2020, 1:20am UTC](https://discourse.julialang.org/t/difference-between-an-ode-with-manifold-projection-and-a-dae/45136/1 "2020-08-18T01:20:35Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [August 18, 2020, 1:38am UTC](https://discourse.julialang.org/t/difference-between-an-ode-with-manifold-projection-and-a-dae/45136/2 "2020-08-18T01:38:47Z")

</div>

This should explain most of it:

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

They are related but not the same.

---

<div class="post-metadata">

**Author:** ![cuihantao](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cuihantao/32/8783_2.png) [@cuihantao](https://discourse.julialang.org/u/cuihantao)\
**Post date:** [August 18, 2020, 2:13am UTC](https://discourse.julialang.org/t/difference-between-an-ode-with-manifold-projection-and-a-dae/45136/3 "2020-08-18T02:13:43Z")

</div>

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:

```julia
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](https://global.discourse-cdn.com/julialang/original/3X/3/5/350f7fef4ecb15c30fb96ff830745c433545589a.png)

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:

```julia
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:

```julia
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.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [August 18, 2020, 1:08pm UTC](https://discourse.julialang.org/t/difference-between-an-ode-with-manifold-projection-and-a-dae/45136/4 "2020-08-18T13:08:40Z")

</div>

> [@cuihantao](#):
>
> It seems to have plotted the intermediate results before projection. With `save_everystep=true` , it works as expected. Did I miss something?

> <https://github.com/SciML/DiffEqCallbacks.jl/blob/master/src/manifold.jl#L72>

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

```julia
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](https://global.discourse-cdn.com/julialang/original/3X/3/2/3276cc9a94570ad2181847329fd10da49b414b98.png)

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

> [@cuihantao](#):
>
> 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:

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.