I am using solve_ivp to solve a series of integrated differential equations. The below code is a simplified representation of a much larger code. I would like to be able to extract and plot/print 'over time', the evolving values of intermediary values, (eg) A_value vs t.
I have tried adding a counter, and adding the current value of 'A_value' to a list every iteration through A(t, pools), and this does work - but the problem is in the full version of this model, this becomes excessively large - I run into storage space/size issues (200+ intermediary calculations I would like to follow, my storage array size is upwards of a million per variable). I'm trying to find out if there is a way to do this more efficiently. I would like something that will allow the final plot A_value vs t.
I also tried this solution: Is there a way to store intermediate values when using solve_ivp? , but I really need all of the intermediary variables within one method (due to size), and I can't quite recreate their plotting results (does it work?).
(Other notes: I know that 'time' doesn't do anything in my A(t, pools) method, I was just re-creating the structure of the bigger model, where t is a required variable in that method. I am sure I am also using 'global' incorrectly, but was just trying anything...)
Any help or guidance would be greatly appreciated! How can I save the data and plot A_value vs t, without using a counter. I'm still a python 'newbie', and trying to learn by doing. Thanks in advance
"""
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
tspan = (0, 10) #set the start and stop time
pools0 = [1, 2] #set the initial values
def A(t, pools):
global A_value
global B_value
global C_value
global D_value
global time
time = t
if (time ==0):
cows = 1.0
C1, C2 = pools #the pools get saved in here
A_value = 2*C1 + 3*C2 - C1**2
B_value = 1*C1 + 4*C2
C_value = 1*C1 + 3*C2 - C2**2
D_value = 2*C1 + 5*C2
print(time, A_value, B_value, C_value, D_value) # essentially, how this prints, I want to be able to save/store this data
return [time, A_value, B_value, C_value, D_value]
time, A_value, B_value, C_value, D_value = A(t, pools)
def model_integ(t, pools):
global A_value
global B_value
global C_value
global D_value
global time
time, A_value, B_value, C_value, D_value = A(t, pools)
dC1 = A_value + B_value
dC2 = C_value + D_value
return [dC1, dC2] # return the differential
pools = solve_ivp(model_integ, tspan, pools0)
#print(A(pools.y[0], pools.y[1]) )
C1, C2 = pools.y
t = pools.t
print(t, C1, C2) # this works fine, pulling the t, C1, C2 values from model_integ() for printing
print(time, A_value, B_value, C_value, D_value) # this only prints out their FINAL values
# C1 pool vs t, storage variables
#---------------------------------------------------------------
plt.scatter(t,C1)
plt.xlabel("t, d") # adding labels and title on the plot
plt.ylabel("C1 pools") # adding labels and title on the plot
#plt.ylim(0.0,2500)
plt.xlim(0,10)
plt.title("C1/C2 change with time") # adding labels and title on the plot
plt.tight_layout() # this takes care of labels overflowing the plot. It's good to use it all the time.
plot_file = "cute_little_plot.png" # output file to save plot
plt.savefig(plot_file) # output file to save plot
plt.show()
plt.close()
# A_value vs t, storage variables
#---------------------------------------------------------------
plt.scatter(t,A_value)
plt.xlabel("t, d") # adding labels and title on the plot
plt.ylabel("A_value") # adding labels and title on the plot
#plt.ylim(0.0,2500)
plt.xlim(0,10)
plt.title("A_value change with time") # adding labels and title on the plot
plt.tight_layout() # this takes care of labels overflowing the plot. It's good to use it all the time.
plot_file = "Avalue_plot.png" # output file to save plot
plt.savefig(plot_file) # output file to save plot
plt.show()
plt.close()
"""
