库兹涅茨曲线增长模型微分方程Python求解及报错修复

基于Maddison项目人均实际GDP数据集，通过最小二乘法推导得到如下方程：
0.012406 + *0.005132*ln(g) - *0.006304*ln(g)²

因需要预测不同经济群组到2050年的人均GDP，参考Tilman等人（文献DOI：10.1073/pnas.1116437108）中将同类关系作为微分方程求解的方法论，参考方程形式为：
dG/dt = G(-0.6284 + 0.157lnG - 0.0093ln(G)²)

按照相同方式将最小二乘结果转化为常微分方程（ODE）形式：
-0.012406*g + g*0.005132*math.log(g) - g*0.006304*math.log(g)**2

在Python中针对多组初始值求解该ODE，用于绘制库兹涅茨曲线并得到2050年的GDP预测值时，代码运行失败。

import numpy as np
import matplotlib.pyplot as plt
import scipy as sp
from scipy.integrate import odeint
from scipy.integrate import solve_ivp 
import math

def solveit(y0):
    def gdp(g, t):
        y = g
        dgdt = [-0.012406*g  + g*0.005132*math.log(g) - g*0.006304*math.log(g)**2]
        return dgdt

#initial conditions
    #y0 = [785.60] 
    t = np.linspace(0, 60000, 1000)
#call integrator
    sol = odeint(gdp, y0, t)
    m = sol[:]
    plt.plot(t,m)
    plt.show()


ys= [[785.60],[1860],[7800]]

fig = plt.figure()
for y_ in ys:
    solveit(y_)

plt.legend(loc='best')
plt.grid()
plt.show()

RuntimeError: The array return by func must be one-dimensional, but got ndim=2.

核心报错原因 #

• 初始值集合ys为嵌套列表结构，每个传入求解器的初始值都是单元素列表（如[785.60]），属于二维输入
• 微分方程定义中，导数计算结果被额外包裹在列表中返回，和二维初始值叠加后返回值维度为2，不符合odeint要求的导数结果、初始值必须为一维结构的规则

其他逻辑问题 #

• 时间范围设置为0到60000，跨度远大于2050年的预测需求，极易引发数值积分溢出、结果发散
• 每次调用solveit都会单独触发plt.show()弹出独立窗口，后续统一设置图例、网格的代码无法作用到这些窗口，无法得到合并曲线的效果图

修复后可运行代码 #

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
import math

def solveit(y0, ax):
    def gdp(g, t):
        # 直接返回一维标量计算结果，不额外嵌套列表
        dgdt = -0.012406*g  + g*0.005132*math.log(g) - g*0.006304*math.log(g)**2
        return dgdt

    # 时间范围调整为基期到2050年的实际跨度，示例设为30年，可根据自身基期年份修改
    t = np.linspace(0, 30, 100)
    # 用flatten强制拉平结果为一维，避免绘图维度异常
    sol = odeint(gdp, y0, t).flatten()
    ax.plot(t, sol, label=f'初始人均GDP: {y0}')
    return sol[-1] # 返回时间序列末端值，即2050年预测结果


# 初始值直接传入标量，不做列表嵌套
ys = [785.60, 1860, 7800]
fig, ax = plt.subplots()

for y_ in ys:
    gdp_2050 = solveit(y_, ax)
    print(f"初始值{y_}对应的2050年人均GDP预测值：{gdp_2050:.2f}")

ax.legend(loc='best')
ax.grid()
ax.set_xlabel('距基期年数')
ax.set_ylabel('人均实际GDP')
plt.show()

内容的提问来源于stack exchange，提问作者ShanksPuff