用matplotlib绘制能量潜力

Question

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

这是返回每组坐标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

所有这些都可以逐个调试：

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

总的来说，你会得到这样的东西：