Dear Julia Users,
I'm trying to speed up some code that employs passing functions as
arguments. One part of the code solves an ODE; if I use CVODE from
Sundials, and rewrite the function to accept callable types, I get a
ReadOnlyMemoryError - as I don't know the Sundials API, can someone help me
with where I'm going wrong? Code below.
Best
Simon
**
using Sundials
function cvode{f}(::Type{f}, y0::Vector{Float64}, t::Vector{Float64};
reltol::Float64=1e-4, abstol::Float64=1e-6)
neq = length(y0)
mem = Sundials.CVodeCreate(Sundials.CV_BDF, Sundials.CV_NEWTON)
flag = Sundials.CVodeInit(mem, cfunction(Sundials.cvodefun, Int32,
(Sundials.realtype, Sundials.N_Vector, Sundials.N_Vector, Ref{Function})),
t[1], Sundials.nvector(y0))
flag = Sundials.CVodeSetUserData(mem, f)
flag = Sundials.CVodeSStolerances(mem, reltol, abstol)
flag = Sundials.CVDense(mem, neq)
yres = zeros(length(t), length(y0))
yres[1,:] = y0
y = copy(y0)
tout = [0.0]
for k in 2:length(t)
flag = Sundials.CVode(mem, t[k], y, tout, Sundials.CV_NORMAL)
yres[k,:] = y
end
Sundials.CVodeFree([mem])
return yres
end
function f(t, y, ydot)
ydot[1] = 0.1*(-72-y[1])+0.1*1.4*exp((y[1]+48)/1.4)+10
ydot[3] = 0.
ydot[2] = 0.
end
immutable J; end
call(::Type{J},t, y, ydot) = f(t, y, ydot)
t = [0.1:0.0001:1]
res = cvode(J, [-60.0, 0.0, 0.0], t);