[WIP/discuss] More strict check that equation is conditionally linear

[mstimberg]

@romainbrette noticed that while we claim that the exponential Euler algorithm requires that equations are conditionally linear (linear in the integrated variable, potentially non-linear in all other variables), our implementation actually allows for equations that are non-linear in the integrated variable itself. This is because we simply separate the equation into a linear and a non-linear part, and the non-linear part might still depend on the integrated variable. So currently, you can integrate the following with exponential Euler:

eqs = 'dv/dt = (-v + exp(-v))/tau : 1

This PR fixes this by checking that the non-linear part does not include the variable.

However, after implementing this I started to have some doubts (also because applying the exponential Euler algorithm to these kind of equations gives somewhat reasonable results). There are certainly better sources, but from Wikipedia, the algorithm requires the equation to decompose as:

The exact solution is then:

And you can use a first-order approximation by considering the non-linear part to be constant over the time step, yielding

which is I think exactly what we are doing. Note that this contains an arbitrary non-linear term
which is a function of y, so it seems to be fine to have a non-linear dependence on the integrated variable. Am I missing something here?

[romainbrette]

This seems to be some kind of generalization of exp Euler. My feeling however is that A should include the linearized term, i.e. N'(y)=0 (otherwise you could choose any arbitrary A).

[mstimberg]

Hmm, I'm still not sure. All descriptions of exponential Euler that I can find seem to state the same thing, e.g. here from some lecture notes:

(the phi function is (e^x - 1)/x). We are restricting things currently by integrating all variables separately, but this does not seem to be relevant here. I agree that you could decompose things incorrectly and include the linear term in the non-linear part, but we try our best to not do it I guess?

[romainbrette]

Well but if we use exprel for example, then exprel'(0) = 1/2, not 0. (In my case I used 1/exprel.)
So we do include part of the linear part in the nonlinear term.
Anyway, the point is to deal with the stiffness so it's probably fine in practice.
Since exp euler is not worse than Euler even if improperly applied, then perhaps we should leave it as is.

[mstimberg]

I still haven't made my mind up about this, @thesamovar what do you think? But I agree that introducing this change now (just before a new release) would potentially break code for some users and we should only do this if we are sure it is the right thing to do.

As a side note, we should probably generalize the exponential Euler updater to work on the whole system of equations and use the matrix exponential. For linear equations it should give the same results as the exact algorithm, but it currently solves all equations independently (inherited from Brian 1's approach), linear dependencies on other variables are therefore not correctly taken into account.

[romainbrette]

Probably we shouldn't introduce the change.
For the side note, that would be very useful!

[thesamovar]

This is an old one, but if you both think it's fine I'm sure it's OK as is. Shall we close this PR?