对于给定的微分方程，如何求解库兹涅茨曲线增长模型

Question

上下文：

使用麦迪森项目的实际人均GDP数据集，我使用最小二乘法推导出以下方程：0.012406+0.005132ln(g)-0.002304ln**2

当我试图预测不同经济群体到2050年的人均GDP时，我引用了本文的方法"；Tilman等人10.1073/pnas.11164317108"；至于他们是如何尝试用微分方程求解的：dG/dt=G(-0.6284+0.157lnG-0.0093ln(G)^2

我用同样的方式转换了我的最小二乘法：-0.012406g+g0.005132math.log(g)-g0.001304*math.log)**2

我试图在python中求解几个初始值的ODE，以获得库兹涅茨曲线并估计2050年。我使用了下面的代码，但我无法解决同样的问题。

Python代码：

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.

---

关于这一点的说明会有所帮助。

Answer 1

当odeint调用用户提供的func参数时，它将传入状态向量y的数组。即使系统是单个标量方程，y也将作为长度为1的一维数组传入。在您的情况下，这意味着g将是长度为1的1-d。你在gdp:中有这个

dgdt = [-0.012406*g  + g*0.005132*math.log(g) - g*0.006304*math.log(g)**2]
return dgdt

问题是dgdt的表达式中有多余的括号。g已经是一维数组，因此-0.012406*g也是。事实证明，math.log(g)实际上返回了一个标量，但一个标量加上一个长度为1的1-d NumPy会达到你的预期，所以在这种情况下，它不会造成问题。因此，完整表达式-0.012406*g + g*0.005132*math.log(g) - g*0.006304*math.log(g)**2的结果也是一维NumPy数组(长度为1)。当你把结果用括号括起来时；"深度"；数据结构的。实际上，您创建了一个二维数组(具有平凡的形状(1，1))。例如，在下面的例子中，w是长度为1:的一维阵列

In [23]: w = np.array([1.0])
In [24]: np.shape(w)
Out[24]: (1,)
In [25]: np.shape([w])
Out[25]: (1, 1)

注意，[w]具有形状(1, 1)，即它是2d-阵列。odeint期望用户的函数返回一个一维数组，但它从函数中得到一个二维数组并引发错误。

修复方法很简单：删除那些括号：

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

(注意，我删除了y = g行，它没有任何作用。)