我已经使用四阶龙格-库塔微分建立了一个钟摆模拟,其中一切都是一步一步完成的:
from scipy import *
from matplotlib.pyplot import *
##A pendulum simulation using fourth order
##Runge-Kutta differentiation
ts=.05 #time step size
td=20 #trial duration
te=int(td/ts) #no of timesteps
mu=0.1 #friction factor
m=1 #mass
g=9.81 #grav. acceleration
l=1 #length
th=[((rand()*2)-1)*pi] #initial angle
om=[0] #initial angular velocity
u=0 #torque
for j in range(te):
#Euler approximation
th.append(th[j] + ts*om[j])
f1 = (-mu*om[j] + m*g*l*sin(th[j]) + u)/(m*(l^2))
om.append(om[j] + ts*f1)
#approximation 1 at mid-interval
th2 = th[j+1] + (ts/2)*om[j+1]
f2 = (-mu*om[j+1] + m*g*l*sin(th[j+1]) + u)/(m*(l^2))
om2 = om[j+1] + (ts/2)*f2
#approximation 2 at mid-interval
th3 = th2 + (ts/2)*om2
f3 = (-mu*om2 + m*g*l*sin(th2) + u)/( m*(l^2))
om3 = om2 + (ts/2)*f3
#approximation at next time step
th4 = th3 + (ts)*om3
f4 = (-mu*om3 + m*g*l*sin(th3) + u)/( m*(l^2))
om4 = om3 + (ts)*f4
dth=(om[j] + 2*om[j+1] + 2*om2 + om3)/6
dom=(f1 + 2*f2 + 2*f3 + f4)/6
th[j+1] = th[j] + ts*dth
om[j+1] = om[j] + ts*dom
subplot(211),plot(th),xlabel('Angle'),ylabel('')
subplot(212),plot(om,'r'),xlabel('Angular velocity'),ylabel('')
show()
这很好用,但我不觉得优雅,这就是为什么我决定使用 odeint() 函数重写它:
from scipy import integrate
from scipy import *
from matplotlib.pyplot import *
##A pendulum simulation
th=((rand()*2)-1)*pi #initial angle
om=0 #initial angular velocity
u=0 #torque
y0 = [om, th] #initial values
t = linspace(0, 20, 400) #
def f(y, t):
mu=0.1 #friction factor
m=1 #mass
g=9.81 #grav. acceleration
l=1 #length
return (y[0],
(-mu*y[0] + m*g*l*sin(y[1]) + u)/(m*(l^2)))
r = integrate.odeint(f, y0, t)
subplot(211),plot(t,r[:,1]),xlabel('Angle'),ylabel('')
subplot(212),plot(t,r[:,0],'r'),xlabel('Angular velocity'),ylabel('')
show()
而这个新版本看起来一点也不像钟摆。谁能告诉我我做错了什么?