Почему мой 2D-интерполант генерирует матрицу с измененными осями в SciPy?

Question

Я решаю дифференциальное уравнение с векторными входами

y' = f(t,y), y(t_0) = y_0

где у0 = у (х)

используя явный метод Эйлера, который говорит, что

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

Результатом является матрица y k,m, где 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.

Вот моя проблема:

• Документация для RectBivariateSpline описывает __call__ Метод с точки зрения х и у:__call__(x, y[, mth]) Оценить сплайн в точках сетки, определенных массивами координат

• Документация matplotlib для plot_surface использует похожую запись:Axes3D.plot_surface(X, Y, Z, *args, **kwargs)с важным отличием, что X и Y являются двумерными массивами, которые генерируются numpy.meshgrid(),

Когда я вычисляю простые примеры, порядок ввода одинаков для обоих, и результат в точности соответствует ожиданиям. Однако в моем явном примере Эйлера начальный порядок 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()

Answer 1

Как я вижу, вы определяете:

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

и работает нормально (если я хорошо понял).