使用Scipy求解微分方程遇TypeError：float与CubicSpline无法做幂运算

问题 #

使用scipy库求解微分方程时出现TypeError，错误提示：TypeError: unsupported operand type(s) for ** or pow(): 'float' and 'CubicSpline'，原因是直接对float类型和CubicSpline对象执行了幂运算操作。

依赖包与数据 #

# 导入依赖包
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
from scipy.interpolate import CubicSpline
from scipy.integrate import odeint
from math import *

# 数据
data = {'day': [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93],
       'soil_temp': [18.15,17.5,19.1,20.3,19.75,17.7,15.2,15.45,14.3,12.45,12.75,14.55,16.55,18.3,19,19,18.8,17.45,17.15,17.4,19.9,19.85,21.4,22.05,21.75,19.9,21.9,23.45,24.65,24.4,25.1,24.75,25.2,25.45,25.75,26.35,26.5,24.8,24.55,25.95,26.35,23.9,22.2,21.2,21.9,23.4,25.45,25.75,25.25,25.65,26.4,25.7,25,26.1,27,26.75,26.95,26.55,25.9,26.2,27.15,28.25,27.95,27.25,26.5,27.45,27.55,27.8,28.4,28.8,28.05,25.05,25.15,25.45,25.3,22.95,22.6,25.1,25.95,26.3,26.55,26.25,27.15,27.75,28.2,25.45,25,25.1,25.15,25.15,26.05,26.2,27.45]}

# 创建DataFrame
df = pd.DataFrame(data)

报错代码 #

# 定义参数
alpha = 52.875
beta = 13.345
gamma = -1.44
delta = 2.29
constant = 60.589
g = 80.64
g2 = 1.04

# 定义模型
def model(a,t,om):
    # 三次样条插值
    day = df['day'].to_numpy()
    temp = df['soil_temp'].to_numpy()
    cubT = CubicSpline(day, temp, bc_type='natural',extrapolate=False)
    d_cubT = CubicSpline.derivative(cubT)

    # 模型参数计算
    bigfrac = (t/(((alpha-(beta*(-delta+(gamma*om))))/constant)*(g -(g2**cubT))))
    smallfrac = (t*(alpha-(beta*(-delta+(gamma*om)))/constant))

    # 导数方程
    k = (0.5**bigfrac) * log(0.5) * (((bigfrac - smallfrac - (g2**cubT) * log(g2) * d_cubT)) / (bigfrac**2))
    dherb_dt = -k*a
    return dherb_dt

# 初始条件
a0 = 4.271 

# 时间区间（天）
t = np.linspace(0, 90) # 90天

om = 0.3
y1 = odeint(model, a0, t, args=(om,))

# 绘图
plt.plot(t,y1, 'r-', linewidth = 2, label = 'om = 0.3')
plt.xlabel("天数")
plt.ylabel("除草剂浓度")
plt.legend()
plt.show()

错误堆栈 #

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Input In [1], in <cell line: 233>()
    230 t = np.linspace(0, 90) # 90 days
    232 om = 0.3
--> 233 y1 = odeint(model, a0, t, args=(om,))
    235 # Plot
    236 plt.plot(t,y1, 'r-', linewidth = 2, label = 'om = 0.3')

File ~\anaconda3\envs\agron893\lib\site-packages\scipy\integrate\_odepack_py.py:241, in odeint(func, y0, t, args, Dfun, col_deriv, full_output, ml, mu, rtol, atol, tcrit, h0, hmax, hmin, ixpr, mxstep, mxhnil, mxordn, mxords, printmessg, tfirst)
    239 t = copy(t)
    240 y0 = copy(y0)
--> 241 output = _odepack.odeint(func, y0, t, args, Dfun, col_deriv, ml, mu,
    242                          full_output, rtol, atol, tcrit, h0, hmax, hmin,
    243                          ixpr, mxstep, mxhnil, mxordn, mxords,
    244                          int(bool(tfirst)))
    245 if output[-1] < 0:
    246     warning_msg = _msgs[output[-1]] + " Run with full_output = 1 to get quantitative information."

Input In [1], in model(a, t, om)
    215 d_cubT = CubicSpline.derivative(cubT)
    217 # Model parameters
