# Use band-storage jacobian with CVODE\_BDF in OrdinaryDiffEq

**URL:** <https://discourse.julialang.org/t/use-band-storage-jacobian-with-cvode-bdf-in-ordinarydiffeq/58068>\
**Category:** General Usage\
**Tags:** diffeq, ode\
**Created:** [March 27, 2021, 7:06am UTC](https://discourse.julialang.org/t/use-band-storage-jacobian-with-cvode-bdf-in-ordinarydiffeq/58068 "2021-03-27T07:06:39Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![Nicholaswogan](https://avatars.discourse-cdn.com/v4/letter/n/f1d935/32.png) [@Nicholaswogan](https://discourse.julialang.org/u/Nicholaswogan)\
**Post date:** [March 27, 2021, 7:06am UTC](https://discourse.julialang.org/t/use-band-storage-jacobian-with-cvode-bdf-in-ordinarydiffeq/58068/1 "2021-03-27T07:06:39Z")

</div>

I am trying to solve a large system of ODEs. The jacobian of these ODEs is a band matrix. Is it possible to feed the CVODE\_BDF algorithm, from Sundials.jl, a function for the jacobian which returns the jacobian in band storage (not a dense matrix). The answer should look something like the following

```julia
function right_hand_side(u,p,t)
    ...
    return dudt
end
function jacobian(u,p,t)
    ...
    return J # must be band matrix of some kind
end

tspan = [0.0,1.0]
func = ODEFunction(right_hand_side, jac=jacobian)
prob = ODEProblem(func,u0,tspan)
sol = solve(prob,CVODE_BDF(linear_solver=:Band,jac_upper=10,jac_lower=10));

```

I have tried making `J` a `BandedMatrix` (from BandedMatrices.jl). This does not work with CVODE\_BDF. I have also tried making `J` a band storage matrix, where each row is a diagonal. This also does not work. Thanks for any advice!

---

<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:** [March 27, 2021, 12:04pm UTC](https://discourse.julialang.org/t/use-band-storage-jacobian-with-cvode-bdf-in-ordinarydiffeq/58068/2 "2021-03-27T12:04:43Z")

</div>

Use the in-place format.

---

<div class="post-metadata">

**Author:** ![Nicholaswogan](https://avatars.discourse-cdn.com/v4/letter/n/f1d935/32.png) [@Nicholaswogan](https://discourse.julialang.org/u/Nicholaswogan)\
**Post date:** [March 27, 2021, 4:27pm UTC](https://discourse.julialang.org/t/use-band-storage-jacobian-with-cvode-bdf-in-ordinarydiffeq/58068/3 "2021-03-27T16:27:51Z")

</div>

I have tried in place format, but I get similar errors. Below is an example solving the diffusion equation in 1 dimension: \frac{du}{dt} = \frac{d^2u}{dz^2}. Replacing spatial derivatives with center differences results in ODEs with a tridiagonal jacobian. In this case, my specified jacobian is not necessary because Sundials can calculate the jacobian with finite differences - this is just for demonstrations sake. How should i format the band-storage jacobian `J`, and which options should be selected, so that Sundials can use it.

```julia
using OrdinaryDiffEq
using Sundials
using ForwardDiff
using BandedMatrices

struct parameters{T<:AbstractFloat}
    v_dep::T # lower boundary condition
    PhiT::T # upper flux boundary condition
end

function rhs(dudt,u,p::parameters,t) # right-hand-side of ODEs
    n = length(u)
    h = 1.0/n # spatial grid size
    for i in 2:(n-1)
        dudt[i] = (1/h^2.0)*(u[i-1] - 2.0*u[i] + u[i+1])
    end

    # lower boundary
    PhiG = -p.v_dep*u[1] # flux into the bottom layer
    dudt[1] = (1.0/h^2.0)*(-u[1] + u[2]) + PhiG/h

    # upper boundary
    dudt[end] = (1.0/h^2.0)*(u[end-1] - u[end]) + p.PhiT/h
    nothing
end;

function jacobian(J,u,p::parameters,t)
    dudt = Array{Float64,1}(undef,length(u))
    rhs1 = (dudt,uu) -> rhs(dudt,uu,p,0.0)
    ForwardDiff.jacobian!(J,rhs1,dudt,u0)
    nothing
end

n = 1000 # number of ODEs.
u0 = zeros(Float64,n) # initial condition
params = parameters(10.0,1.0) # boundary conditions
tspan = [0.0,10.0]

# Solve with NO SPECIFIED JACOBIAN
prob1 = ODEProblem(rhs,u0,tspan,params)
sol1 = solve(prob1,CVODE_BDF(linear_solver=:Band,jac_upper=1,jac_lower=1)); # this works

# Solving WITH specified jacobian
J = BandedMatrix(zeros(Float64,(n,n)), (1,1))
jacobian(J,u0,params,0.0)
f = ODEFunction(rhs, jac=jacobian,jac_prototype=J)
prob2 = ODEProblem(f,u0,tspan,params)
sol2 = solve(prob2,Rodas5()); # this works
sol3 = solve(prob2,CVODE_BDF(linear_solver=:Band,jac_upper=1,jac_lower=1)) # Method error. CVODE can't accept BandedMatrix

```

---

<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:** [March 27, 2021, 4:53pm UTC](https://discourse.julialang.org/t/use-band-storage-jacobian-with-cvode-bdf-in-ordinarydiffeq/58068/4 "2021-03-27T16:53:55Z")

</div>

Oh yes, I don’t know if we wrapped the hook for analytical Jacobians of banded matrices. Open an issue on Sundials.jl about that.

---

<div class="post-metadata">

**Author:** ![Nicholaswogan](https://avatars.discourse-cdn.com/v4/letter/n/f1d935/32.png) [@Nicholaswogan](https://discourse.julialang.org/u/Nicholaswogan)\
**Post date:** [March 28, 2021, 2:07am UTC](https://discourse.julialang.org/t/use-band-storage-jacobian-with-cvode-bdf-in-ordinarydiffeq/58068/5 "2021-03-28T02:07:21Z")

</div>

Thanks! Good to know. I opened [an issue.](https://github.com/SciML/Sundials.jl/issues/302)