# DiffEqOperators and gradients

**URL:** <https://discourse.julialang.org/t/diffeqoperators-and-gradients/70121>\
**Category:** Numerics\
**Created:** [October 20, 2021, 11:39pm UTC](https://discourse.julialang.org/t/diffeqoperators-and-gradients/70121 "2021-10-20T23:39:34Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![salazardetroya](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/salazardetroya/32/30093_2.png) [@salazardetroya](https://discourse.julialang.org/u/salazardetroya)\
**Post date:** [October 20, 2021, 11:39pm UTC](https://discourse.julialang.org/t/diffeqoperators-and-gradients/70121/1 "2021-10-20T23:39:34Z")

</div>

I am trying to solve the following problem

\begin{aligned} \nabla \cdot \nabla \phi\_{1} &= \frac{\partial\left(\phi\_{1}-\phi\_{2}\right)}{\partial t}\\ \nabla \cdot \nabla \phi\_{2} &= -\tau p\_1p\_2\frac{\partial\left(\phi\_{1}-\phi\_{2}\right)}{\partial t} \end{aligned}

in a one dimensional x\in[0, 1] domain, where p\_1 and p\_2 are parameters. The boundary conditions are

\begin{aligned} \phi\_1(0, t)& = t \\ \nabla\phi\_1(1, t) &= 0.0 \\ \phi\_2(1, t) &= 0.0 \\ \nabla\phi\_2(0, t) &= 0.0 \end{aligned}

I was successful in implementing the forward solver using the `DAEProblem` class.

```julia
using DiffEqOperators, OrdinaryDiffEq, Plots, DifferentialEquations, Zygote

nknots = 10
h = 1.0 / (nknots + 1)
knots = range(h, step=h, length=nknots)
ord_deriv = 2
ord_approx = 2

const Δ = CenteredDifference(ord_deriv, ord_approx, h, nknots)

τ = 0.1
ξ = 0.5
t0 = 0.0
tlimit = 1.0 / ξ
t1 = 2.0 * tlimit
vlimit = 1.0
p = [1.0, 2.0]

function heat!(res, du, u, p, t)
    rbc_1 = RobinBC((1.0, 0.0, t), (0.0, 1.0, 0.0), h, 1);
    rbc_2 = RobinBC((0.0, 1.0, 0.0), (1.0, 0.0, 0.0), h, 1);
    res[:, 1] = du[:, 1] - du[:, 2] - 1.0 * Δ * rbc_1 * u[:, 1]
    res[:, 2] = τ * p[1] * p[2] * (du[:, 1] - du[:, 2]) + 1.0 * Δ * rbc_2 * u[:, 2]
end

u0 = zeros(nknots, 2)
diff_var = trues(nknots, 2)
du0 = zeros(nknots, 2)
prob = DAEProblem{true}(heat!, du0, u0, (t0, t1)) # , differential_vars=diff_var)
sol = solve(prob, DABDF2(), abstol=1e-3, reltol=1e-3)

```

Now I’d like to get gradients with respect to p\_1 and p\_2, with the following lines:

```julia
function sum_of_solution(u0, p)
    _prob = remake(prob, u0=u0, p=p)
    sum(solve(_prob, DABDF2(), rtol=1e-6, atol=1e-6, saveat=0.1, sensealg=ReverseDiffAdjoint()))
end
du01, dp1 = Zygote.gradient(sum_of_solution, u0, p)

```

but I am getting the error:

```julia
ERROR: LoadError: MethodError: no method matching DiffEqOperators.BoundaryPaddedVector(::ReverseDiff.TrackedReal{Float64, Float64, Nothing}, ::ReverseDiff.TrackedReal{Float64, Float64, Nothing}, ::Vector{ReverseDiff.TrackedReal{Float64, Float64, ReverseDiff.TrackedArray{Float64, Float64, 2, Matrix{Float64}, Matrix{Float64}}}})
Closest candidates are:
  DiffEqOperators.BoundaryPaddedVector(::T, ::T, ::T2) where {T, T2<:AbstractVector{T}} at /Users/salazardetro1/.julia/packages/DiffEqOperators/FnmOt/src/boundary_padded_arrays.jl:14

```

I believe that my current equations might be solved using another method (maybe an ODEProblem with a mass matrix?), but my next step will be to include a nonlinear term multiplying the time derivatives as in

\begin{aligned} \nabla \cdot \nabla \phi\_{1} &= c(\phi\_1, \phi\_2) \frac{\partial\left(\phi\_{1}-\phi\_{2}\right)}{\partial t}\\ \nabla \cdot \nabla \phi\_{2} &= -\tau p\_1p\_2c(\phi\_1, \phi\_2)\frac{\partial\left(\phi\_{1}-\phi\_{2}\right)}{\partial t} \end{aligned}

and I think that I can only solve those problems with a `DAEProblem`. With that said, how can I compute the gradients? Thanks.

---

<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:** [October 21, 2021, 2:19pm UTC](https://discourse.julialang.org/t/diffeqoperators-and-gradients/70121/2 "2021-10-21T14:19:29Z")

</div>

> [@salazardetroya](#):
>
> ```julia
> ERROR: LoadError: MethodError: no method matching DiffEqOperators.BoundaryPaddedVector(::ReverseDiff.TrackedReal{Float64, Float64, Nothing}, ::ReverseDiff.TrackedReal{Float64, Float64, Nothing}, ::Vector{ReverseDiff.TrackedReal{Float64, Float64, ReverseDiff.TrackedArray{Float64, Float64, 2, Matrix{Float64}, Matrix{Float64}}}})
> 
> ```

That’s not an error with DAEProblem gradients, but differentiation of DiffEqOperators.jl. I assume this would happen on ODEProblems as well. It would be good to isolate this down to a simple example and open an issue on DiffEqOperators.jl

---

<div class="post-metadata">

**Author:** ![salazardetroya](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/salazardetroya/32/30093_2.png) [@salazardetroya](https://discourse.julialang.org/u/salazardetroya)\
**Post date:** [October 21, 2021, 4:00pm UTC](https://discourse.julialang.org/t/diffeqoperators-and-gradients/70121/3 "2021-10-21T16:00:05Z")

</div>

FWIW, there is an issue open already: [https://github.com/SciML/DiffEqOperators.jl/issues/333](https://github.com/SciML/DiffEqOperators.jl/issues/333) I will leave it here for future reference.