Difference between an ODE with manifold projection and a DAE

Specific Domains Numerics
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.

1 Like