Цель:
У меня есть 2 массива vx
и vy
, представляющие компоненты скорости. Я хочу написать алгоритм оптимизации:
seed
)seed
Проблема/Вопрос:
Сначала я написал алгоритм Эйлера, который очень плохо решал мою проблему. Мне посоветовали рассматривать мою проблему как обыкновенное дифференциальное уравнение (ОДУ), где dx/dt = v_x(t) и dy/dt = v_y(t). Я могу интерполировать свои скорости, но мне сложно решить ОДУ с помощью Scipy. Как я мог это сделать?
Самодельный алгоритм:
У меня есть 2 массива vx
и vy
, представляющие компоненты скорости. Когда у одного есть NaN, у другого тоже есть. У меня есть точка, с которой я начинаю, точка seed
. Я хочу отслеживать, через какие ячейки прошла эта точка, основываясь на компонентах скорости.
Я интерполирую компоненты скорости vx
и vy
, чтобы ввести их в решатель ОДУ.
Пример:
Этот код проверяет алгоритм для массива скоростей 10x11. Я заблокирован в решателе ODE.
import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import RegularGridInterpolator
from scipy.integrate import odeint
# Create coordinates
x = np.linspace(0, 10, 100)
y = np.linspace(11, 20, 90)
Y, X = np.meshgrid(x, y)
# Create velocity fields
vx = -1 - X**2 + Y
vy = 1 + X - Y**2
# Seed point
J = 5
I = 14
# Interpolate the velocity fields
interpvx = RegularGridInterpolator((y,x), vx)
interpvy = RegularGridInterpolator((y,x), vy)
# Solve the ODE to get the point's path, but I don't know what to put for the parameter t
#solx = odeint(interpvx, interpvx((I,J)), np.linspace(0,1,501))
#soly = odeint(interpvy, interpvx((I,J)), np.linspace(0,1,501))
Кто-то, с кем я работаю, пришел с решением, я не беру на себя ответственность за ответ. Но, учитывая, что алгоритма оптимизации для Python не существует, надеюсь, он кому-нибудь здесь пригодится:
import numpy as np
# Create coordinates
x = np.linspace(0, 1000, 100)
y = np.linspace(0, 1000, 90)
Y, X = np.meshgrid(x, y)
# Create velocity fields
vx = -1 - X**2 + Y
vy = 1 + X - Y**2
# Seed point
J = 5
I = 14
# Interpolate the velocities
from scipy.interpolate import RegularGridInterpolator
X, Y = np.meshgrid(x,y)
fx = RegularGridInterpolator((y, x), vx, bounds_error=False, fill_value=None)
fy = RegularGridInterpolator((y, x), vy, bounds_error=False, fill_value=None)
# define the velocity function to be integrated:
def f(t, y):
return np.squeeze([fy(y), fx(y)])
# Solve for a seed point
from scipy.integrate import solve_ivp
sol = solve_ivp(f, [0, 100], [J,I], t_eval=np.arange(0,100,1))
Я пытался использовать его, но не смог получить доступ к индексам ячеек, пересекаемых определенными исходными точками. Он рисует, но я не могу получить доступ к промежуточной информации
matplotlib.pyplot.streamplot
существует и работает с теми же входными данными. Только то, что он не производит только одну линию тока, а заполняет всю область графика.