Home > database >  Plot a well known double well surface
Plot a well known double well surface

Time:06-15

I am trying to plot a well known energy landscape from physics, the Muller Brown potential. Taken from the literature (enter image description here enter image description here The authors' plot:

With my contour plot, however, I cannot see the two wells, it looks as if it was just a single Gaussian. Am I doing something wrong?

import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mp

# Parameters of Muller Brown Potential (cf. Bonfati & Cob 2017)
A = np.array([-200, -100, -170, 15])
a = np.array([-1, -1, -6.5, 0.7])
b = np.array([0, 0, 11, 0.6])
c = np.array([-10, -10, -6.5, 0.7])
x_m = np.array([1, 0, -0.5, -1])
y_m = np.array([0, 0.5, 1.5, 1])

x = np.linspace(-1.5, 1, 1000)
y = np.linspace(-0.5, 2, 1000)
XX, YY = np.meshgrid(x, y)


Z =A[0]*np.exp( a[0]*(XX-x_m[0])**2   b[0]*(XX-x_m[0])*(YY-y_m[0])   c[0]*(YY-y_m[0])**2 )
   A[1]*np.exp( a[1]*(XX-x_m[1])**2   b[1]*(XX-x_m[1])*(YY-y_m[1])   c[1]*(YY-y_m[1])**2 )
   A[2]*np.exp( a[2]*(XX-x_m[2])**2   b[2]*(XX-x_m[2])*(YY-y_m[2])   c[2]*(YY-y_m[2])**2 )
   A[3]*np.exp( a[3]*(XX-x_m[3])**2   b[3]*(XX-x_m[3])*(YY-y_m[3])   c[3]*(YY-y_m[3])**2 )



fig, ax = plt.subplots()

c=ax.contourf(XX, YY, Z)
plt.colorbar(c)
ax.set_xlabel('x')
ax.set_ylabel('y')

enter image description here

edit: Setting the limits of the plot region to the ones the authors does not seem to help. If I use

x = np.linspace(-2, 0, 1000)
y = np.linspace(0, 2, 1000)

I see this: enter image description here

CodePudding user response:

You forgot the parantheses around the expression for Z so that you effectively evaluated only the first summand.

x = np.linspace(-1.5, 0.2, 1000)
y = np.linspace(0, 1.9, 1000)
XX, YY = np.meshgrid(x, y)

Z = (A[0]*np.exp( a[0]*(XX-x_m[0])**2   b[0]*(XX-x_m[0])*(YY-y_m[0])   c[0]*(YY-y_m[0])**2 )
     A[1]*np.exp( a[1]*(XX-x_m[1])**2   b[1]*(XX-x_m[1])*(YY-y_m[1])   c[1]*(YY-y_m[1])**2 )
     A[2]*np.exp( a[2]*(XX-x_m[2])**2   b[2]*(XX-x_m[2])*(YY-y_m[2])   c[2]*(YY-y_m[2])**2 )
     A[3]*np.exp( a[3]*(XX-x_m[3])**2   b[3]*(XX-x_m[3])*(YY-y_m[3])   c[3]*(YY-y_m[3])**2 ))

enter image description here

  • Related