Get back the original discretized system when using Euler with 5y steps

[natema]

As noted here in #79, using the (supposedly) same method than that in World3 messes up many trajectories:

Here's an example:

using WorldDynamics, .World3
using ModelingToolkit, DifferentialEquations 
##
@variables t

nr_parameters_7_10 = World3.NonRenewable.getparameters()
nr_parameters_7_10[:nri] = 2e12
##
system = World3.historicalrun(nonrenewable_params=nr_parameters_7_10)
##
sys = structural_simplify(system)
timespan = (1900, 2100)
prob = ODEProblem(sys, [], timespan)
solution = ModelingToolkit.solve(prob, Euler(), dt=5.0)
##
plotvariables(solution, (t, 1900, 2100), World3.variables_7(); title="Fig. 7.10")

There should be some choices that we made when we derived the differential equations associated to their discrete system, or some auxiliary functions, such that considering the Euler method now is not equivalent.
It would be good to understand where the problem lies.

[paulobruno]

Do you know if at some point in the past Euler method was working properly?
If so, we could try to compare the current version with the old one, but I've never checked this before.

[natema]

No, since I tried the first time it was already oscillating.

[natema]

In an offline discussion today we were mentioning that the Vensim implementation does consider a continuous version, which (AFAIU) is solved using the Euler discretization, so it will be useful to compare our sources.

[natema]

Do you know if at some point in the past Euler method was working properly?

I started investigating the issue by comparing the resulting figure of the tutorial on implementing a new model
when the system is solved with WorldDynamics.solve(nrs_run(), (0, 200), solver=Euler(), dt=5.0).
There everything works well: cfr

with

[natema]

I did a test substituting the resource equation with the interpolate function, providing as table the previous solution: resource ~ interpolate(t, resourcesol, (0., 200.)) with

resourcesol = (1000.0, 975.0, 947.56875, 917.4773545312501, 884.4769486340473, 848.2972918828095, 808.645338703934, 765.2038048417348, 717.7992755943341, 666.1623725143835, 610.3695014863697, 550.444466741117, 487.2970675974564, 421.60214568012793, 355.0076487677415, 288.93030633877447, 225.21948626168356, 168.31131155056394, 119.45699755129705, 80.46824236236091, 53.86274373248295, 37.816035376247534, 28.339094094673595, 22.62140581011455, 19.044241915838743, 16.72268172363289, 15.167571867388984, 14.098944970685153, 13.349842609103222, 12.816644846215786, 12.432699599252006, 12.153799966678665, 11.949870097525709, 11.800021154477228, 11.689504353127594, 11.60777076868378, 11.54719938163122, 11.502241837656044, 11.468834929013973, 11.44398975260486, 11.425500240005112)

It works fine, so the problem is elsewhere.

[natema]

Interesting news: Euler with dt=1.5 works decently.
Instability deteriorates quickly when raising the step above that value.

[natema]

In fact, in Vensim they set dt to less than 0.008, so it seems that there's no issue w.r.t. other implementations.
I'm then closing this and suggest using Euler with dt=0.008.

[paulobruno]

@natema Nice discovery!
Would you have any suggestion to include this information in the docs (or even in the code), so that we don't forget about it in the future?

[natema]

We'll have to update the solver that gets used when we call the methods to compute the runs for the different models. The docstrings of such methods look like a good place.

[paulobruno]

It sounds good to me as well.