--> 218 bigfrac = (t/(((alpha-(beta*(-delta+(gamma*om))))/constant)*(g -(g2**cubT))))
    219 smallfrac = (t*(alpha-(beta*(-delta+(gamma*om)))/constant))
    221 # Derivative equation

TypeError: unsupported operand type(s) for ** or pow(): 'float' and 'CubicSpline'

解决方案 #

错误核心 #

CubicSpline是插值对象，不是数值数组，不能直接参与算术运算。必须调用该对象传入时间t，才能得到对应时刻的插值结果；导数对象同理，需要传入t获取对应时刻的导数值。另外，每次调用model函数都重新创建样条插值会降低效率，应该提前创建。

修改步骤 #

1. 提前创建样条插值对象：将cubT和d_cubT的定义移到model函数外部，避免重复计算。
2. 调用插值对象获取数值：在模型中使用cubT(t)获取对应时刻的土壤温度，d_cubT(t)获取温度的导数值，再参与运算。
3. 处理外插问题：原数据从第1天开始，若时间包含第0天，设置extrapolate=True允许外插，避免出现NaN。

修改后的完整代码 #

# 导入依赖包
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
from scipy.interpolate import CubicSpline
from scipy.integrate import odeint
from math import *

# 数据
data = {'day': [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93],
       'soil_temp': [18.15,17.5,19.1,20.3,19.75,17.7,15.2,15.45,14.3,12.45,12.75,14.55,16.55,18.3,19,19,18.8,17.45,17.15,17.4,19.9,19.85,21.4,22.05,21.75,19.9,21.9,23.45,24.65,24.4,25.1,24.75,25.2,25.45,25.75,26.35,26.5,24.8,24.55,25.95,26.35,23.9,22.2,21.2,21.9,23.4,25.45,25.75,25.25,25.65,26.4,25.7,25,26.1,27,26.75,26.95,26.55,25.9,26.2,27.15,28.25,27.95,27.25,26.5,27.45,27.55,27.8,28.4,28.8,28.05,25.05,25.15,25.45,25.3,22.95,22.6,25.1,25.95,26.3,26.55,26.25,27.15,27.75,28.2,25.45,25,25.1,25.15,25.15,26.05,26.2,27.45]}

# 创建DataFrame
df = pd.DataFrame(data)

# 定义参数
alpha = 52.875
beta = 13.345
gamma = -1.44
delta = 2.29
constant = 60.589
g = 80.64
g2 = 1.04

# 提前创建三次样条插值及导数对象
day = df['day'].to_numpy()
temp = df['soil_temp'].to_numpy()
cubT = CubicSpline(day, temp, bc_type='natural', extrapolate=True)  # 允许外插
d_cubT = cubT.derivative()  # 更简洁的导数获取方式

# 定义模型
def model(a, t, om):
    # 获取当前时刻的温度及导数数值
    T = cubT(t)
    dT_dt = d_cubT(t)
    
    # 预计算常量部分，避免重复计算
    const_term = alpha - beta * (-delta + gamma * om)
    const_over_constant = const_term / constant
    
    # 模型参数计算
    bigfrac = t / (const_over_constant * (g - (g2 ** T)))
    smallfrac = t * const_over_constant

    # 导数方程
    g2_pow_T = g2 ** T
    numerator = bigfrac - smallfrac - g2_pow_T * log(g2) * dT_dt
    k = (0.5 ** bigfrac) * log(0.5) * (numerator / (bigfrac ** 2))
    dherb_dt = -k * a
    return dherb_dt

# 初始条件
a0 = 4.271 

# 时间区间（天）
t = np.linspace(0, 90) # 90天

om = 0.3
y1 = odeint(model, a0, t, args=(om,))

# 绘图
plt.plot(t, y1, 'r-', linewidth=2, label='om = 0.3')
plt.xlabel("天数")
plt.ylabel("除草剂浓度")
plt.legend()
plt.show()

额外优化说明 #

• 将常量计算移到模型外部，减少每次调用时的重复计算。
• 使用cubT.derivative()替代CubicSpline.derivative(cubT)，代码更简洁。
• 拆分复杂表达式，提高代码可读性。

内容的提问来源于stack exchange，提问作者s_o_c_account