Интеграция векторного поля (массив numpy) с использованием scipy.integrate

Я собирался предложить matplotlib.pyplot.streamplot, который поддерживает аргумент ключевого слова start_points по состоянию на версия 1.5.0, однако она непрактична и очень неточна.

Ваши примеры кода меня немного сбивают с толку: если у вас есть vx, vy координаты векторного поля, то у вас должно быть две сетки: x и y. Используя их, вы действительно можете использовать scipy.interpolate.griddata для получения гладкого векторного поля для интегрирования, однако, когда я пытался это сделать, казалось, что это занимает слишком много памяти. Вот похожее решение на основе scipy.interpolate.interp2d:

import numpy as np
import matplotlib.pyplot as plt
import scipy.interpolate as interp
import scipy.integrate as integrate

#dummy input from the streamplot demo
y, x = np.mgrid[-3:3:100j, -3:3:100j]
vx = -1 - x**2 + y
vy = 1 + x - y**2

#dfun = lambda x,y: [interp.griddata((x,y),vx,np.array([[x,y]])), interp.griddata((x,y),vy,np.array([[x,y]]))]
dfunx = interp.interp2d(x[:],y[:],vx[:])
dfuny = interp.interp2d(x[:],y[:],vy[:])
dfun = lambda xy,t: [dfunx(xy[0],xy[1])[0], dfuny(xy[0],xy[1])[0]]

p0 = (0.5,0.5)
dt = 0.01
t0 = 0
t1 = 1
t = np.arange(t0,t1+dt,dt)

streamline=integrate.odeint(dfun,p0,t)

#plot it
plt.figure()
plt.plot(streamline[:,0],streamline[:,1])
plt.axis('equal')
mymask = (streamline[:,0].min()*0.9<=x) & (x<=streamline[:,0].max()*1.1) & (streamline[:,1].min()*0.9<=y) & (y<=streamline[:,1].max()*1.1)
plt.quiver(x[mymask],y[mymask],vx[mymask],vy[mymask])
plt.show()

Обратите внимание, что я сделал сетку интегрирования более плотной для большей точности, но в данном случае это мало что изменило.

Результат:

xyarr = np.array(zip(x.flatten(),y.flatten()))
dfun = lambda p,t: [interp.griddata(xyarr,vx.flatten(),np.array([p]))[0], interp.griddata(xyarr,vy.flatten(),np.array([p]))[0]]

Если у кого-то есть выражение поля, я использовал краткую версию ответа Андраса без маски и векторов:

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import ode

vx = lambda x,y: -1 - x**2 + y
vy = lambda x,y: 1 + x - y**2

y0, x0 = 0.5, 0.6

def f(x, y):
    return vy(x,y)/vx(x,y)

r = ode(f).set_integrator('vode', method='adams')
r.set_initial_value(y0, x0)

xf =  1.0
dx = -0.001

x, y = [x0,], [y0,]
while r.successful() and r.t <= xf:
    r.integrate(r.t + dx)
    x.append(r.t + dx)
    y.append(r.y[0])

#plot it
plt.figure()
plt.plot(x, y)
plt.axis('equal')
plt.show()

Я надеюсь, что это полезно для кого-то с такими же потребностями.