使用scipy.integrate集成矢量字段(numpy数组)

我有兴趣使用scipy.integrate库为给定的初始点集成矢量场(即查找流线).由于向量字段是在计算网格上定义的numpy.ndarray对象,因此必须对网格点之间的值进行插值.有任何集成商可以处理吗？也就是说,例如,如果我要尝试以下操作

import numpy as np
import scipy.integrate as sc
vx = np.random.randn(10,10)
vy = np.random.randn(10,10)
def f(x,t):
    return [vx[x[0],x[1]], vy[x[0],x[1]]] # which obviously does not work if x[i] is a float
p0 = (0.5,0.5)
dt = 0.1
t0 = 0
t1 = 1
t = np.arange(t0,t1+dt,dt)
sc.odeint(f,p0,t)

编辑：

我需要返回周围网格点的矢量场的插值：

def f(x,t):
    im1 = int(np.floor(x[0]))
    ip1 = int(np.ceil(x[1]))
    jm1 = int(np.floor(x[0]))
    jp1 = int(np.ceil(x[1]))
    if (im1 == ip1) and (jm1 == jp1):
        return [vx[x[0],x[1]], vy[x[0],x[1]]]
    else:
        points = (im1,jm1),(ip1,jm1),(im1,jp1),(ip1,jp1)
        values_x = vx[im1,jm1],vx[ip1,jm1],vx[im1,jp1],vx[ip1,jp1]
        values_y = vy[im1,jm1],vy[ip1,jm1],vy[im1,jp1],vy[ip1,jp1]
        return interpolated_values(points,values_x,values_y) # how ?

最后的return语句只是一些伪代码.但这基本上就是我想要的.

编辑：

scipy.interpolate.griddata函数似乎是可行的方法.是否可以将其合并到自身的功能中？符合以下内容：

    def f(x,t):
        return [scipy.interpolate.griddata(x,vx),scipy.interpolate.griddata(x,vy)]

解决方法:

我本来建议matplotlib.pyplot.streamplot,从1.5.0版开始支持关键字参数start_points,但是它不切实际,而且也很不准确.

您的代码示例让我有些困惑：如果您有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()

请注意,我使集成网格更加密集以提高精度,但是在这种情况下,它并没有太大变化.

结果：

更新资料

在发表评论后,我重新审视了我最初的基于网格数据的方法.这样做的原因是,虽然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兼容.它为odeint给定的每个p点计算插值.此解决方案不会消耗过多的内存,但是使用上述参数运行会花费更长的时间.这可能是由于对odeint中的dfun进行了许多评估,远远超过了从作为输入提供的100个时间点中得出的结论.

但是,即使两种方法都使用默认的线性插值方法,最终的流线也比使用interp2d获得的流线要平滑得多：

标签：scipy,numerical-integration,pde,python,interpolation
来源： https://codeday.me/bug/20191119/2035856.html