3D打印预测一阶积分微分方程数值求解问题求助

本人是工程师，刚学习Python。此前已使用solve_ivp求解出μ为常数时的常微分方程（ODE），现需求解μ随ϕ动态变化的一阶积分微分方程，编写的代码出现TypeError: only size-1 arrays can be converted to Python scalars错误，求助该方程的数值求解方法。

常数μ时的ODE及代码 #

常数定义：

kI_0 = 0.0015
μ = 3.74
z = np.linspace(0, 1.15, 100)

求解代码：

import numpy as np
import math as ma
import matplotlib.pyplot as plt
import scipy as sp
from scipy.integrate import odeint
from scipy.integrate import solve_ivp

# 需要用np.exp()否则无法运行
def dϕdt(t, ϕ, kI_0, μ, z):
    return kI_0*(1-ϕ)*np.exp(-μ*z)

kI_0 = 0.0015
μ = 3.74

# 初始条件设为100个0组成的向量，求解器需要数组而非标量
ϕ = np.zeros(100)

# z从0到1.15取100个点
z = np.linspace(0, 1.15, 100) 

# t_span不指定步长，让求解器自动选择合适步长
# 用args传递ODE函数的额外参数
sol_m1 = solve_ivp(dϕdt, y0=ϕ, t_span=(0, 3600), args=(kI_0, μ, z))

# 取y的最后一列所有行，得到最终时刻的转化率结果
plt.plot(z, sol_m1.y[:,-1])
plt.ylabel("转化率")
plt.xlabel("固液界面移动距离 (mm)")
plt.title("液相转化率随深度变化\n(数值解)")
print(sol_m1)

变量μ的积分微分方程及尝试代码 #

μ的表达式：

μ = μ_0*(1 - ϕ) + μ_infny*ϕ

常数定义：

μ_0 = 1
μ_infny = 5
kI_0 = 0.0015

尝试代码：

import numpy as np
import math as ma
import matplotlib.pyplot as plt
import sympy as sy
import scipy as sp
from scipy.integrate import odeint
from scipy.integrate import solve_ivp
from scipy.integrate import quad
from scipy.integrate import cumulative_trapezoid

μ_0 = 1
μ_infny = 5
kI_0 = 0.0015


def μ(ϕ, μ_0, μ_infny):
    return μ_0*(1 - ϕ) + μ_infny*ϕ 


# 初始条件设为100个0组成的向量
ϕ = np.zeros(100)

# z从0到1.15取100个点
z = np.linspace(0, 1.15, 100) 

# 定义积分函数
f = lambda z : -μ(ϕ, μ_0, μ_infny)
Int = quad(f, 0, 1.15)

# 需要用np.exp()否则无法运行
def dϕdt(t, ϕ, kI_0, Int):
    return kI_0*(1-ϕ)*np.exp(Int)


# 调用solve_ivp求解
sol_m1 = solve_ivp(dϕdt, y0=ϕ, t_span=(0, 3600), args=(kI_0, Int))
  

# 绘制结果
plt.plot(z, sol_m1.y[:,-1])
plt.ylabel("转化率")
plt.xlabel("固液界面移动距离 (mm)")
plt.title("液相转化率随深度变化\n(数值解)")
print(sol_m1)

报错信息：TypeError: only size-1 arrays can be converted to Python scalars（类型错误：仅长度为1的数组可转换为Python标量）

错误原因 #

1. quad返回值不匹配：quad返回的是包含积分结果和误差估计的元组，不是单纯的数值数组，直接传入np.exp()会触发类型错误。
2. 积分逻辑错误：μ随ϕ（随时间、z位置变化）动态改变，不能用初始ϕ提前计算一次全局积分，必须在每个时刻、每个z点实时计算累积积分。
3. 维度不兼容：μ对应每个z点的ϕ值，是数组类型，全局积分无法匹配每个z点的计算需求，导致维度冲突。

修正后的代码 #

假设积分微分方程形式为：$\frac{d\phi}{dt} = kI_0(1-\phi) \exp\left( -\int_0^z \mu(\phi(z',t)) dz' \right)$，修正代码如下：

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

μ_0 = 1
μ_infny = 5
kI_0 = 0.0015

# 计算μ的函数
def calc_mu(phi, mu0, mu_inf):
    return mu0*(1 - phi) + mu_inf*phi

# 定义ODE右侧函数
def dϕdt(t, phi, kI0, mu0, mu_inf, z):
    # 实时计算每个z点对应的μ
    mu = calc_mu(phi, mu0, mu_inf)
    # 用梯形法计算从0到每个z的累积积分
    integral = np.cumsum(mu * np.diff(z, prepend=0))
    # 返回每个z点的dphi/dt
    return kI0*(1 - phi)*np.exp(-integral)

# 初始条件：每个z点的初始phi为0
phi0 = np.zeros(100)
# z的网格点
z = np.linspace(0, 1.15, 100)

# 调用solve_ivp求解
sol = solve_ivp(dϕdt, y0=phi0, t_span=(0, 3600), args=(kI_0, μ_0, μ_infny, z), method='RK45')

# 绘制最终时刻的结果
plt.plot(z, sol.y[:, -1])
plt.ylabel("转化率")
plt.xlabel("固液界面移动距离 (mm)")
plt.title("液相转化率随深度变化\n(数值解，μ随ϕ动态变化)")
plt.show()
print(sol)

关键修正点 #

• 用np.cumsum结合np.diff实现累积积分，对每个z点计算从0到当前位置的积分，符合方程物理意义。
• 去掉错误的quad调用，改为在每个求解时刻实时计算μ和对应积分，适配μ随ϕ动态变化的特性。
• 全程使用numpy数组操作，保持维度一致，满足solve_ivp对向量输入的要求（y数组对应每个z点的ϕ值）。

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