# DifferentialEquations.jl inplace implicit ODE problem

I am trying to combine the ODE solvers in `DifferentialEquations.jl` with the PDE solvers in the `Gridap.jl` library.

It is mandatory for our applications to use the inplace version of the `DifferentialEquations.jl` package. This way, we can pre-allocate the required arrays in our `Gridap` solvers.

On the other hand, even when considering ODE systems, we want to use an implicit ODE expression of the solver, i.e., `A(du,u,p,t)=0`, since our mass matrices are not just identity matrices, and can be fairly complicated or even nonlinear.

Finally, what we would provide are a `residual(res,du,u,p,t)` and `jacobian(J,du,u,p,t)` functions.

As a result, I think that the way to go is to consider an IIP problem in `DifferentialEquations.jl`. But it seems that only the `DAE` constructor allows such syntax.

However, for implicit ODE systems I cannot use an `ODE` algorithm, e.g., Backward-Euler, even though it is acceptable for our problems when the mass matrix is non-singular (thus not a `DAE`). And it seems we must provide an initial value for `du` that does not have mathematical sense in our problems.

Summarizing, is there a way to create an `ODE` solver that takes `residual!(res,du,u,p,t)` and `jacobian!(J,du,u,p,t)` inplace functions, and allows one to use an `ODE` algorithm?

I have written the following silly linear example with what I would like to have

``````tspan = (0.0,1.0)
u0 = zeros(1)
u0[1] = 1/2

function residual(res,du,u,p,t)
res[1] = - 1.01u[1] + du[1]
end

function jacobian(jac,du,u,p,gamma,t)
jac[1,1] = gamma - 1.01
nothing
end

f_iip = ODEFunction{true}(residual;jac=jacobian)
prob_iip = ODEProblem{true}(f_iip,u0,tspan)
sol_iip = solve(prob_iip, ImplicitEuler(), reltol=1e-8, abstol=1e-8)
``````

but it does not work, because it seems the library is looking for an `ODE`-style inplace version of `residual`, i.e., `residual(du,u,p,t)`, based on the error:

``````MethodError: no method matching residual(::Array{Float64,1}, ::Array{Float64,1}, ::DiffEqBase.NullParameters, ::Float64)
Closest candidates are:
residual(::Any, ::Any, ::Any, ::Any, !Matched::Any) at /.../DifEqWrappers.jl:20
...
``````

If you can write it as a matrix, you could use the state-dependent mass matrix part of the interface, and it can be more efficient.

``````function update_func(A,u,p,t)
A[1,1] = cos(t)
A[2,1] = sin(t)*u[1]
A[3,1] = t^2
A[1,2] = cos(t)*sin(t)
A[2,2] = cos(t)^2 + u[3]
A[3,2] = sin(t)*u[2]
A[1,3] = sin(t)
A[2,3] = t^2
A[3,3] = t*cos(t) + 1
end
dependent_M = DiffEqArrayOperator(ones(3,3),update_func=update_func)
prob = ODEProblem(ODEFunction(f, mass_matrix=mm_A), u0, tspan))
``````

Itâ€™s just used in the consistent initialization algorithms. You can set it to all zero if you want. The interface for `DAEProblem` should probably make that `du0` optional in the future, but for now itâ€™s required but sometimes ignored by solvers/initialization algorithms.

No. You can use a state-dependent mass matrix with an algorithm for ODEs, or you can use a DAEProblem with â€śimplicit ODEâ€ť DAEs.

`Gridap.jl` is nice! I would like to use it with PseudoArcLengthContinuation.jl. I will see how you expose your jacobian in this thread

1 Like

Thanks for the quick reply, I will consider this option.

In any case, following

I am creating a DAE problem

``````tspan = (0.0,1.0)
u0 = zeros(1)
u0[1] = 1/2

function residual(res,du,u,p,t)
res[1] = - 1.01u[1] + du[1]
end

function jacobian(jac,du,u,p,gamma,t)
jac[1,1] = gamma - 1.01
nothing
end

f_iip = DAEFunction{true}(residual;jac=jacobian)
prob_iip = DAEProblem{true}(f_iip,u0,u0,tspan)
sol_iip = solve(prob_iip, DImplicitEuler(), reltol=1e-8, abstol=1e-8)
``````

but it returns an error,

``````type ODEIntegrator has no field duprev
``````

which is weird, since I would expect a `DAEIntegrator` being created.

It would be great if you could tell me what I am doing wrong.

Itâ€™s us thatâ€™s doing something wrong. The new DAE solvers donâ€™t support user Jacobians quite yet:

Bare with us: these new integrators are missing the polish right now, so Iâ€™d still recommend using `IDA()` for production for the time being (though these are improving daily!).

1 Like