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:

- 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")
```

It seems to have plotted the intermediate results before projection. With `save_everystep=true`

, it works as expected. Did I miss something?

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