Search code examples
pythonnumpyodeint

I don't know why I get the multiple graphs of W versus PT using the attached python code

I want to get three graphs (W vs. X, W vs. y, W vs PT). But I can get two proper graphs (W vs. X and W vs. y). Unfortunately, I finally got multiple graphs of W vs PT (green lines). I don't know how to handle it. Anything you could do for me would be highly appreciated.

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

def PBR(X, W):
    a = 9.8*10**(-5)
    y = (1-a*W)**0.5
    PH2 = PT0*(1.5-X)*y
    PB = PT0 * X * y
    PT = PT0 * (1 - X)*y
    r = -k * PT * PH2 / (1 + KB*PB + KT *PT)
    dXdW = -r/FT0
    return dXdW

W = np.linspace(0, 10000)
KT = 1.038
KB = 1.39
FT0 = 50
k = 0.00087
PT0 = 12
X0 = 0
a = 9.8*10**(-5)
y = (1-a*W)**0.5
PT = PT0 * (1 - X)*y
X = odeint(PBR, X0, W)
plt.plot(W, PT, 'g', linewidth=0.5)
plt.plot(W, X,'r', linewidth=3.0)
plt.plot(W, y,'b', linewidth=5.0)

enter image description here

Solution

  • Are you looking for output like this. If yes - small issue in your script. Look at line marked as #**. You want to multiply with first dimension of X array.

    enter image description here

    import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

def PBR(X, W):
    a = 9.8*10**(-5)
    y = (1-a*W)**0.5
    PH2 = PT0*(1.5-X)*y
    PB = PT0 * X * y
    PT = PT0 * (1 - X)*y
    r = -k * PT * PH2 / (1 + KB*PB + KT *PT)
    dXdW = -r/FT0
    return dXdW

W = np.linspace(0, 10000)
KT = 1.038
KB = 1.39
FT0 = 50
k = 0.00087
PT0 = 12
X0 = 0
a = 9.8*10**(-5)
y = (1-a*W)**0.5
X = odeint(PBR, X0, W)
PT = PT0 * np.multiply(1 - X[:,0],y) #**
print(PT)
plt.plot(W, PT, 'g', linewidth=0.5)
plt.plot(W, X,'r', linewidth=3.0)
plt.plot(W, y,'b', linewidth=5.0)