Mismatch on Figure 7.38

[universmile]

In Figure 7.38, starting from around year 2050, the plot of the variable pop is under the plot of the variable iopc (cf. right subfigure attached below). Instead, following the book figure, pop should be over iopc (cf. left subfigure attached below).

This issue persists after considering different solvers.

In the following, we provide a minimum reproducing example.

using WorldDynamics
using DifferentialEquations

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


system = world3_historicalrun()
sol = WorldDynamics.solve(system, (1900, 2100))


@named pop = World3.Pop15.population()
@named br = World3.Pop15.birth_rate()
@named dr = World3.Pop15.death_rate()
@named is = World3.Capital.industrial_subsector()
@named ss = World3.Capital.service_subsector()
@named ld = World3.Agriculture.land_development()
@named nr = World3.NonRenewable.non_renewable()
@named pp = World3.Pollution.persistent_pollution()
@named se = World3.SupplementaryEquations.supplementary_equations()
@named ai = World3.Agriculture.agricultural_inputs()
@named lfd = World3.Agriculture.land_fertility_degradation()

@variables t

fig_7_38_variables = [
    (nr.nrfr,  0, 1,    "nrfr"),
    (is.iopc,  0, 1000, "iopc"),
    (ld.fpc,   0, 1000, "fpc"),
    (pop.pop,  0, 16e9, "pop"),
    (pp.ppolx, 0, 32,   "ppolx"),
    (br.cbr,   0, 50,   "cbr"),
    (dr.cdr,   0, 50,   "cdr"),
]

cap_tables_7_38 = World3.Capital.gettables()
cap_tables_7_38[:isopc2] = (60, 450, 960, 1500, 1830, 2175, 2475, 2700, 3000)

agr_tables_7_38 = World3.Agriculture.gettables()
agr_tables_7_38[:ifpc2] = (345, 720, 1035, 1275, 1455, 1605, 1725, 1815, 1875)

pop_parameters_7_38 = World3.Pop15.getparameters()
pop_parameters_7_38[:zpgt] = 1975

cap_parameters_7_38 = World3.Capital.getparameters()
cap_parameters_7_38[:iet] = 1990
cap_parameters_7_38[:iopcd] = 320
cap_parameters_7_38[:alic2] = 21
cap_parameters_7_38[:alsc2] = 30

nr_parameters_7_38 = World3.NonRenewable.getparameters()
nr_parameters_7_38[:nruf2] = 0.125

pol_parameters_7_38 = World3.Pollution.getparameters()
pol_parameters_7_38[:ppgf21] = 0.25

system = world3_historicalrun(pop_params=pop_parameters_7_38, capital_params=cap_parameters_7_38, nonrenewable_params=nr_parameters_7_38, pollution_params=pol_parameters_7_38, capital_tables=cap_tables_7_38, agriculture_tables=agr_tables_7_38)
sol_7_38 = WorldDynamics.solve(system, (1900, 2100), solver = AutoVern9(Rodas5()))

plotvariables(sol_7_38, (t, 1900, 2100), fig_7_38_variables, name="Fig. 7.38", showlegend=true, showaxis=true, colored=true)

For sanity check, we also provide herafter the DYNAMO code snippet from the book chapter.

[natema]

We should test what we get if we solve the system like them, i.e. with a simple Euler method with fixed step size of 5 years.

[universmile]

The solution returned by the Euler solver is noisy.

The changes have been made under solvesystems.jl in the following manner (corresponding to the first plot below) :

function solve(system::ODESystem, timespan; solver=Euler())
    sys = structural_simplify(system)

    prob = ODEProblem(sys, [], timespan)
    sol = ModelingToolkit.solve(prob, solver, dt = 5)

    return sol
end

Fixed step size of 5 :

Fixed step size of 50 :

[natema]

Yes, @piluc and I also verified this independently.
Now @piluc and @paulobruno are looking into this.

[natema]

There must be an error in the book, since our Figure agrees with the one of pyworld3.

[paulobruno]

To illustrate what @natema mentioned, here is the comparison of fig 7.38.

World3	WorldDynamics.jl	PyWorld3

[paulobruno]

Since the general agreement is that our code is now returning the correct plot, I'm closing this issue.

[natema]

For future reference, the same mismatch is present in Figure 7.39 as well.