用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)
答
return -s0 * (1.0 - (hn/y)**3)/(1.0 - (hc/y)**3)
是从在顶部的方程不同。
dh/dx = - s0 * (1 - (hn-h)^3)/(1 - (hc-h)^3)
谢谢您的通知。我的Python代码中的公式是正确的。 –