CODEWITHSUNDEEP
pythonnumpymatplotlibscipyipython
Stephen Bosch
Stephen Bosch

Reputation: 273

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:

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

Upvotes: 4

Views: 399

Answers (1)

George
George

Reputation: 5691

As I can see you define :

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

Change this to :

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

and 

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

to

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

and it works fine (if I understood well ).

Upvotes: 3

Related Questions