Harrison Harrison - 12 days ago 3
Python Question

Plotting a conditional function using matplotlib in Python

I'm trying to create a line or scatter plot of this algorithm and it gives me the error


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()
."


I have looked up possible solutions to this error but I don't think any apply to me.

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


Is there a way to plot this function without separating each condition?

Answer

Numpy has a vectorize function that might be what you are looking for.

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

Vectorize is essentially a for loop, so wont speed up your code.

If you are able to use Pandas then I reckon you could also try using the apply function:

import pandas as pd

xd = pd.Series(xr)

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

plt.plot(xr, yr)

I believe the reason you were getting that specific error was because your condition was in effect asking if the array xr was greater than a specific number. Some items in the array were, some were not, so the result was ambiguous. The error message was asking if you wanted the question to be "is anything in xr greater than this number" OR " is everything in xr greater than this number".

Note: you should return np.nan rather than a string in your first condition. Also, the function integ.quad needs something that is callable as the first argument.

Comments