InexactError in ode23s

[lruthotto]

Hey folks,

I'm getting an inexact error in ode23s when solving the below problem:

Does anybody have an idea how to fix it?

Thanks,
Lars

using ODE

function f(t,y)
    fy = zeros(3)
    fy[1] = - 0.04*y[1] + 10^4*y[2]*y[3]
    fy[2] =   0.04*y[1] - 10^4*y[2]*y[3] - 3*10^7*y[2]^2
    fy[3] = 3*10^7*y[2]^2
    return fy
end

function J(t,y)
    Jy = [-0.04 10^4*y[3] 10^4*y[2]; 
        0.04 -6*10^7*y[2]-10^4*y[3] -10^4*y[2];
        0    6*10^7*y[2] 0
        ]
    return Jy
end

y0 = [1;0;0.0;]
@time T,Y = ode23s(f,y0,[0 100],jacobian=J)

[acroy]

Try @time T,Y = ode23s(f,y0, [0. 100.], jacobian=J) (note the floats in the time span array). Currently we use the eltype of tspan to construct the array, which holds the output times. But since the intermediate times are floats, we cant store them in an array of integers ...

[lruthotto]

Thanks for the quick response. It helped.

I must say, however, that it might be a bit inconvenient and unexpected for most users that the code breaks because of that. Would it be a better idea to construct tspan as floats by default?

[acroy]

Yes, this definitely needs to be fixed. I think at some point we already discussed this issue and agreed that tspan should be converted to floats. It's just that nobody got around to do it yet ... (PRs are always greatly appreciated :-)

[tomasaschan]

Unfortunately it's a little more complicated than just unconditionally wrapping the input array in float(), since we might get inputs of a bunch of other cool number types that don't support float conversion but do support everything we require to solve ODE's. The entire discussion is in #26 if you want to read it from the start :)

[acroy]

@tlycken : That is true, but it seems to me that the solution in #26 consists of 3 parts: promotion of elements in tspan to a common floating point type, using this type to set up tout storage and yout storage. The more complicated part is probably the first one (although Steven has already given a sketch of the implementation). But if we skip that one and just use typeof(float(eltype(tspan))) for the other two parts (instead of eltype(tspan)) we would probably fix the most common issues like the one reported above. We should still address the promotion, but this can be done separately if necessary.