用 scipy 终止 ODE 积分的代数约束答案

【问题标题】：algebraic constraint to terminate ODE integration with scipy用 scipy 终止 ODE 积分的代数约束
【发布时间】：2014-07-28 15:46:12
【问题描述】：

我正在使用 Scipy 14.0 求解一个常微分方程组，该方程组描述了由于浮力而在静止流体中垂直（沿 z 方向）上升的气泡的动力学。特别是，我有一个方程将上升速度 U 表示为气泡半径 R 的函数，即 U=dz/dt=f(R)，还有一个方程将半径变化表示为 R 和 U 的函数，即 dR/dT =f(R,U)。下面代码中出现的所有其余部分都是材料属性。我想实现一些东西来解释对 z 的物理约束，这显然受到液体高度 H 的限制。因此，我实现了一种 z

提前致谢

埃米

from scipy.integrate import ode 

Db0 = 0.001 # init bubble radius
y0, t0 = [ Db0/2 , 0. ], 0. #init conditions
H = 1

def y_(t,y,g,p0,rho_g,mi_g,sig_g,H):
    R = y[0] 
    z = y[1]
    z_ = ( R**2 * g * rho_g ) / ( 3*mi_g )  #velocity
    R_ = ( R/3 * g * rho_g * z_ ) / ( p0 + rho_g*g*(H-z) + 4/3*sig_g/R ) #R dynamics    
    return [R_, z_]

def z_constraint(t,y):
    H = 1  #should rather be a variable..
    z = y[1] 
    if z >= H:
        flag = -1 
    else:
        flag = 0
    return flag

r = ode( y_ )
r.set_integrator('dopri5')
r.set_initial_value(y0, t0)
r.set_f_params(g, 5*1e5, 2000, 40, 0.31, H)
r.set_solout(z_constraint)

t1 = 6
dt = 0.1

while r.successful() and r.t < t1:
    r.integrate(r.t+dt)

【问题讨论】：

标签： python scipy

【解决方案1】：

您遇到了this issue。要使set_solout 正常工作，必须在set_integrator 之后、set_initial_value 之前调用它。如果您将此修改引入您的代码（并为g 设置一个值），集成将在z >= H 时终止，如您所愿。

要找到气泡到达表面的确切时间，您可以在积分被solout 终止后更改变量，并相对于z（而不是t）积分回@987654332 @。描述该技术的论文是M. Henon, Physica 5D, 412 (1982)；您可能还会发现 this discussion 很有帮助。这是一个非常简单的示例，在给定dy/dt = -y 的情况下，找到了y(t) = 0.5 的时间t：

import numpy as np
from scipy.integrate import ode

def f(t, y):
    """Exponential decay: dy/dt = -y."""
    return -y

def solout(t, y):
    if y[0] < 0.5:
        return -1
    else:
        return 0

y_initial = 1
t_initial = 0

r = ode(f).set_integrator('dopri5')
r.set_solout(solout)
r.set_initial_value(y_initial, t_initial)

# Integrate until solout constraint violated
r.integrate(2)

# New system with t as independent variable: see Henon's paper for details.
def g(y, t):
    return -1.0/y

r2 = ode(g).set_integrator('dopri5')
r2.set_initial_value(r.t, r.y)

r2.integrate(0.5)

y_final = r2.t
t_final = r2.y

# Error: difference between found and analytical solution
print t_final - np.log(2)

【讨论】：