Fluctuations of the Crude Birth Rate variable

[universmile]

In few figures of Chapter 2 (Population), the plot corresponding to the cbr variable fluctuates. This instability might be related to Issue #59 (e.g., iphst parameter).

In the following, we provide a minimal reproducible example.

Source code :

using WorldDynamics

include("../scenarios/pop1_historicalrun.jl")

@named pop = WorldDynamics.World3.Pop1.population()
@named br = WorldDynamics.World3.Pop1.birth_rate()
@named dr = WorldDynamics.World3.Pop1.death_rate()
@named so = WorldDynamics.World3.Pop1.service_output()
@named io = WorldDynamics.World3.Pop1.industrial_output()
@named f = WorldDynamics.World3.Pop1.food()

@variables t

fig_2_88_variables = [
    (so.sopc, 0,   1000,   "sopc"),
    (io.iopc, 0,   1000,   "iopc"),
    (f.fpc,   0,   1000,   "fpc"),
    (pop.pop, 0,   1.6e10, "pop"),
    (br.cbr,  0,   50,     "cbr"),
    (dr.cdr,  0,   50,     "cdr"),
    (dr.le,   0,   80,     "le"),
    (dr.fpu,  0,   1,      "fpu"),
    (br.fce,  0.5, 1,      "fce"),
]

parameters_2_88 = WorldDynamics.World3.Pop1.getparameters()
parameters_2_88[:iphst] = 4000
parameters_2_88[:lt2] = 1900
parameters_2_88[:cio] = 1000
parameters_2_88[:cso] = 1500
parameters_2_88[:cfood] = 2500

system = pop1_historicalrun(params=parameters_2_88)
sol_2_88 = WorldDynamics.solve(system, (1900, 2100))

plotvariables(sol_2_88, (t, 1900, 2100), fig_2_88_variables, name="Fig. 2.88a", showlegend=true, showaxis=true, colored=true)

Resulting plot :

[universmile]

The issue is related to the method used by the solver. We should find the adequate method.

In the following, we show the plots returned by different methods, namely Tsit5 (default), BS5 (for non-stiff ODE), and Rosenbrock23, TRBDF2, Rodas5 (for stiff ODE).

Cf. https://docs.juliahub.com/DifferentialEquations/UQdwS/6.17.1/solvers/ode_solve/#ode_solve

[natema]

I had a look at the Unknown Stiffness section in the docs, and calling

sol_2_88 = WorldDynamics.solve(system, (1900, 2100), solver=AutoVern9(Rodas5()))

gives

[paulobruno]

We set the default solver to AutoVern9(Rodas5()), but I'm still getting the fluctuations.

[natema]

We set the default solver to AutoVern9(Rodas5()), but I'm still getting the fluctuations.

This persists on julia 1.8.1?

[paulobruno]

We set the default solver to AutoVern9(Rodas5()), but I'm still getting the fluctuations.

This persists on julia 1.8.1?

Yes, it does.

[universmile]

I do not get the fluctuations when I move (from Julia 1.7.2) to Julia 1.8.1.

Side information : I am using DifferentialEquations 7.3.0 and ModelingToolkit 8.16.0.

[paulobruno]

This is very strange, I'm also using those versions.

(WorldDynamics) pkg> status
Project WorldDynamics v0.1.0
Status `/home/pdesousa/Documents/git/WorldDynamics/Project.toml`
  [35d6a980] ColorSchemes v3.19.0
  [3da002f7] ColorTypes v0.11.4
⌅ [0c46a032] DifferentialEquations v7.3.0
  [bd48cda9] GraphRecipes v0.5.12
  [615f187c] IfElse v0.1.1
⌅ [961ee093] ModelingToolkit v8.16.0
  [f0f68f2c] PlotlyJS v0.18.9

[natema]

@paulobruno on Julia 1.8.5 AutoVern9(Rodas5()) works fine (I instantiated a fresh environment).
If you still get the issue, we can compare our Manifest files.
Otherwise, this should be closed as soon as we solve #114 .

[paulobruno]

I confirm there are no fluctuations anymore.

[natema]

Since #114 has been closed and the oscillation disappeared, I close this issue.