Есть ли способ легко интегрировать систему дифференциальных уравнений по полной сетке точек? - PullRequest
Новая
       23

Есть ли способ легко интегрировать систему дифференциальных уравнений по полной сетке точек?

0 голосов
/ 18 мая 2019

Проблема в том, что я хотел бы иметь возможность интегрировать дифференциальные уравнения, начиная сразу для каждой точки сетки, вместо того, чтобы обходить интегратор scipy для каждой координаты. (Я уверен, что есть простой способ)

В качестве фона для кода я пытаюсь решить траектории потока Куэтта, чередуя направление скорости каждый определенный период, это хорошо известная динамическая система, которая производит хаос. Я не думаю, что остальная часть кода действительно имеет значение как часть интеграции с scipy и моего использования функции meshgrid numpy. 

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation, writers
from scipy.integrate import solve_ivp

start_T = 100
L = 1
V = 1
total_run_time = 10*3
grid_points = 10

T_list = np.arange(start_T, 1, -1)
x = np.linspace(0, L, grid_points)
y = np.linspace(0, L, grid_points)
X, Y = np.meshgrid(x, y)
condition = True
totals = np.zeros((start_T, total_run_time, 2))
alphas = np.zeros(start_T)
i = 0

for T in T_list:
    alphas[i] = L / (V * T)
    solution = np.array([X, Y])
    for steps in range(int(total_run_time/T)):
        t = steps*T
        if condition:
            def eq(t, x):
                return V * np.sin(2 * np.pi * x[1] / L), 0.0
            condition = False
        else:
            def eq(t, x):
                return 0.0, V * np.sin(2 * np.pi * x[1] / L)
            condition = True
        time_steps = np.arange(t, t + T)
        xt = solve_ivp(eq, time_steps, solution)
        solution = np.array([xt.y[0], xt.y[1]])
        totals[i][t: t + T][0] = solution[0]
        totals[i][t: t + T][1] = solution[1]
    i += 1

np.save('alphas.npy', alphas)
np.save('totals.npy', totals)

Ошибка: 

ValueError: y0 must be 1-dimensional.

И это происходит от функции 'solve_ivp' scipy, потому что она не принимает формат функции numpy meshgrid. Я знаю, что могу запустить несколько циклов и справиться с этим, но я предполагаю, что должен быть «хороший» способ сделать это, используя numpy и scipy. Я также принимаю совет для остальной части кода.

...