extract values from function used by odeint scipy python

Question

I have the following script to calculate dRho using odeint.

P_r = 10e5
rho_r = 900
L = 750
H = 10
W = 150
A = H * W
V = A * L
fi = 0.17

k = 1.2e-13
c = 12.8e-9
mu = 2e-3

N = 50
dV = V/N
dx = L/N

P_in = P_r
rho_in = rho_r

P_w = 1e5    
rho_w = rho_r* np.exp(c*(P_w-P_r))

# init initial case
P = np.empty(N+1)*10e5
Q = np.ones(N+1)
out = np.empty(N+1)

P[0] = P_w
Q[0] = 0
out[0] = 0

def dRho(rho_y, t, N):

    P[1:N] = P_r + (1/c) * np.log(rho_y[1:N]/rho_r)
    P[N] = P_r + (1/c) * np.log(rho_y[N]/rho_r)


    Q[1:N] = (-A*k/mu)*((P[1-1:N-1] - P[1:N])/dx)
    Q[N] = (-A*k/mu)*((P[N]-P_r)/dx)


    out[1:N] = ((Q[1+1:N+1]*rho_y[1+1:N+1] - Q[1:N]*rho_y[1:N])/dV) 
    out[N] = 0

    return out

t0 = np.linspace(0,1e9, int(1e9/200))
rho0 = np.ones(N+1)*900
ti = time.time()
solve = odeint(dRho, rho0, t0, args=(N,))
plt.plot(t0,solve[:,1:len(rho0)], '-', label='dRho')
plt.legend(loc='upper right')
plt.show()

P and Q are calculated within the function dRho, they P acts and input for Q and both P, Q and rho_y act as input for out. The function returns "out". I can plot out without any issues, however, I am interested in plotting P and Q as well.

I have tried various approaches to achieve this like: Recalculating P and Q after the integration method but this increased the runtime of the script. So since the calculation is done within dRho I was wondering if and how I could access it from outside to plot it.

I have also tried to add P and Q together with rho0 as input for odeint but both P and Q were taken in the integration which resulted in wrong outcome when returned by the function.

A simplified version:

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
def dY(y, x):
    a = 0.001
    yin = 1
    C = 0.01
    N = 1
    dC = C/N
    b1 = 0
    y_diff = -np.copy(y)
    y_diff[0] += yin
    y_diff[1:] += y[:-1]
    print(y)
    return (a/dC)*y_diff+b1*dC

x = np.linspace(0,20,1000)
y0 = np.zeros(4)
res = odeint(dY, y0, x)
print(res)
plt.plot(x,res, '-')
plt.show()

in this simplified example, I would like to create an additional plot of ydiff.

Here another simple case:

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

def func(z,t):
    x, y=z
    xnew = x*2
    print(xnew)
    ynew = y*0.5
#     print y
    return [x, y]    

z0=[1,3]
t = np.linspace(0,10)
xx=odeint(func, z0, t)
plt.plot(t, xx[:,0],t,xx[:,1])
plt.show()

I am interested in plotting all xnew and ynew values.

Another example:

xarr = np.ones(4)
def dY(y, x):
    a = 0.001
    yin = 1
    C = 0.01
    N = 1
    dC = C/N
    b1 = 0
    xarr[0] = 0.25
    xarr[1:] = 2 
    mult = xarr*2
    out = mult * y
    print(mult)
    return out

x = np.linspace(0,20,1000)
y0 = np.zeros(4)+1.25
res = odeint(dY, y0, x)
dif = np.array([dY(y,x) for y in res])
print(dif)
plt.plot(x,res, '-')
plt.show()

I would like to plot the mult values against x

I think you would increase your chances of getting help by providing a [mcve] of the issue, i.e. we don't need that complex systems to get to the point of identifying the problem and solve it; all it does is making the whole thing much too complicated. — ImportanceOfBeingErnest, Oct 18 '17 at 14:51
Per the request of @ImportanceOfBeingErnest see edited post. — MD', Oct 18 '17 at 15:03
`y_diff` is changing for each iteration. Do you want to plot the last `y_diff` that is used? — ImportanceOfBeingErnest, Oct 18 '17 at 16:26
@hpaulj I have tried it but it unfortunately did not yield the correct output. — MD', Oct 18 '17 at 16:40
@ImportanceOfBeingErnest I want to plot all changes y_diff encounters during all time steps against x, where x = time steps — MD', Oct 18 '17 at 16:42
`y_diff` is an array. So you have one such array per timestep. You would want to plot n such arrays for n timesteps? How would that plot look like? — ImportanceOfBeingErnest, Oct 18 '17 at 16:45
@ImportanceOfBeingErnest I am not sure what your question exactly means. I have added another two simplified examples which hopefully elaborates what I am looking for. However in my original case, I am calculating the delta Rho which is dependent on Q and P. Rho changes based on those two changing, I would like to visualise the changes in Q and P in graphs for further studying. I am able to plot P and Q by simply hardcoding their equations after the odeint. — MD', Oct 18 '17 at 17:07
I know this is not an "answer". In general, when solving an ODE problem the "costly" part is to obtain an accurate solution. Once that is obtained, it is relatively "cheap" in terms of CPU time to recompute "P" and "Q". Consider whether you need to store them while solving the ODE or not based on performance. — Pierre de Buyl, Oct 18 '17 at 20:37
@PierredeBuyl Your remark is definitly valid and worth noting. I might even run both options side by side to study the differences based on performance. — MD', Oct 19 '17 at 21:19

score 4 · Accepted Answer · answered Oct 18 '17 at 17:08

The following could be what you want. You could store the intermediate values in a list and later plot that list. That would require to store the x values as well.

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

xs = []
yd = []

def dY(y, x):
    a = 0.001
    yin = 1
    C = 0.01
    N = 1
    dC = C/N
    b1 = 0
    y_diff = -np.copy(y)
    y_diff[0] += yin
    y_diff[1:] += y[:-1]
    xs.append(x)
    yd.append(y_diff)
    return (a/dC)*y_diff+b1*dC

x = np.linspace(0,20,1000)
y0 = np.zeros(4)
res = odeint(dY, y0, x)

plt.plot(x,res, '-')

plt.gca().set_prop_cycle(plt.rcParams['axes.prop_cycle'])
plt.plot(np.array(xs),np.array(yd), '-.')

plt.show()

Dotted lines are the respective y_diff values for the res solutions of the same color.

This is definitely the answer I am looking. At first sight, the implementation in my own code seems to output the right results. — MD', Oct 19 '17 at 21:17

extract values from function used by odeint scipy python

1 Answers1

Linked