stiff ode integrator gives wrong result

[rajeev2010]

The stiff ode integrator ode4s_kr gives wrong answers instead of failing. The code given below is an example of this behavior

load("ode.jl")

function func(t, y)
    no_osc = int (size(y,1) / 2 );
    dy = zeros(no_osc * 2);
    dy[1:no_osc]        = y[no_osc+1:2*no_osc];
    dy[no_osc+1:2*no_osc] = -y[1:no_osc] - 1000. * y[no_osc+1:2*no_osc] .*  (y[1:no_osc].*y[1:no_osc]-1.);
    return dy
end

no_obs = 101
tspan = linspace(0.,2*10^3,no_obs)
y_0   = rand(4)

t, y = ode4s_kr(func, tspan, y_0)

[pao: formatting]

[pao]

This is probably down to the estimation of the Jacobian. The implementation is a straight MATLAB port of some code I wrote in a simulation class, and I'm not surprised it's not good enough. See also #615; I actually never figured out how to make the API work for providing a different Jacobian function. Maybe after a few months not working on it I'll be able to sort it out.

Better algorithms for estimating the Jacobian are also welcome.

[pao]

I've just gone back to take a look at ode.jl and I see that I apparently did add a signature for a user-provided Jacobian:

ode4s_kr(F, G, tspan, x0) = oderosenbrock(F, G, tspan, x0, kr4_coefficients...)

Where G is a function of the form (t, x) -> J (i.e., given the time and the state vector, return the Jacobian.) If your system has an analytic Jacobian matrix please try this form and see if there's still a problem. If so, there's still a bug in the solver (I don't think there is after #615), but if not, then I can change this to "implement a non-toy Jacobian estimator by default."

[rajeev2010]

It doesn't seem to work even when the jacobian is provided. Below is the jacobian for the above example -

function jac(t, y)
    no_osc = int (size(y,1) / 2 )
    jy = zeros(no_osc*2, no_osc*2)

    for ii = 1 : no_osc
        jy[ii,no_osc+ii] = 1
        jy[no_osc+ii,ii] = -1 - 2000. * y[ii] * y[no_osc+ii] 
        jy[no_osc+ii,no_osc+ii] = -1000. * ( y[ii] * y[ii] - 1.0 )
    end

    return jy
end

[pao]

Just noticed your timesteps are pretty large. The stiff solver we have is not adaptive, unlike MATLAB's ode15s (which is why it can't "fail", there's no metric to test like ode45 has using it's embedded RKF step.) The ode23 and ode45 solvers are adaptive. You note that GSL's rkf45 works, but it looks like it uses the Fehlberg coefficients. That's available in Julia as ode45_fb (The Dormand-Prince coefficients are used by default); that may still not work but at least it will be apples-to-apples. Alternatively, crank down the timesteps with ode4s/ode4s_kr.

[rajeev2010]

ode45_fb is also not working. i get the following error -

Step size grew too small. t=0.0, h=7.52026195675261e-29, x=[0.302402, 0.936122, 0.264298, 0.60962]

I was getting the same error in ode45 also. I can not crank down the timesteps for ode4s etc because the period of the system is very large and it will make the output a very large array which would be a problem as my system of interest is already a very large one.

[pao]

For the sake of understanding whether the stiff method is implemented correctly, could you run a smaller timestep over a shorter period of time to see if it draws the solution closer?

The ode45 code is adapted from Octave; I wonder if the original implementation there has the same problem adapting the timestep?

It's sounding likely that none of the current algorithms are suitable for your problem, which makes this similar to #75.

As I said before, more or improved solvers are certainly a welcome contribution.

[pao]

housekeeping: see SciML/ODE.jl#1 for continuation.