 # Use band-storage jacobian with CVODE_BDF in OrdinaryDiffEq

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

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!

Use the in-place format.

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.

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 # flux into the bottom layer
dudt = (1.0/h^2.0)*(-u + u) + 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


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.

1 Like

Thanks! Good to know. I opened an issue.