繪製在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需要可調用的東西作爲第一個參數。