用 python 求解 ODE答案

【问题标题】：Solving ODE with python用 python 求解 ODE
【发布时间】：2015-12-01 06:30:25
【问题描述】：

我尝试使用 scipy.integrate.ode 求解一阶 ode 方程。

dh/dx = - s0 * (1 - (hn/h)^3)/(1 - (hc/h)^3)

初始条件是 x = 0 和 h = 10。
当 x = 15000 时，s0 = 0.0005。

hn = (f*q^2/8*g*s0)^(1/3)
hc = (q^2/g)^(1/3)

f、q、g 是常数。

我使用的方法是node bdf，但是得到的结果和matlab解决的答案不一样。答案应该是这样的： https://dl.dropboxusercontent.com/u/18438495/result.png

谁能看出问题所在？

import numpy as np
from scipy.integrate import ode
import matplotlib.pyplot as plt

def waterdepth(t, y):
    if t < 15000:
        s0 = 0.001
    elif t >= 15000:
        s0 = 0.0005
    q = 3.72
    f = 0.03
    g = 9.81
    hn = (f*q*q/8*g*s0)**(1.0/3.0)
    hc = (q*q/g)**(1.0/3.0)
    return -s0 * (1.0 - (hn/y)**3)/(1.0 - (hc/y)**3)

y0 = 10.0
t0 = 0.0

solver = ode(waterdepth).set_integrator('node', method = 'bdf') 
solver.set_initial_value(y0, t0)
dt = 100.0
t1 = 25000

x = []
h = []
while solver.successful() and solver.t < t1:
    x.append(solver.t)
    solver.integrate(solver.t + dt)
    h.append(solver.y)

    plt.plot (x, h)

【问题讨论】：

标签： python ode

【解决方案1】：

return -s0 * (1.0 - (hn/y)**3)/(1.0 - (hc/y)**3)

与上面的等式不同。

dh/dx = - s0 * (1 - (hn-h)^3)/(1 - (hc-h)^3)

【讨论】：

感谢您的通知。我的 python 代码中的等式是正确的。