Конфликт между векторизацией / трансляцией и решением ODE с solve_ivp

Question

Используя массив NumPy и векторизацию, я пытаюсь создать популяцию из n разных индивидуумов, каждый из которых имеет три свойства: alpha, beta и phenotype (phenotypeрассчитывается как стационарное состояние дифференциального уравнения, включающего alpha и beta).Итак, я хочу, чтобы у каждого человека был свой фенотип.

Однако мой код выдает один и тот же фенотип для каждого человека.Более того, это нежелательное поведение возникает только в том случае, если в массиве solve_ivp y0 (который здесь [0, 1]) находится ровно n записей, в противном случае возникает ошибка широковещания:

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 должен быть return sol.y[:, -1].Это возвращает последний столбец массива решения (т.е. все фенотипы стационарного состояния).