Why is my 2D interpolant generating a matrix with swapped axes in SciPy?

I solve a differential equation with vector inputs y' = f(t,y), y(t_0) = y_0 where y0 = y(x) using the explicit Euler method, which says that y_(i+1) = y_i + h*f(t_i, y_i) where t is a time vector, h is the step size, and f is the right-hand side of the differential equation. The python code for the method looks like this: ``` for i in np.arange(0,n-1): y[i+1,...] = y[i,...] + dt*myode(t[i],y[i,...]) ``` The result is a k,m matrix y, where k is the size of the t dimension, and m is the size of y. The vectors y and t are returned. t, x, and y are passed to `scipy.interpolate.RectBivariateSpline(t, x, y, kx=1, ky=1)`: ``` g = scipy.interpolate.RectBivariateSpline(t, x, y, kx=1, ky=1) ``` The resulting object `g` takes new vectors `ti,xi ( g(p,q) )` to give `y_int`, which is y interpolated at the points defined by ti and xi. Here is my problem: The documentation for RectBivariateSpline describes the `__call__` method in terms of x and y: `__call__(x, y[, mth])` Evaluate spline at the grid points defined by the coordinate arrays The matplotlib documentation for plot_surface uses similar notation: `Axes3D.plot_surface(X, Y, Z, *args, **kwargs)` with the important difference that X and Y are 2D arrays which are generated by `numpy.meshgrid()`. When I compute simple examples, the input order is the same in both and the result is exactly what I would expect. In my explicit Euler example, however, the initial order is `ti,xi`, yet the surface plot of the interpolant output only makes sense if I reverse the order of the inputs, like so: ``` ax2.plot_surface(xi, ti, u, cmap=cm.coolwarm) ``` While I am glad that it works, I'm not satisfied because I cannot explain why, nor why (apart from the array geometry) it is necessary to swap the inputs. Ideally, I would like to restructure the code so that the input order is consistent. Here is a working code example to illustrate what I mean: ``` # 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() ```