用RK 45方法求解耦合方程组

Question

，我试图求解耦合微分方程组，但我得到了以下误差，任何一个都能解决这个问题。

dY[1] = (1/M*C_p)((m_c*C_p*(Y[0]-Y[1]))-(U_l*(A_s_2)*(Y[1]-Ta))) TypeError: 'float' object is not callable

import math
import scipy.integrate as spi
import numpy as np
import pandas as pd
def odes(t,Y):                                                                           
    m_c=0.9                                          
    C_p =4200
    Tc_O=91                                                                                   
    U_l =0.8                                          
    D_st=2                                         
    L_s=1
    rho=1000
    V=0.5
    M=(rho*V)/3                                            
    A_s_1=((math.pi*pow(D_st,2))/4)+((math.pi*D_st*L_s)/3)
    A_s_2=((math.pi*D_st*L_s))/3
    A_s_3=((math.pi*pow(D_st,2))/4)+((math.pi*D_st*L_s)/3)
    Ta =20 
    dY = np.zeros((3))
    dY[0] = (1/M*C_p)*((m_c*C_p*(Tc_O-Y[0]))-(U_l*(A_s_1)*(Y[0]-Ta)))
    dY[1] = (1/M*C_p)((m_c*C_p*(Y[0]-Y[1]))-(U_l*(A_s_2)*(Y[1]-Ta)))
    dY[2] = (1/M*C_p)((m_c*C_p*(Y[1]-Y[2]))-(U_l*(A_s_3)*(Y[2]-Ta)))
    return dY
t_start, t_end = 0, 3600.0
Y = np.array([91,89,75]); 
Yres = spi.solve_ivp(odes, [t_start, t_end], Y, method='RK45', max_step=60)
#---- Results ---------------------
yy = pd.DataFrame(Yres.y).T
tt = np.linspace(t_start,t_end,yy.shape[0])
print(yy)

Answer 1

问题在于：

dY[1] = (1/M*C_p)((m_c*C_p*(Y[0]-Y[1]))-(U_l*(A_s_2)*(Y[1]-Ta)))
                 ^
                 |
              No multiplication sign!

Python不像数学。如果两组括号相邻，没有运算符，这并不意味着乘法。它的意思是‘呼叫功能’

所以如果你这么做

print((3)(4))

您将得到以下错误：

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
TypeError: 'int' object is not callable

在这条线下面，你也有同样的问题。

Answer 2

你可以让它变得更干净，例如：

def odes(t, Y):
    Tc_O = 91
    D_st = 2
    L_s = 1
    A_s_1 = A_s_3 = ((math.pi * pow(D_st, 2)) / 4) + ((math.pi * D_st * L_s) / 3)
    A_s_2 = ((math.pi * D_st * L_s)) / 3
    dY = np.zeros((3))
    dY[0] = calc(Tc_O, Y[0], A_s_1)
    dY[1] = calc(Y[0], Y[1], A_s_2)
    dY[2] = calc(Y[1], Y[2], A_s_3)
    return dY


def calc(Y1, Y2, A_s):
    m_c = 0.9
    C_p = 4200
    U_l = 0.8
    rho = 1000
    V = 0.5
    M = (rho * V) / 3
    Ta = 20
    x = (1 / M * C_p) * ((m_c * C_p * (Y1 - Y2)) - (U_l * (A_s) * (Y2 - Ta)))
    return x


t_start, t_end = 0, 3600.0
Y = np.array([91, 89, 75])
Yres = spi.solve_ivp(odes, [t_start, t_end], Y, method="RK45", max_step=60)
# ---- Results ---------------------
yy = pd.DataFrame(Yres.y).T
tt = np.linspace(t_start, t_end, yy.shape[0])
print(yy)