用python解决ODE问题

我试着用scipy.integrate.ode解决一个1阶的ode方程。用python解决ODE问题

DH/DX = - S 0 *（1 - （HN /小时）^ 3）/（1 - （HC /小时）^ 3）

的初始条件为x = 0和h = 10。
S0 = 0.001当x < 15000， S0 = 0.0005当x> = 15000。

HN =（F * q^2/8 * G * S0）^（1/3）
HC =（（1/3）

f，q，g是恒定的。

我使用的方法是节点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)

来源

2015-09-05 John Yang

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

是从在顶部的方程不同。

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

来源

2015-09-05 03:40:39

谢谢您的通知。我的Python代码中的公式是正确的。 –

用python解决ODE问题

回答

相关问题