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

Question

我正在使用Scipy 14.0来解决一个普通微分方程系统，这个系统描述了一个气泡在静止液体中因浮力而垂直上升的动态情况。具体来说，我有一个方程来表示气泡的上升速度U与气泡半径R之间的关系，也就是U=dz/dt=f(R)，还有一个方程表示半径的变化与R和U之间的关系，也就是dR/dT=f(R,U)。代码中出现的其他内容都是材料属性。

我想实现一个功能来考虑z的物理限制，显然这个限制是由液体高度H决定的。因此，我实现了一种z<=H的约束，以便在需要时提前停止积分：我使用了set_solout来做到这一点。现在的情况是，代码运行正常并给出了不错的结果，但set_solout根本没有起作用（看起来z_constraint根本没有被调用……）。你知道这是为什么吗？

有没有人有更聪明的想法，或许可以在z=H时准确中断（也就是一个最终值问题）？这是正确的方法/工具吗，还是我应该重新表述这个问题？

提前谢谢你

Emi

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)

Answer 1

你遇到了这个问题。要让set_solout正常工作，你需要在调用set_integrator之后，set_initial_value之前立刻调用它。如果你在代码中做这个修改（并为g设置一个值），当z >= H时，积分就会结束，这正是你想要的结果。

为了找到气泡到达水面时的确切时间，你可以在solout结束积分后，改变变量，重新以z为变量进行积分，直到z = H。有一篇论文描述了这种技术，链接是M. Henon, Physica 5D, 412 (1982)；你也可以参考这个讨论，可能会对你有帮助。下面是一个非常简单的例子，展示如何找到时间t使得y(t) = 0.5，给定dy/dt = -y：

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)