用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問題

回答

相關問題