我有15个节点方程需要同时求解，我想用solve_ivp来求解它们。
对于T、co2和q，每一个都有5个状态。初始条件是T=20，co2 = 0，q=0
我试着把它们分成3个列表，一个是T，一个是CO2，一个是Q。
我不知道如何解决这个bug，已经工作了几个小时了。真的很感谢你的帮助！
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)以下是错误消息
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,)
2条答案
按热度按时间bsxbgnwa1#
我给出了一个基于示例的答案，因为我重新创建了确切的错误消息。所以这里的问题是，你没有遵循Numpy广播规则。基本上说，数组可以被广播（给定特定操作），如果它们的维数是相同或等于1。
让我们看看下面的错误和解决方案示例：
数组_1具有形体**（3，5）**数组_2具有形体（3，）
如果您尝试使用以下代码添加它们（将在内部广播）
这将是错误
因为我们没有遵循上面提到的规则。要解决这个问题，我们需要如下重新整形array_2
现在，如果您检查，array_2形状为
现在如果我们尝试将两个数组
它将工作，因为现在两个规则都满足，两个数组都有维度或维度值为1（array_2 shape（3，1））
blmhpbnm2#
https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html
solve_ivp的fun参数指定为y是init_cond = np.array([20, 0, 0])，形状为（3，）这意味着
deriv必须返回一个（3，）结果。我不会下载并运行您的代码，但
deriv的结尾是创建3个形状的（5，）数组，并将它们作为元组返回。这样就得到了一个（3，5）数组。
我的猜测是错误点
f0来自于func的试验运行，因此是（3，5），而scale是从init_cond导出的（3，）。未能阅读或遵循规范是这些
scipy求解函数（以及优化和集成）出现问题的常见原因。