在Python中使用matplotlib绘制条件函数

Question

我试着画出这个算法的一条线或散点图，它给了我错误

Traceback (most recent call last): File "/Users/itstest/Documents/workspace/Practice/src/PlutoModel.py", line 73, in module plt.plot(xr, P(xr)) File "/Users/itstest/Documents/workspace/Practice/src/PlutoModel.py", line 55, in P if x > r: ValueError: "The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()."

我已经找到了解决这个错误的可能办法，但我不认为有什么适合我的。在

import numpy as np
import scipy.integrate as integ
import matplotlib.pyplot as plt

rho = 1860
rhos = 250 #Assuming Nitrogen Ice
rhom = 1000 #Assuming Water
rhoc = 3500 #Assuming a mix of Olivine and Pyroxene

def rhof(x):
    if x > r:
        return "Point is outside of the planet"
    elif x < r and x > rm:
        return rhos
    elif x < rm and x > rc:
        return rhom
    else:
        return rhoc

r = 1.187e6
rc = 8.5e5 #Hypothesized
rm = 9.37e5 #Estimated based on crustal thickness of 250 km

Ts = 44
B = 0.15
G = 6.67e-11

m = 1.303e22
mc = (4*np.pi*rhoc*rc**3)/3
mm = (4*np.pi*rhom*((rm**3) - (rc**3)))/3
ms = (4*np.pi*rhos*((r**3) - (rm**3)))/3

Ic = 0.4*mc*rc**2
Im = 0.4*mm*rm**2
Is = 0.4*ms*r**2
Itot = Is + Im + Ic

def gi(x):
    if x == r:
        return G*m/r**2
    elif x > r:
        return "Point is outside of the planet"
    elif x > rc and x < rm:
        return (G*mc/rc**2) + (G*mm/((x**2) - (rc**2)))
    elif x > rm and x < r:
        return (G*mc/rc**2) + (G*mm/((rm**2) - (rc**2))) + (G*ms/((x**2) - (rm**2)))
    else:
        return G*((3*rhoc)/4*np.pi*x**3)/x**2   

def Psmb(z):
    return rhos*G*(4.0/3.0)*np.pi*(1/z**2)*(rhom*(rm**3) + rhos*(z - rm**3))
def Pmcb(z):
    return rhom*G*(4.0/3.0)*np.pi*(1/z**2)*(rhoc*(rc**3) + rhom*(z - rc**3))
def P(x):
    if x > r:
        return "The point is outside of the planet"
    elif x == r:
        return 1
    elif x > rm and x < r:
        return (integ.quad(1000*gi(x), x, r))[0]
    elif x == rm:
        return (integ.quad(Psmb, x, r))[0]
    elif x > rc and x < rm:
        return (integ.quad(1000*gi(x), x, rm) + P(rm))[0]
    elif x == rc:
        return (integ.quad(Pmcb, x, rm) + P(rm))[0]
    elif x < rc and x != 0:
        return (integ.quad(1000*gi(x), x, rc) + P(rc))[0]
    else:
        return ((2.0/3.0)*np.pi*G*(rhoc**2)*r**2)

xr = np.linspace(0, 1187000, 1000)
plt.plot(xr, P(xr))

print("Mass = " + str(m) + " kg")
print("Radius = " + str(r) + " m")
print("Density = " + str(rho) + " kg/m^3")
print("Moment of Inertia = " + str(Itot) + " kgm^2")
print("Mean Surface Temperature = " + str(Ts) + " K")
print("Magnetic Field = " + str(B) + " nT")
print("Surface Gravity = " + str(gi(r)) + " m/sec^2")
print("Pressure at Surface = " + str(P(r)) + " Pa")
print("Pressure at Crust-Mantle Boundary = " + str(P(rm)) + " Pa")
print("Pressure at Mantle-Core Boundary = " + str(P(rc)) + " Pa")
print("Pressure at the Center = " + str(P(0)) + " Pa")

有没有办法在不分离每个条件的情况下绘制这个函数？在

Answer 1

Numpy有一个vectorize函数，这可能就是您要找的。在

pFunc = np.vectorize(P)
plt.plot(xr, pFunc(xr))

向量化本质上是一个for循环，因此不会加快代码的速度。在

如果您能够使用Pandas，那么我想您也可以尝试使用apply函数：

^{pr2}$

我相信你得到那个特定错误的原因是因为你的条件实际上是在询问数组xr是否大于一个特定的数字。数组中的某些项是，有些不是，因此结果不明确。错误消息询问您是否希望问题是“isanythisin xr greater this number”或“iseverythisin xr greater this number”。在

注意：你应该回来np.nan公司而不是第一个条件下的字符串。还有，函数集成电路需要一个可以调用的参数作为第一个参数。在