使用 matplotlib 绘制能量势图

Question

我想绘制引力势能以突出其极值（两个天体周围的拉格朗日点）。

这里的函数是returns每组坐标x和y的势：

def gravitational_potential(M,m,R,x,y):
    G = 6.674*10**(-11)
    omega2 = G*(M+m)/(R**3)
    r = np.sqrt(x**2+y**2)
    r2 = R*m/(M+m)
    r1 = R-r2

    phi = -G*(M/abs(r-r1)+m/abs(r-r2))-1/2*omega2*(x**2+y**2)

    return phi

我想用meshgrid和plot_surface来绘制势的表面和等高线，但它不起作用。

我做错了什么？

PS：我设法用 WolframAlpha 绘制了潜力，所以我知道数学是有效的。

from mpl_toolkits.mplot3d import axes3d
import matplotlib.pyplot as plt
from matplotlib import cm
import numpy as np

def gravitational_potential(M,m,R,x,y):
    G = 6.674*10**(-11)
    omega2 = G*(M+m)/(R**3)
    r = np.sqrt(x**2+y**2)
    r2 = R*m/(M+m)
    r1 = R-r2

    phi = -G*(M/abs(r-r1)+m/abs(r-r2))-1/2*omega2*(x**2+y**2)

    return phi


fig = plt.figure()
ax = fig.gca(projection='3d')

X, Y = np.meshgrid(np.arange(-20, 20, 0.5), np.arange(-20, 20, 0.5))

M = 10
m = 1
R = 10
Z = gravitational_potential(M,m,R,X,Y)



ax.plot_surface(X, Y, Z, rstride=8, cstride=8, alpha=0.9)
cset = ax.contour(X, Y, Z, zdir='z', offset=-40, cmap=cm.coolwarm)
cset = ax.contour(X, Y, Z, zdir='x', offset=-20, cmap=cm.coolwarm)
cset = ax.contour(X, Y, Z, zdir='y', offset=20, cmap=cm.coolwarm)

ax.set_xlabel('X')
ax.set_xlim(-20, 20)
ax.set_ylabel('Y')
ax.set_ylim(-20, 20)
ax.set_zlabel('Z')
ax.set_zlim(-40, 40)


plt.show()

当我执行它时，我得到以下信息：

runfile('C:/Users/python/Google Drive/lagrangepoint_maths/potential/gravitational_potential.py', wdir='C:/Users/python/Google Drive/lagrangepoint_maths/potential')
C:/Users/python/Google Drive/lagrangepoint_maths/potential/gravitational_potential.py:13: RuntimeWarning: divide by zero encountered in divide
  phi = -G*(M/abs(r-r1)+m/abs(r-r2))-1/2*omega2*(x**2+y**2)

这不是我真正想要的。 Z有问题。我想要这样的东西：

from mpl_toolkits.mplot3d import axes3d
import matplotlib.pyplot as plt
from matplotlib import cm

fig = plt.figure()
ax = fig.gca(projection='3d')
X, Y, Z = axes3d.get_test_data(0.05)
ax.plot_surface(X, Y, Z, rstride=8, cstride=8, alpha=0.9)
cset = ax.contour(X, Y, Z, zdir='z', offset=-100, cmap=cm.coolwarm)
cset = ax.contour(X, Y, Z, zdir='x', offset=-40, cmap=cm.coolwarm)
cset = ax.contour(X, Y, Z, zdir='y', offset=40, cmap=cm.coolwarm)

ax.set_xlabel('X')
ax.set_xlim(-40, 40)
ax.set_ylabel('Y')
ax.set_ylim(-40, 40)
ax.set_zlabel('Z')
ax.set_zlim(-100, 100)

plt.show()

Answer 1

所有这些都是可以一一调试的东西：

如果分母小于分母，python 2 中的整数除法结果为 0。您可以 from __future__ import division 或更正您的代码以除以浮点数。
如果要显示介于 -2 x 10^-8 和 +2 x 10^-8 之间的数字，将 z_limits 设置为 -40 到 40 是没有用的。
如果要在绘图中显示小特征，则不应将绘图分辨率粗略设置为rstride=8, cstride=8。

总的来说你会得到这样的结果：

使用 matplotlib 绘制能量势图

Plot an energy potential with matplotlib

python

matplotlib

mplot3d