Matplotlib cmap not working as expected

Question

I have the following code for Python 3 which generates a wave function over time and plot the result in 3D. Notice that the schroedinger1D(...) function returns two numpy arrays, each of shape (36,1000).

import numpy as np
from scipy.integrate import fixed_quad
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm

# Initial conditions, outside for plotting.
f_re = lambda x: np.exp(-(x-xc)**2.0/s)*np.cos(2.0*np.pi*(x-xc)/wl)
f_im = lambda x: np.exp(-(x-xc)**2.0/s)*np.sin(2.0*np.pi*(x-xc)/wl)

def schroedinger1D(xl, xr, yb, yt, M, N, xc, wl, s):
    """
        Schrödinger Equation Simulation (no potential)
    """
    f = lambda x: f_re(x)**2 + f_im(x)**2
    area = fixed_quad(f, xl, xr, n=5)[0]
    f_real = lambda x: f_re(x)/area
    f_imag = lambda x: f_im(x)/area
    # Boundary conditions for all t
    l = lambda t: 0*t
    r = lambda t: 0*t
    # "Diffusion coefficient"
    D = 1
    # Step sizes and sigma constant
    h, k = (xr-xl)/M, (yt-yb)/N
    m, n = M-1, N
    sigma = D*k/(h**2)
    print("Sigma=%f" % sigma)
    print("k=%f" % k)
    print("h=%f" % h)
    # Finite differences matrix
    A_real = np.diag(2*sigma*np.ones(m)) + np.diag(-sigma*np.ones(m-1),1) + np.diag(-sigma*np.ones(m-1),-1)
    A_imag = -A_real
    # Left boundary condition u(xl,t) from time yb
    lside = l(yb+np.arange(0,n)*k)
    # Right boundary condition u(xr,t) from time tb
    rside = r(yb+np.arange(0,n)*k)
    # Initial conditions
    W_real = np.zeros((m, n))
    W_imag = np.zeros((m, n))
    W_real[:,0] = f_real(xl + np.arange(0,m)*h)
    W_imag[:,0] = f_imag(xl + np.arange(0,m)*h)
    # Solving for imaginary and real part
    for j in range(0,n-1):
        b_cond = np.concatenate(([lside[j]], np.zeros(m-2),[rside[j]]))
        W_real[:,j+1] =  W_real[:,j] + A_real.dot(W_imag[:,j]) - sigma*b_cond
        W_imag[:,j+1] =  W_imag[:,j] + A_imag.dot(W_real[:,j]) + sigma*b_cond
    return np.vstack([lside, W_real, rside]), np.vstack([lside, W_imag, rside])


xl, xr, yb, yt, M, N, xc, wl, s = (-9, 5, 0, 4, 35, 1000, -5, 4.0, 3.0)
W_real, W_imag = schroedinger1D(xl, xr, yb, yt, M, N, xc, wl, s)
[X, T] = np.meshgrid(np.linspace(xl, xr, M+1), np.linspace(yb, yt,N))
# Plot results
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.set_xlabel("$x$", fontsize=20)
ax.set_ylabel("$t$", fontsize=20)
ax.set_zlabel("$\Psi(x,t)$", fontsize=20)
print(X.shape)
print(T.shape)
print(W_real.T.shape)
surface = ax.plot_surface(X, T, W_real.T, cmap=cm.jet, linewidth=0, 
                          antialiased=True, rstride=10, cstride=10)
fig.colorbar(surface)
plt.tight_layout()
plt.show()

The output has the right shape for X, T and W_real.T (1000,36), however the surface has, as far as I understand, wrong assigned colors. I was expecting that the color varies in Z axis, but here I can't tell what is measuring:

This may be a real problem, a bug, or some problem with the input data; nobody can know without a minimal reproducible example of the issue. — ImportanceOfBeingErnest
– ImportanceOfBeingErnest, Commented Apr 29, 2017 at 9:32

ImportanceOfBeingErnest · Accepted Answer · 2017-04-29 14:22:33Z

2

It may be easier to understand what is happening when reducing the number of grid points

xl, xr, yb, yt, M, N, xc, wl, s = (-9, 5, 0, 4, 10, 10, -5, 4.0, 3.0)

and also using rstride and cstride of 1.

In this case the plot may look something like this

Now it's easy to spot the problem: The frequency of oscillations in your wavefunction is larger than the resolution of the grid. That means one single patch of the surface plot may start at a very low value and go up to a very high value. In this case, its color can be anything since it's determined by the value of one single edge of the patch. (If the patch starts at high value and goes down to a low value, it redder than if it starts at low value and goes up to a high value.)

The only possible solution is to make the grid more dense than the frequency of oscillations. I.e. try to visualize a wavefuction which is more smoothly variing, or visualize only part of the wavefunction but with a denser grid.

answered Apr 29, 2017 at 14:22

ImportanceOfBeingErnest

342k61 gold badges739 silver badges771 bronze badges

Sign up to request clarification or add additional context in comments.

1 Comment

Alejandro Sazo Over a year ago

Got it, I changed the wavefunction to a smoother one and over a small domain (just to feel the changes) and I started to work with the grid and now I can see the colour variation correctly. Thanks

Collectives™ on Stack Overflow

Matplotlib cmap not working as expected

1 Answer 1

1 Comment

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

1 Comment

Your Answer

Sign up or log in

Post as a guest

Related