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

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

0 голосов
/ 09 октября 2018

Я пытаюсь решить следующую проблему связанных ODE, используя 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 ... но не знаю, где.

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

1 Ответ

0 голосов
/ 09 октября 2018

Вы пытаетесь представить неизвестную функцию двумя массивами внутри списка.Это должен быть одномерный массив.Таким образом, вместо 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
...