使用Python solve_ivp求解气泡声场运动耦合微分方程时u_b项异常指数增长问题求助

我正在用Python的scipy.integrate.solve_ivp()函数求解气泡在声场中运动的耦合微分方程，但其中的气泡速度项u_b出现了无法解释的指数增长，直接导致求解器崩溃。我需要求解的变量包括气泡半径R、半径变化率R_dot、位置z以及气泡速度u_b，但u_b会迅速增长到完全不切实际的数值，让整个求解过程无法进行下去。

以下是我使用的代码：

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

def solver_ode(t, y):
    rho_l = 997.4  # Liquid density (kg/m^3)
    mu_l = 8.9e-4  # Dynamic viscosity (Pa·s)
    sigma = float(0.072)  # Surface tension (N/m)
    P_l0 = 1e5  # Far-field pressure (Pa)
    gamma = 1.4  # Polytropic index
    R0 = 0.00025  # Initial bubble radius (m)
    f = 600  # Acoustic frequency (Hz)
    omega = 2 * np.pi * f  # Angular frequency (rad/s)
    c = 1500  # Speed of sound in medium (m/s)
    lamda = c / f  # Acoustic wavelength (m)
    k = 2 * np.pi / lamda  # Wavenumber (1/m)
    g = 9.81  # Gravity (m/s^2)
    dPa = 60000  # Acoustic amplitude (Pa)
    rho_b = 1.2  # Bubble density (kg/m^3)

    R, R_dot, z, u_b = y  # Unpacking condition


    # Gas pressure

    P_v = 0
    Pg0 = ((2 * sigma) / R0) + P_l0 - P_v
    Pg = Pg0 * (R0 / R) ** (3 * gamma)

    # Acoustic pressure
    P_ac = dPa * np.cos(omega * t) * np.sin(k*z)
    P_inf = P_ac + P_l0
    P_b = Pg + P_v

    # Radial acceleration
    R_ddot = (1 / (rho_l * R)) * (P_b - P_inf) - \
             (3 / (2 * R)) * (R_dot ** 2) - \
             ((4 * mu_l) / (rho_l * (R ** 2))) * R_dot - \
             (2 * sigma / (rho_l * (R ** 2)))
    
    V = (4 / 3) * np.pi * R ** 3
    V_dot = 4 * np.pi * R**2 * R_dot
    m_b = rho_b * ((4 / 3) * np.pi * R0 ** 3)    
    A = 4 * np.pi * R ** 2
    

    # Forces
    u_l = ((-k * dPa) / (omega * rho_l)) * np.sin(omega * t) * np.cos(2 * np.pi * k*z)
    u_l_dot = ((-k * dPa) / (rho_l)) * np.cos(omega * t) * np.cos(2 * np.pi * k*z)

    F_bj = -V * (dPa * k * np.cos(omega * t) * np.cos(2 * np.pi * k*z))
    V_dot = 0

    u_b_dot = (1/(m_b + 0.5*rho_l*V))*(F_bj - 0.5*rho_l*V_dot*(u_b - u_l) +
                0.5*rho_l*V*(-k*(dPa/(omega*rho_l)*(omega*np.cos(omega*t)*np.cos(k*z) - k*np.sin(omega*t)*np.sin(k*z)*u_b))) - 0.5*rho_l*(u_b - u_l)**2*A*(27*(abs(mu_l/(2*rho_l*(u_b - u_l)*R0)))**0.78) + V*(rho_l - rho_b)*g)


    return [R_dot, R_ddot, u_b, u_b_dot]


# Initial conditions: [R, R_dot, pos, u_b]
y0 = [0.00025, 0, 0.001, 0.00001]  # Start at a quarter wavelength

t_span = (0, 1)  # Time span (s)
t_eval_1 = np.linspace(0, 1, 1000)  # Time steps

# Solve the system

solution = solve_ivp(
    solver_ode, t_span, y0,  method='RK45', max_step=0.0001,
    rtol=1e-6, atol=1e-9
)

print(solution)

提前感谢各位的帮助！

我已经尝试过调整参数，也试过简化方程——固定气泡半径只求解位置和速度，但问题依然存在，u_b还是会指数增长最终导致求解器崩溃。我也换了solve_ivp()提供的其他求解器，结果还是一样。

备注：内容来源于stack exchange，提问作者Lucas Leal Abadi