Интегриране на векторно поле (масив 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()

Имайте предвид, че направих интеграционната мрежа по-плътна за допълнителна прецизност, но в този случай не се промени много.

Резултат:

Актуализация

След някои бележки в коментарите преразгледах моя оригинален подход, базиран на griddata. Причината за това е, че докато interp2d изчислява интерполант за цялата решетка с данни, griddata изчислява интерполиращата стойност само в дадените му точки, така че в случай на няколко точки последният трябва да бъде много по-бърз.

Поправих грешките в предишния си опит griddata и измислих

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]]

който е съвместим с odeint. Той изчислява интерполираните стойности за всяка p точка, дадена му от odeint. Това решение не консумира излишна памет, но отнема много много повече време, за да работи с горните параметри. Това вероятно се дължи на много оценки на dfun в odeint, много повече от това, което би било очевидно от 100-те времеви точки, дадени му като вход.

Въпреки това, получената рационализация е много по-плавна от тази, получена с interp2d, въпреки че и двата метода използват метода за интерполация по подразбиране linear:

В случай, че някой има израз на полето, използвах кратка версия на отговора на Andras без маската и векторите:

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()

Надявам се да е полезно за някой със същите нужди.