I've tried to solve a 1-st order ode equation using scipy.integrate.ode.
dh/dx = - s0 * (1 - (hn/h)^3)/(1 - (hc/h)^3)
The initial condition is x = 0 and h = 10.
s0 = 0.001 when x < 15000,
s0 = 0.0005 when x >= 15000.
hn = (f*q^2/8*g*s0)^(1/3)
hc = (q^2/g)^(1/3)
f, q, g are constant.
The method I used is node bdf, but the result I got is different than the answer solved by matlab. The answer should be like this: https://dl.dropboxusercontent.com/u/18438495/result.png
Can anyone see the problem?
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)
