请选择 进入手机版 | 继续访问电脑版
爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 1614|回复: 1

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

[复制链接]

新浪微博达人勋

发表于 2023-2-13 10:15:35 | 显示全部楼层 |阅读模式

登录后查看更多精彩内容~

您需要 登录 才可以下载或查看,没有帐号?立即注册 新浪微博登陆

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

求解受重力摩擦作用的摆锤角θ的二阶微分方程示例:
  1. from mipylib.numeric.integrate import odeint

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

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

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

odeint_pendulum.png

Lorenz吸引子常微分方程组求解:
  1. from mipylib.numeric import integrate

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

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

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

odeint_lorentz.png
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2023-2-13 12:11:43 | 显示全部楼层
王老师太强了
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表