使用 scipy.integrate 积分向量场（一个 numpy 数组）

Question

我有兴趣使用 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)

编辑：

我需要return周围网格点向量场的插值：

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

Answer 1

我打算建议 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兼容。它计算 odeint 给它的每个 p 点的内插值。此解决方案不会消耗过多的内存，但是使用上述参数需要 多得多 到运行。这可能是由于在 odeint 中对 dfun 进行了大量评估，远远超过从作为输入给它的 100 个时间点可以明显看出的结果。

但是，生成的流线比使用 interp2d 获得的流线平滑得多，即使这两种方法都使用默认的 linear 插值方法：

Answer 2

如果有人对该字段有表达，我使用了没有掩码和向量的简明版本的 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()

希望对有相同需求的人有用

使用 scipy.integrate 积分向量场（一个 numpy 数组）

Integrating a vector field (a numpy array) using scipy.integrate

python

interpolation

scipy

numerical-integration

pde

更新