python - 矢量化/广播与使用solve_ivp求解ODE之间的冲突

Question

使用 NumPy 数组和矢量化，我尝试创建一个由 n 个不同个体组成的群体，每个个体具有三个属性:alpha、beta 和 表型(表型被计算为涉及alpha和beta<的微分方程的稳态)。所以，我希望每个人都有自己的表型。

但是，我的代码为每个人产生相同的表型。此外，只有当 solve_ivp 的 y0 数组(这里是 [ 0, 1]) -- 否则，会产生广播错误:

ValueError: operands could not be broadcast together with shapes (2,) (3,)

代码如下:

import numpy as np
from scipy.integrate import solve_ivp

def create_population(n):
    """creates a population of n individuals"""
    pop = np.zeros(n, dtype=[('alpha','<f8'),('beta','<f8'),('phenotype','<f8')])
    pop['alpha'] = np.random.randn(n)
    pop['beta'] = np.random.randn(n) + 5
    def phenotype(n):
        """creates the phenotype"""
        def pheno_ode(t_ode, y):
            """defines the ode for the phenotype"""
            dydt = 0.123 - y + pop['alpha'] * (y ** pop['beta'] / (1 + y ** pop['beta']))
            return dydt
        t_end = 1e06
        sol = solve_ivp(pheno_ode, [0, t_end], [0, 1], method='BDF')
        return sol.y[0][-1] # last entry is assumed to be the steady state
    pop['phenotype'] = phenotype(n)
    return pop

popul = create_population(3)
print(popul)

相反，如果表型是通过“简单”方程根据 alpha 和 beta 计算出来的，则矢量化效果很好:

def phenotype(n):
    """creates the phenotype"""
    phenotype_simple = 2 * pop['alpha'] + pop['beta']
    return phenotype_simple

Answer 1

我发现有两个问题:

首先，您将 ODE 的初始条件设置为 [0, 1] 。设置 solve_ivp 的向量解的大小为 2，无论 n 的值如何。然而，数组pop['alpha']和pop['beta']有长度n ，并在您的脚本中调用 create_population与 n设置为 3。因此 dydt 公式中的数组形状不匹配:y长度为 2，但是 pop['alpha']和pop['beta']长度为 3。这会导致您看到的错误。

您可以使用 np.ones(n) 来解决此问题而不是[0, 1]作为调用 solve_ivp 的初始条件.

第二个问题出现在声明return sol.y[0][-1]中在函数 phenotype(n) 中。 sol.y有形状(n, num_points) ，其中num_points是 solve_ivp 计算得出的点数。所以sol.y[0]只是解决方案的第一个组成部分，并且 sol.y[0][-1]是第一个分量的解的最后一个值。它是一个标量，因此当您执行 pop['phenotype'] = phenotype(n) 时，您为所有表型分配相同的值(第一个组件的稳定状态)。

返回语句应该是 return sol.y[:, -1] 。这将返回解决方案数组的最后列(即所有稳态表型)。