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

【问题标题】：Plotting a conditional function using matplotlib in Python在 Python 中使用 matplotlib 绘制条件函数
【发布时间】：2016-11-22 22:56:39
【问题描述】：

我正在尝试创建此算法的线图或散点图，但它给了我错误

Traceback（最近一次调用最后一次）：模块中的文件“/Users/itstest/Documents/workspace/Practice/src/PlutoModel.py”，第 73 行 plt.plot(xr, P(xr)) 文件“/Users/itstest/Documents/workspace/Practice/src/PlutoModel.py”，第 55 行，在 P 如果 x > r： ValueError: "具有多个元素的数组的真值不明确。使用a.any() 或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")

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

【问题讨论】：

总是显示完整的错误信息（回溯）。还有其他有用的信息 - 即。哪一行有问题。
什么是r？ np.linspace 创建一个数组。您的意思是要将您的条件应用于 x 中的每个数字吗？另外，matplotlib 应该如何绘制返回字符串？
更新了完整的回溯以及完整的代码。希望这能让它更清楚。 r 是一个常数。是的，我希望 x 中的每个值都通过算法运行并绘制该值。

标签： python if-statement matplotlib plot

【解决方案1】：

Numpy 有一个 vectorize 函数，这可能是您正在寻找的。p>

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

Vectorize 本质上是一个 for 循环，因此不会加速您的代码。

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

import pandas as pd

xd = pd.Series(xr)

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

plt.plot(xr, yr)

我相信您收到该特定错误的原因是因为您的条件实际上是在询问数组 xr 是否大于特定数字。数组中的一些项目是，一些不是，所以结果是模棱两可的。错误消息询问您是否希望问题是“xr 中的 任何事物 是否大于此数字”或“xr 中的 所有事物 是否大于此数字”。 p>

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

【讨论】：