Little bump fluctuation at the initialization

[universmile]

For some figures of Chapter 2 and Chapter 7, there is an incorrect bump during the initialization. This bump concerns most of the time the variable cbr.

While some major fluctuations of the cbr variable have been fixed by changing the solver (cf. #60), such initial bump resist and seems unrelated to the choice of the solver.

Examples :

Example affecting other variables in addition to the cbr variable :

[natema]

This initial bump can be a simple consequence of the fact that we solve the system in continuous time: you can check that the behavior is the same after approx. 5 years after the start.
The bump should then disappear when we solve the system in the same way, namely with Euler() and stepsize dt=5.0.
I will test it now.

A related point is that the current solve interface doesn't allow to pass the stepsize dt keyword argument when choosing Euler() as integration method.
We might open an issue regarding enabling the possibility to pass to our solve any argument that is supported by ModelingToolkit.

[natema]

I'm puzzled... using the (supposedly) same method than that in World3 messes up many trajectories:

Here's the source:

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")

Maybe I'm forgetting about some choices that we made when we derived the differential equations associated to their discrete system, and considering the Euler method now is not equivalent.
This should be further investigated...

[paulobruno]

With the changes in #94, @piluc decided to remove the artificial initialization of frsn that used a clip function. frsn directly affects the values of several other variables, including "crude birth rate" (cbr).
As a result, there is no bump anymore, but the initialization starts lower from the beginning, as you can see in the figure below.

As @natema pointed out, this is not necessarily an error in our side, but a discrepancy between the discretization used in the book and our solution with ModelingToolkit.
If everyone is on the same page, I think we can close this issue.