2015-09-05 166 views
2

我試着用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) 

回答

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

是從在頂部的方程不同。

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

謝謝您的通知。我的Python代碼中的公式是正確的。 –