增加求解常微分方程的函数 odeint

MeteoInfo

参考scipy的odeint函数，在MeteoInfo 3.5.5版本中增加了求解常微分方程的函数 odeint 。调用了Apache Commons Math库中的相关功能，并用Jython进行了封装。

求解受重力摩擦作用的摆锤角θ的二阶微分方程示例：

1. from mipylib.numeric.integrate import odeint
   
2. 
   
3. def pend(y, t, b, c):
   
4.     theta, omega = y
   
5.     dydt = [omega, -b*omega - c*np.sin(theta)]
   
6.     return dydt
   
7. 
   
8. b = 0.25
   
9. c = 5.0
   
10. y0 = [np.pi - 0.1, 0.0]
    
11. t = np.linspace(0, 10, 101)
    
12. sol = odeint(pend, y0, t, args=(b, c))
    
13. 
    
14. plot(t, sol[:, 0], 'b', label='theta(t)', linewidth=2)
    
15. plot(t, sol[:, 1], 'g', label='omega(t)', linewidth=2)
    
16. legend(loc='lower right')
    
17. xlabel('t')
    
18. grid()

from mipylib.numeric.integrate import odeint<br /> <br /> def pend(y, t, b, c):<br />     theta, omega = y<br />     dydt = [omega, -b*omega - c*np.sin(theta)]<br />     return dydt<br /> <br /> b = 0.25<br /> c = 5.0<br /> y0 = [np.pi - 0.1, 0.0]<br /> t = np.linspace(0, 10, 101)<br /> sol = odeint(pend, y0, t, args=(b, c))<br /> <br /> plot(t, sol[:, 0], 'b', label='theta(t)', linewidth=2)<br /> plot(t, sol[:, 1], 'g', label='omega(t)', linewidth=2)<br /> legend(loc='lower right')<br /> xlabel('t')<br /> grid()

Lorenz吸引子常微分方程组求解：

1. from mipylib.numeric import integrate
   
2. 
   
3. def lorenz(p,t,s,r,b):
   
4.     x,y,z = p         
   
5.     return s*(y-x), x*(r-z)-y, x*y-b*z   # dx/dt,dy/dt,dz/dt
   
6. 
   
7. t = np.arange(0, 30, 0.01)
   
8. track1 = integrate.odeint(lorenz, (0.0,1.00,0.0), t, args=(10.0,28.0,2.6))
   
9. track2 = integrate.odeint(lorenz, (0.0,1.01,0.0), t, args=(10.0,28.0,2.6))
   
10. 
    
11. axes3d()
    
12. plot3(track1[:,0], track1[:,1], track1[:,2], linewidth=2, color='r')            
    
13. plot3(track2[:,0], track2[:,1], track2[:,2], linewidth=2, color='g')   

from mipylib.numeric import integrate<br /> <br /> def lorenz(p,t,s,r,b):<br />     x,y,z = p          <br />     return s*(y-x), x*(r-z)-y, x*y-b*z   # dx/dt,dy/dt,dz/dt<br /> <br /> t = np.arange(0, 30, 0.01)<br /> track1 = integrate.odeint(lorenz, (0.0,1.00,0.0), t, args=(10.0,28.0,2.6))<br /> track2 = integrate.odeint(lorenz, (0.0,1.01,0.0), t, args=(10.0,28.0,2.6))<br /> <br /> axes3d()<br /> plot3(track1[:,0], track1[:,1], track1[:,2], linewidth=2, color='r')             <br /> plot3(track2[:,0], track2[:,1], track2[:,2], linewidth=2, color='g')   

贫道敬孔

王老师太强了