MeteoInfo的个人博客分享 http://blog.sciencenet.cn/u/MeteoInfo

博文

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

已有 1151 次阅读 2023-2-13 10:12 |个人分类:MeteoInfo|系统分类:科研笔记

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

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

from mipylib.numeric.integrate import odeint

def pend(y, t, b, c):
   theta, omega = y
   dydt = [omega, -b*omega - c*np.sin(theta)]
   return dydt

b = 0.25
c = 5.0
y0 = [np.pi - 0.1, 0.0]
t = np.linspace(0, 10, 101)
sol = odeint(pend, y0, t, args=(b, c))

plot(t, sol[:, 0], 'b', label='theta(t)', linewidth=2)
plot(t, sol[:, 1], 'g', label='omega(t)', linewidth=2)
legend(loc='lower right')
xlabel('t')
grid()

odeint_pendulum.png

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

from mipylib.numeric import integrate

def lorenz(p,t,s,r,b):
   x,y,z = p          
   return s*(y-x), x*(r-z)-y, x*y-b*z   # dx/dt,dy/dt,dz/dt

t = np.arange(0, 30, 0.01)
track1 = integrate.odeint(lorenz, (0.0,1.00,0.0), t, args=(10.0,28.0,2.6))
track2 = integrate.odeint(lorenz, (0.0,1.01,0.0), t, args=(10.0,28.0,2.6))

axes3d()
plot3(track1[:,0], track1[:,1], track1[:,2], linewidth=2, color='r')            
plot3(track2[:,0], track2[:,1], track2[:,2], linewidth=2, color='g')  

odeint_lorentz.png



https://blog.sciencenet.cn/blog-611158-1376037.html

上一篇:发布排放源清单处理系统EMIPS
下一篇:FY4A AGRI L1真彩色图合成
收藏 IP: 106.120.73.*| 热度|

0

该博文允许注册用户评论 请点击登录 评论 (0 个评论)

数据加载中...

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2024-4-16 22:49

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部