Use band-storage jacobian with CVODE_BDF in OrdinaryDiffEq

General Usage
,
1

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!

2

Use the in-place format.

3

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[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
4

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
5

Thanks! Good to know. I opened an issue.