Dear Dan,
Changing the function to (say) g doesn't help. I want cvode to use J as a
function (as I'm overloading call); this is really just doing what
FastAnonymous does under the hood. I think it may be because Sundials
passes the derivatives as an argument to the function, and modifies them in
place.
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 g(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) = g(t, y, ydot)
t = [0.1:0.0001:1]
res = Sundials.cvode(g, [-60.0, 0.0, 0.0, t); # this works, passing a
Function type
res = cvode(J, [-60.0, 0.0, 0.0], t); # this gives ReadOnlyMemoryError
On Saturday, December 19, 2015 at 12:52:23 AM UTC, Dan wrote:
>
> The parameter for the `cvode` function is `f` and so is the function you
> want to use. These get confused, and it tries to use the "function" `J`
> instead. Changing the parameter name to something other than `{f}` should
> work
>
> On Wednesday, December 16, 2015 at 4:26:41 PM UTC+2, Simon Frost wrote:
>>
>> 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);
>>
>
