参数随状态变量变化的ODE求解实现方法及电容案例修改咨询

当ODE的参数随状态变量动态变化时，核心是在微分方程的计算函数内部，实时根据当前状态值计算动态参数。如果参数与状态变量存在隐式关系（比如你的例子中R依赖电流i，而i又依赖R），则需要先对每个状态值求解出当前的参数，再代入微分方程计算导数。

针对你给出的RC充电电路，当R满足 R(t) = R₁·sin(v(t)) + R₂·cos(i(t))（其中 i(t) = (V_S - v(t))/R(t)），可以按以下步骤修改代码：

步骤1：明确隐式关系，求解动态参数 #

由于R的表达式中包含i，而i又依赖R，这是一个隐式方程：
$$R = R_1 \cdot \sin(v) + R_2 \cdot \cos\left( \frac{V_S - v}{R} \right)$$
我们需要对每个当前的v值，求解这个方程得到对应的R。

步骤2：修改ODE求解代码 #

import numpy as np
from scipy.integrate import odeint
from scipy.optimize import root

# 定义求解动态R的函数
def solve_R(v, R1, R2, V_S):
    # 构造隐式方程：F(R) = R - R1*sin(v) - R2*cos( (V_S - v)/R ) = 0
    def F(R):
        if R == 0:
            return np.inf  # 避免除以0
        return R - R1 * np.sin(v) - R2 * np.cos( (V_S - v)/R )
    # 用初始猜测值（比如恒定R的初始值）求解
    res = root(F, x0=10)
    return res.x[0]

# 定义微分方程
def f(v, t, R1, R2, C, V_S):
    # 实时求解当前v对应的R
    R = solve_R(v, R1, R2, V_S)
    # 计算dv/dt
    dvdt = (V_S - v) / (R * C)
    return dvdt

# 设定参数
C = 500e-6  # Farad
R1 = 5
R2 = 15
V_S = 10
t = np.linspace(0, 0.04, 999)
v0 = 0

# 求解ODE
v = odeint(f, v0, t, args=(R1, R2, C, V_S))

# 计算每个时刻的R和i
R_list = []
i_list = []
for v_val in v:
    R_val = solve_R(v_val, R1, R2, V_S)
    i_val = (V_S - v_val)/R_val
    R_list.append(R_val)
    i_list.append(i_val)
R_array = np.array(R_list)
i_array = np.array(i_list)

关键说明 #

• 动态参数的计算完全嵌入到ODE的导数函数f中，确保每一步求解时都使用当前状态对应的参数值。
• 对于隐式参数关系，使用数值求解器（如scipy.optimize.root）来得到当前状态下的参数值，初始猜测值可以选择恒定参数的初始值，保证求解收敛。
• 如果动态参数是状态变量的显式函数（比如R(t)=R1*sin(v(t))），则无需额外求解，直接在f中计算即可，代码会更简洁。

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