I have 15 ode equations need to be solved simultaneously and I want to solve them using solve_ivp.
There are each 5 states for T, co2, and q. The initial conditions are T=20, co2 = 0, q=0
I tried to separate them into 3 lists, one for T, one for co2, and one for q.
I am not sure how to resolve this bug and have worked on it for couple hours. Really appreciate your help!
import math
from scipy.integrate import solve_ivp
from bokeh.core.property.instance import Instance
from bokeh.io import save
from bokeh.layouts import column
from bokeh.model import Model
from bokeh.models import CustomJS, Slider, Callback
from bokeh.plotting import ColumnDataSource, figure, show
import numpy as np
# three plots , co2 as y and z as x
############### User generated - Slider initial value ###############
V= 100.0 # volume
r = 5.0
T = 20.0
c_co2_0 = 5.0 # concentration
episl_r = 0.3 # void
v0 = 2.0 # initial vilocity
############### --- Static Parameters --- ###############
b0 = 93.0 * (10**(-5))
deltH_0 = 95.3 # calculate b
Tw = -5.0 # room temperature
T0 = 353.15 # temeperature
t0 = .37 # heterogeneity constant, in paper is denoted as t_h0
alpha = 0.33
chi = 0.0
q_s0 = 3.40
R = 8.314
kT = 3.5*(10**3) #calculate rA
ρs = 880.0
deltH_co2 = 75.0 # calculate temeprature change
# ------------------ For Equation 4 : Enegergy Ballance --------------
ρg = 1.87 # ?
h = 13.8
Cp_g = 37.55 # J/molK
Cp_s = 1580.0 # J/molK
############### ----- Parameters depend on input ----- ###############
L = V / (math.pi * r**2)
deltZ = L / 5.0 # 5 boxes in total
p_co2 = R * T * c_co2_0
a_s = deltZ / r
theta = (1-episl_r) * ρs * Cp_s episl_r * ρg * Cp_g
# Equations are calclulated in order
def b(T):
b = ( b0 ** ( (deltH_0/ (R * T0) ) * (T0/T - 1) ) )
return b
def t_h(T):
return ( t0 alpha * (1 - T0 / T) )
def q_s(T):
return ( q_s0 ** ( chi * (1 - T / T0)) )
# Calculate rco2_n (not ode)
# change it to q
def R_co2(T, c_co2, q):
b_var = b(T)
t_var = t_h(T)
qs_var = q_s(T)
# print(qs_var)
r_co2 = kT * ( R * T * c_co2 * ( (1- ( (q / qs_var)**t_var) )**(1/t_var) ) - q / (b_var*qs_var) )
# print(r_co2)
return r_co2
# ODE Part
# Repetitive shortcut
# Equation 2
ener_balan_part1 = v0 * ρg* Cp_g
def ener_balan(theta, deltZ): # replace v0 * ρg* Cp_g / (theta * deltZ)
return(ener_balan_part1/ (theta*deltZ) )
def ener_balan2(episl_r):
return( (1-episl_r) * ρs * deltH_co2)
def ener_balan3(a_s, Tw, T0):
return (a_s * h *(Tw-T0))
# Equation 1 Mass Balance : find co2_n
def mass_balan(episl_r, deltZ):
return ( v0/ (episl_r * deltZ) )
def masss_balan2(episl_r, ρs):
return( (1-episl_r ) * ρs )
def deriv(t, y):
T_n, co2_n, q_n = y
# rco2_ first, rate of generation
T1 = -ener_balan(theta, deltZ) * T_n ener_balan(theta, deltZ) * T0 ener_balan2(episl_r)* (R_co2(T_n, co2_n, q_n)) ener_balan3(a_s, Tw, T0)
co2_1 = -mass_balan(episl_r, deltZ) * co2_n mass_balan(episl_r, deltZ) * c_co2_0 - (R_co2(T_n, co2_n, q_n)) * masss_balan2(episl_r, ρs)
q_1 = R_co2(T_n, co2_n, q_n)
T2 = -ener_balan(theta, deltZ) * T_n ener_balan(theta, deltZ) * T0 ener_balan2(episl_r)* (R_co2(T_n, co2_n, q_n)) ener_balan3(a_s, Tw, T0)
co2_2 = -mass_balan(episl_r, deltZ) * co2_n mass_balan(episl_r, deltZ) * c_co2_0 - (R_co2(T_n, co2_n, q_n)) * masss_balan2(episl_r, ρs)
q_2 = R_co2(T_n, co2_n, q_n)
T3 = -ener_balan(theta, deltZ) * T_n ener_balan(theta, deltZ) * T0 ener_balan2(episl_r)* (R_co2(T_n, co2_n, q_n)) ener_balan3(a_s, Tw, T0)
co2_3 = -mass_balan(episl_r, deltZ) * co2_n mass_balan(episl_r, deltZ) * c_co2_0 - (R_co2(T_n, co2_n, q_n)) * masss_balan2(episl_r, ρs)
q_3 = R_co2(T_n, co2_n, q_n)
T4 = -ener_balan(theta, deltZ) * T_n ener_balan(theta, deltZ) * T0 ener_balan2(episl_r)* (R_co2(T_n, co2_n, q_n)) ener_balan3(a_s, Tw, T0)
co2_4 = -mass_balan(episl_r, deltZ) * co2_n mass_balan(episl_r, deltZ) * c_co2_0 - (R_co2(T_n, co2_n, q_n)) * masss_balan2(episl_r, ρs)
q_4 = R_co2(T_n, co2_n, q_n)
T5 = -ener_balan(theta, deltZ) * T_n ener_balan(theta, deltZ) * T0 ener_balan2(episl_r)* (R_co2(T_n, co2_n, q_n)) ener_balan3(a_s, Tw, T0)
co2_5 = -mass_balan(episl_r, deltZ) * co2_n mass_balan(episl_r, deltZ) * c_co2_0 - (R_co2(T_n, co2_n, q_n)) * masss_balan2(episl_r, ρs)
q_5 = R_co2(T_n, co2_n, q_n)
T_ls = np.array([T1, T2, T3, T4, T5])
co2_ls = np.array([co2_1, co2_2, co2_3, co2_4, co2_5])
q_ls = np.array([q_1, q_2, q_3, q_4, q_5])
return T_ls, co2_ls, q_ls
t0, tf = 0, 10
############# initial condition
T_initial = 20
c_co2_0 = 0
q0 = 0
init_cond = np.array([20, 0, 0])
N=5
soln = solve_ivp(deriv, (t0, tf), init_cond)
Here is the error message
helper.py:293: RuntimeWarning: divide by zero encountered in double_scalars
r_co2 = kT * ( R * T * c_co2 * ( (1- ( (q / qs_var)**t_var) )**(1/t_var) ) - q / (b_var*qs_var) )
Traceback (most recent call last):
File "helper.py", line 350, in <module>
soln = solve_ivp(deriv, (t0, tf), init_cond)
File "/Users/cocochen/.local/share/virtualenvs/py-HkKPxrQC/lib/python3.8/site-packages/scipy/integrate/_ivp/ivp.py", line 546, in solve_ivp
solver = method(fun, t0, y0, tf, vectorized=vectorized, **options)
File "/Users/cocochen/.local/share/virtualenvs/py-HkKPxrQC/lib/python3.8/site-packages/scipy/integrate/_ivp/rk.py", line 96, in __init__
self.h_abs = select_initial_step(
File "/Users/cocochen/.local/share/virtualenvs/py-HkKPxrQC/lib/python3.8/site-packages/scipy/integrate/_ivp/common.py", line 104, in select_initial_step
d1 = norm(f0 / scale)
ValueError: operands could not be broadcast together with shapes (3,5) (3,)
CodePudding user response：
I am giving an example based answer as I recreated the exact error message. So what is the problem here is that you are not following Numpy broadcasting rules. Which basically says arrays can be broadcasted (given certain operation) if their dimensions are same or equal to 1.
Lets see an example below with error & solution:
import numpy as np
array_1 = np.arange(15).reshape(3,5)
array_2 = np.arange(3)array_1 has shape (3, 5) array_2 has shape (3,)
If you trying to add them (internally it will be broadcasted) using below code
array_1 array_2
This will be the error
ValueError: operands could not be broadcast together with shapes (3,5) (3,)
As we are not following above mentioned rules. To resolve this we need to reshape array_2 like this
array_2 = array_2.reshape(-1,1)
Now if you check, array_2 shape is
(3, 1)
And now if we try to add both arrays
array_1 array_2
It will work because now both of the rules are satisfied, both array either have dimensions or have dimension value as 1 (array_2 shape (3,1))
CodePudding user response：
https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html
The fun argument to solve_ivp is specified as
The calling signature is fun(t, y). Here t is a scalar, and there are two options for the ndarray y:
It can either have shape (n,); then fun must return array_like with shape (n,).
y is init_cond = np.array([20, 0, 0]), shape (3,)
That means deriv must return a (3,) result. What is
deriv(0, init_cond)
I'm not going download and run your code, but the end of deriv is
T_ls = np.array([T1, T2, T3, T4, T5])
co2_ls = np.array([co2_1, co2_2, co2_3, co2_4, co2_5])
q_ls = np.array([q_1, q_2, q_3, q_4, q_5])
return T_ls, co2_ls, q_ls
That creates 3 shape (5,) array, and returns them as tuple. That makes a (3,5) array.
My guess is that the error point
norm(f0 / scale)
f0 is from a trial func run, and hence (3,5), while scale is (3,) derived from the init_cond.
Failure to read or follow the specifications is a common cause of problems with these scipy solve functions (also the optimize and integrate).