我需要使用python代码来解决这个问题,该代码基于重力受空气阻力影响的假设y"=g-(1/70)*(y')^2,g=9.8,y(0)=0,y'(0)=0
我已经用Euler Heun方法解决了类似的问题y"=y-(1/10)*y',y(0)=0,y'(0)=1,[0,无穷大]
我该如何修改我必须解决的第一个问题,即空气阻力下降?有什么建议吗?
h= 0.0001
y = 0.0
t = 0.0
v = 1.0
while y >= 0:
ym = v*h+y
vm = v - h*(y +v*0.1)
y = y + 0.5*h*(v+vm)
v = v - 0.5*h*(y + ym + 0.1*(v + vm))
t += h
error = abs(1.0 - y/ym)
if error > 1.0e-8:
h*=0.1
if error < 1.0e-10:
h*=1.1
print y, v, t
注意:我也尝试了一个简单的欧拉方法,但的误差不够准确
delta_t = .000001 #time step size
y = 0. #initial height
g = 9.8 #gravitational acceleration
t = 0. #initial time
y_prime = 0. #initial velocity
y_2_prime = g #initial acceleration
while y < 100:
y = y_prime * delta_t + y
y_prime = y_2_prime * delta_t + y_prime
y_2_prime = g - 0.014285714* y_prime**2
t = t + delta_t
print t
您可以使用龙格-库塔方法,它基本上是一种更高阶的欧拉方法。你可以查看维基百科了解更多关于这方面的详细信息。它真的很常见,所以你在谷歌上搜索它并找到一些关于如何实现它的好例子应该不会有任何困难
还有其他方法,但大多数人(至少是那些开始ODE集成的人)都喜欢Runge-Kutta。