# DifferentialEquations.jl inplace implicit ODE problem

**URL:** <https://discourse.julialang.org/t/differentialequations-jl-inplace-implicit-ode-problem/37495>\
**Category:** General Usage\
**Tags:** package\
**Created:** [April 13, 2020, 9:16am UTC](https://discourse.julialang.org/t/differentialequations-jl-inplace-implicit-ode-problem/37495 "2020-04-13T09:16:53Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![santiagobadia](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/santiagobadia/32/7792_2.png) [@santiagobadia](https://discourse.julialang.org/u/santiagobadia)\
**Post date:** [April 13, 2020, 9:16am UTC](https://discourse.julialang.org/t/differentialequations-jl-inplace-implicit-ode-problem/37495/1 "2020-04-13T09:16:53Z")

</div>

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

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

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

```

---

<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:** [April 13, 2020, 11:07am UTC](https://discourse.julialang.org/t/differentialequations-jl-inplace-implicit-ode-problem/37495/2 "2020-04-13T11:07:47Z")

</div>

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

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.

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

```

> [@santiagobadia](#):
>
> And it seems we must provide an initial value for `du` that does not have mathematical sense in our problems.

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.

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

No. You can use a state-dependent mass matrix with an algorithm for ODEs, or you can use a DAEProblem with “implicit ODE” DAEs.

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [April 13, 2020, 1:06pm UTC](https://discourse.julialang.org/t/differentialequations-jl-inplace-implicit-ode-problem/37495/3 "2020-04-13T13:06:32Z")

</div>

`Gridap.jl` is nice! I would like to use it with [PseudoArcLengthContinuation.jl](https://github.com/rveltz/PseudoArcLengthContinuation.jl). I will see how you expose your jacobian in this thread

---

<div class="post-metadata">

**Author:** ![santiagobadia](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/santiagobadia/32/7792_2.png) [@santiagobadia](https://discourse.julialang.org/u/santiagobadia)\
**Post date:** [April 14, 2020, 11:29pm UTC](https://discourse.julialang.org/t/differentialequations-jl-inplace-implicit-ode-problem/37495/4 "2020-04-14T23:29:51Z")

</div>

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

In any case, following

[https://docs.sciml.ai/latest/features/performance\_overloads/#](https://docs.sciml.ai/latest/features/performance_overloads/#)

I am creating a DAE problem

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

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

---

<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:** [April 15, 2020, 5:12am UTC](https://discourse.julialang.org/t/differentialequations-jl-inplace-implicit-ode-problem/37495/5 "2020-04-15T05:12:24Z")

</div>

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

[https://github.com/SciML/OrdinaryDiffEq.jl/issues/1109](https://github.com/SciML/OrdinaryDiffEq.jl/issues/1109)

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