python - 为什么我的2D内插器在SciPy中生成带有交换轴的矩阵？

我用矢量输入解一个微分方程
y'=f（t，y），y（t_0）=y_0
式中y0=y（x）
使用显式的euler方法，也就是说
y_（i+1）=y_i+h*f（t_i，y_i）
其中t是时间向量，h是步长，f是微分方程的右侧。
方法的python代码如下所示：

for i in np.arange(0,n-1):
        y[i+1,...] = y[i,...] + dt*myode(t[i],y[i,...])

结果是k，m矩阵y，其中k是t维的大小，m是y的大小。
返回向量y和t。
t、x和y被传递到scipy.interpolate.RectBivariateSpline(t, x, y, kx=1, ky=1)：

g = scipy.interpolate.RectBivariateSpline(t, x, y, kx=1, ky=1)

得到的对象g采用新的向量ti,xi ( g(p,q) )给出y_int，这是在TI和XI定义的点上进行Y插值的。
这是我的问题：
rectBivariateSpline的文档用x和y描述了__call__方法：
__call__(x, y[, mth])在坐标数组定义的网格点处计算样条曲线
绘图表面的matplotlib文档使用类似的符号：
Axes3D.plot_surface(X, Y, Z, *args, **kwargs)
重要的区别是X和Y是由numpy.meshgrid()生成的二维数组。
当我计算简单的例子时，两者的输入顺序都是一样的，结果正是我所期望的然而，在我的显式Euler示例中，初始顺序是ti,xi，但插值输出的曲面图只有在我反转输入顺序时才有意义，如下所示：

ax2.plot_surface(xi, ti, u, cmap=cm.coolwarm)

虽然我很高兴它能工作，但我不满意，因为我无法解释为什么，也无法解释为什么（除了阵列几何）必须交换输入。理想情况下，我希望重新构造代码，使输入顺序一致。
下面是一个工作代码示例来说明我的意思：

# Heat equation example with explicit Euler method
import numpy as np
import matplotlib.pyplot as mplot
import matplotlib.cm as cm
import scipy.sparse as sp
import scipy.interpolate as interp
from mpl_toolkits.mplot3d import Axes3D
import pdb

# explicit Euler method
def eev(myode,tspan,y0,dt):
    # Preprocessing
    # Time steps
    tspan[1] = tspan[1] + dt
    t = np.arange(tspan[0],tspan[1],dt,dtype=float)
    n = t.size
    m = y0.shape[0]
    y = np.zeros((n,m),dtype=float)
    y[0,:] = y0

    # explicit Euler recurrence relation
    for i in np.arange(0,n-1):
        y[i+1,...] = y[i,...] + dt*myode(t[i],y[i,...])

    return y,t

# generate matrix A
# u'(t) = A*u(t) + g*u(t)
def a_matrix(n):
    aa = sp.diags([1, -2, 1],[-1,0,1],(n,n))
    return aa

# System of ODEs with finite differences
def f(t,u):
    dydt = np.divide(1,h**2)*A.dot(u)
    return dydt

# homogenous Dirichlet boundary conditions
def rbd(t):
    ul = np.zeros((t,1))
    return ul

# Initial value problem -----------

def main():
    # Metal rod
    # spatial discretization
    # number of inner nodes
    m = 20
    x0 = 0
    xn = 1
    x = np.linspace(x0,xn,m+2)
    # Step size
    global h
    h = x[1]-x[0]

    # Initial values
    u0 = np.sin(np.pi*x)

    # A matrix
    global A
    A = a_matrix(m)

    # Time
    t0 = 0
    tend = 0.2
    # Time step width
    dt = 0.0001
    tspan = [t0,tend]

    # Test r for stability
    r = np.divide(dt,h**2)
    if r <= 0.5:
        u,t = eev(f,tspan,u0[1:-1],dt)
    else:
        print('r = ',r)
        print('r > 0.5. Explicit Euler method will not be stable.')

    # Add boundary values back
    rb = rbd(t.size)
    u = np.hstack((rb,u,rb))

    # Interpolate heat values
    # Create interpolant. Note the parameter order
    fi = interp.RectBivariateSpline(t, x, u, kx=1, ky=1)

    # Create vectors for interpolant
    xi = np.linspace(x[0],x[-1],100)
    ti = np.linspace(t0,tend,100)

    # Compute function values from interpolant
    u_int = fi(ti,xi)

    # Change xi, ti in to 2D arrays
    xi,ti = np.meshgrid(xi,ti)

    # Create figure and axes objects
    fig3 = mplot.figure(1)
    ax3 = fig3.gca(projection='3d')
    print('xi.shape =',xi.shape,'ti.shape =',ti.shape,'u_int.shape =',u_int.shape)

    # Plot surface. Note the parameter order, compare with interpolant!
    ax3.plot_surface(xi, ti, u_int, cmap=cm.coolwarm)
    ax3.set_xlabel('xi')
    ax3.set_ylabel('ti')

main()
mplot.show()

最佳答案

正如我所见，你定义了：

# Change xi, ti in to 2D arrays
    xi,ti = np.meshgrid(xi,ti)

将此更改为：

ti,xi = np.meshgrid(ti,xi)

和

ax3.plot_surface(xi, ti, u_int, cmap=cm.coolwarm)

到

ax3.plot_surface(ti, xi, u_int, cmap=cm.coolwarm)

如果我听得懂的话，它也能很好地工作。