Solving ODE based on experimental data (DifferentialEquations.jl)


I am solving a simple ODE system with 5 state variables.
It is essentially made of 5 coupled RC sytems with a few (current) excitations.
These excitations are 6000 experimental data measured every 10 minutes.
I am using DifferentialEquations.jl and wanted to test the various solvers.
The experimental data are provided to the model by LinearInterpolation.

To my big surprise, almost all solvers failed.
For example:

  • for Euler I got this warning: “Warning: Instability detected. Aborting”

  • for Tsit5: “Warning: Interrupted. Larger maxiters is needed.”
    I checked that the model was accessed 6000001 times !
    Even stranger, the solver only explored the 80 first time steps of the data.

  • the same for RK4, …

However, TRBDF2 solves the problem very easily and accessed the model only 134187 times, and of course accessed all the data.

Any suggestion?



Seems like all the solvers that are failing are based on explicit time-stepping methods. Is the system of ODEs fairly stiff?

1 Like

These equations should not be stiff.
It is a very simple model of heat conduction in a building.
Below are the equations.
The heat sources (i1, …) are measurements averaged over the time step (10 min).
These sources are a bit noisy, but that should not matter a lot.
The time constants are all rather large, except maybe for one of them.

function building(dT,T,p,t)
    T1, T2, T3, T4, T5 = T
    τ1, τ2, τ3, τ4, τ5, U1, U2, U3, U4, U5, k1, AgS = abs.(p)
    T0    = te(t)
    i2    = i21(t) + AgS*i22(t)
    U1k   = 1/(1/U1 + k1*sh(t))
    dT[1] = (i1(t) + U1k*(T0 - T1)  + U2*(T2 - T1)  ) /(τ1*U1)
    dT[2] = (i2 + U2*(T1 - T2)  + U3*(T3 - T2) + U4*(T4 - T2) + U5*(T5 - T2)) /(τ2*U2)
    dT[3] = (i3(t)  + U3*(T2 - T3)  ) /(τ3*U3)
    dT[4] = (i4(t) + U4*(T2 - T4) ) /(τ4*U4)
    dT[5] = (i5(t) + U5*(T2 - T5)  ) /(τ5*U5)

That would likely lead to a form of stiffness in the definition. You might want to use DataInterpolations to generate a regression spline instead of doing a direct interpolation if that’s the case.

Those models are usually stiff.

All of this is pointing to the equation either being stiff, or, the other issue is that

So if you have a lot of data points, then your f is discontinuous at every data point, in which case the error estimators will have a very difficult time handling your problem. You probably want to use a higher order interpolation to smooth it out.

Just observed that ImplicitEuler worked very well.
I overlooked the fact that there is a small time constant in this problem.

Thanks ChrisRackauckas
Thanks jlchan

Further suggestions welcome …