绘制在Python

Question

使用matplotlib条件函数

我试图创建这个算法的行或散点图，它给我的错误

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)) 

矢量化本质上是一个为循环，因此不会加速你的代码。

如果你能够使用熊猫然后我想你也可以使用apply功能尝试：

import pandas as pd 

xd = pd.Series(xr) 

yr = xd.apply(lambda x: P(x)) 

plt.plot(xr, yr) 

我相信你得到的特定错误的原因是因为你的条件实际上是在询问是否数组xr大于特定数字。阵列中的一些项目，有些不是，所以结果不明确。错误消息是问你是否想要问题是“是在xr中的任何东西大于这个数字”OR“是所有在xr大于这个数字”。

注意：您应该返回np.nan而不是第一个条件中的字符串。另外，函数integ.quad需要可调用的东西作为第一个参数。