0

I'm trying to plot this power law equation: Power law equation

This is my code. However, when n < 1, the plot only plots one side of the equation (the real side not the imaginary side). How do I plot the imaginary side (when n < 1)?

import numpy as np
from matplotlib import pyplot as plt

n = 0.4
k = 0.479862
p = 9470
r = np.linspace(-0.0127, 0.0127, 500)
d = 0.0254
u = (n/n+1)*(p/(2*k))** (1/n)*(d**((n+1)/n) - r**((n+1)/n))
plt.xlabel('Radius m')
plt.ylabel('Velocity l/min')
plt.plot(r, u)

n = 0.4 plot

n = 0.4 plot

It should look similar to this n=1 or newtonian plot

n=1 or newtonian plot

JohanC
  • 71,591
  • 8
  • 33
  • 66

1 Answers1

0

The problem has nothing to do with n being smaller than 1 or not. The problem arises whenever (n+1)/n isn't a perfect integer. Then, r**((n+1)/n) isn't well-defined for floating point numbers, for negative r, as explained e.g. in this post. Just try n = 2, and you'll notice it also doesn't have a value for any negative r.

Please also note that your python formula starts with (n/n+1) while your maths formula starts with (n/(n+1)).

To draw a symmetric curve, you might just copy and mirror the positive part.

import numpy as np
from matplotlib import pyplot as plt

n = 0.4
k = 0.479862
p = 9470
r = np.linspace(0, 0.0127, 100)
d = 0.0254
u = (n/(n+1))*(p/(2*k))** (1/n)*(d**((n+1)/n) - r**((n+1)/n))

plt.xlabel('Radius m')
plt.ylabel('Velocity l/min')
plt.plot(np.concatenate([-r[::-1], r]), np.concatenate ([u[::-1], u]))
plt.show()

Here is the resulting curve with the first factor of the formula adapted (n/(n+1)), for n=0.4. resulting plot

Alternatively, you can convert everything to complex numbers, and plot their real and the imaginary part separately. Therefore, you convert r using r.astype(np.complex). And np.imag(u) to get the imaginary part of u.

import numpy as np
from matplotlib import pyplot as plt

n = 0.4
k = 0.479862
p = 9470
r = np.linspace(-0.02, 0.02, 500)
r = r.astype(np.complex)
d = 0.0254
u = n / (n + 1) * (p / (2 * k)) ** (1 / n) * (d ** ((n + 1) / n) - r ** ((n + 1) / n))

plt.xlabel('Radius m')
plt.ylabel('Velocity l/min')
plt.plot(r, u, color='dodgerblue', label='real part')
plt.plot(r, np.imag(u), color='crimson', label='imaginary part')
plt.legend()
plt.show()

plot with real and imaginary part

JohanC
  • 71,591
  • 8
  • 33
  • 66