Необъяснимый выигрыш в моем коде при решении связанных дифференциальных уравнений с использованиемsolve_ivp() в Python ⇐ Python

[Anonymous]

Я пытаюсь решить два связанных дифференциальных уравнения в Python, используя функциюsolve:ivp() из scipy.integrate, но получаю необъяснимый выигрыш в одном из моих членов, который нарушает уравнение. Речь идет о движении пузырьков в акустическом поле, и я пытаюсь определить радиус, скорость изменения радиуса, положение (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(), но происходит то же самое.

Подробнее здесь: https://stackoverflow.com/questions/793 ... s-using-so