Система двух многомерных связанных ОДУ

Я пытаюсь решить следующую проблему связанных ОДУ, используя odeint() из scipy. Система выглядит так:

X'_k = среднее (Y_k) + F

Y'_{k,j} = X_k - Y_{k,j}

Это система с 3 переменными X, и для каждой переменной X есть еще 3 переменные Y.

Из того, что я прочитал в документации, и примеры здесь и здесь, я могу передать систему уравнений как список. И это то, что я пробовал в следующем примере:

import numpy as np
from scipy.integrate import odeint

def dZdt(Z, t):

    X = Z[0]
    Y = Z[1]

    F = 4

    d_x = np.zeros(3)
    d_y = np.zeros(3*3).reshape(3,3)

    # Compute the Y values
    for k in range(3):
        for j in range(3):
            d_y[k][j] = X[k] - Y[k][j]

        # X values
        d_x[k] = Y[k].mean() + F

    d = [d_x, d_y]

    return d

# Initial conditions
X0 = np.random.uniform(size=3)
Y0 = np.random.uniform(size = 3*3).reshape(3,3)
Z0 = [X0, Y0]

t = range(20)

Z = odeint(dZdt, Z0, t)

Где k, j = (1,2,3) и Z = [X,Y]

Но я боюсь, что получаю следующую ошибку:

ValueError: could not broadcast input array from shape (3,3) into shape (3)

Моя реальная проблема сложнее, потому что j и k могут быть больше 3 (они идут от 1 до j_max и K_max соответственно), поэтому я не могу записать 12 переменных одну за другой.

Я предполагаю, что где-то в коде переменные Y пытаются заполнить форму X... но не знаю, где.

Есть идеи, что я делаю неправильно?


python scipy odeint
person Xbel    schedule 09.10.2018    source источник

Ответы (1)


arrow_upward
2
arrow_downward

Вы пытаетесь представить неизвестную функцию двумя массивами внутри списка. Это должен быть одномерный массив. Таким образом, вместо 3 X-переменных и 9 Y-переменных это должен быть плоский список из 12 переменных. Как это:

def dZdt(Z, t):

    X = Z[:3]
    Y = Z[3:].reshape(3, 3)
    F = 4

    d_x = np.zeros(3)
    d_y = np.zeros((3, 3))

    # Compute the Y values
    for k in range(3):
        for j in range(3):
            d_y[k, j] = X[k] - Y[k, j]

        # X values
        d_x[k] = Y[k].mean() + F

    d = np.concatenate((d_x.ravel(), d_y.ravel()))

    return d

# Initial conditions
X0 = np.random.uniform(size=3)
Y0 = np.random.uniform(size=(3, 3))
Z0 = np.concatenate((X0.ravel(), Y0.ravel()))

t = range(20)

Z = odeint(dZdt, Z0, t)

Массивы NumPy индексируются как Y[k, j], а не Y[k][j]. И есть широкие возможности векторизации, которые устранили бы циклы в вычислении dZdt. Как это:

def dZdt(Z, t):

    X = Z[:3]
    Y = Z[3:].reshape(3, 3)
    F = 4

    d_y = X[:, None] - Y 
    d_x = Y.mean(axis=1) + F

    d = np.concatenate((d_x.ravel(), d_y.ravel()))

    return d
person Community    schedule 09.10.2018