zhailiangjun 发表于 2007-12-24 20:09

在运用ode4求解微分方程組時出現NAN是什么問題啊

我正在運用ODE45求解一個微分方程組時,老是出現NAN的解,不知道那里出了問題,請各位高手指點一下子,小弟在這里先謝謝了。另外如果要畫出這個方程組的poincare截面,程序该怎么写啊?
function dy=zhzfun(t,y)
global w
global b
global n
global r
global s
global a
global a12
dy=zeros(4,1)
dy(1)=w*y(2)+s*y(4)
dy(2)=-(a+a12)*(n^2)*b*exp(-2*b*(y(1)-r))...
    +(a+a12)*(n^2)*b*(2-exp(-b*(y(1)-r)))*exp(-b*(y(1)-r))...
   -a12*(n^2)*(b*exp(-2*b*(y(1)-r))*(2-exp(-b*(y(3)-r)))*exp(-b*(y(3)-r))...
   -b*(2-exp(-b*(y(1)-r)))*exp(-b*(y(1)-r))*(2-exp(-b*(y(3)-r)))...
   *exp(-b*(y(3)-r)))...
   *((2-exp(-b*(y(1)-r)))*exp(-b*(y(1)-r))*(2-exp(-b*(y(3)-r)))...
   *exp(-b*(y(3)-r)))^(1/2)...
      -1/4*1.0571*(n^2)*(-2*b*exp(-b*(y(1)-r))+2*b*exp(-b*(y(1)-r)-b*(y(3)-r))...
      -(b*exp(-2*b*(y(1)-r))*(2-exp(-b*(y(3)-r)))*exp(-b*(y(3)-r))...
      -b*(2-exp(-b*(y(1)-r)))*exp(-b*(y(1)-r))*(2-exp(-b*(y(3)-r)))*exp(-b*(y(3)-r)))...
      *((2-exp(-b*(y(1)-r)))*exp(-b*(y(1)-r))*(2-exp(-b*(y(3)-r)))*exp(-b*(y(3)-r)))^(-1/2))
dy(3)=w*y(4)+s*y(2)
dy(4)=-(a+a12)*(n^2)*b*exp(-2*b*(y(3)-r))...
    +(a+a12)*(n^2)*b*(2-exp(-b*(y(3)-r)))*exp(-b*(y(3)-r))...
   -a12*(n^2)*(b*exp(-2*b*(y(3)-r))*(2-exp(-b*(y(1)-r)))*exp(-b*(y(1)-r))...
   -b*(2-exp(-b*(y(1)-r)))*exp(-b*(y(1)-r))*(2-exp(-b*(y(3)-r)))...
   *exp(-b*(y(3)-r)))...
   *((2-exp(-b*(y(1)-r)))*exp(-b*(y(1)-r))*(2-exp(-b*(y(3)-r)))...
   *exp(-b*(y(3)-r)))^(1/2)...
      -1/4*1.0571*(n^2)*(-2*b*exp(-b*(y(3)-r))+2*b*exp(-b*(y(3)-r)-b*(y(1)-r))...
      -(b*exp(-2*b*(y(3)-r))*(2-exp(-b*(y(1)-r)))*exp(-b*(y(1)-r))...
      -b*(2-exp(-b*(y(3)-r)))*exp(-b*(y(3)-r))*(2-exp(-b*(y(1)-r)))*exp(-b*(y(1)-r)))...
      *((2-exp(-b*(y(3)-r)))*exp(-b*(y(3)-r))*(2-exp(-b*(y(1)-r)))*exp(-b*(y(1)-r)))^(-1/2))

happy 发表于 2007-12-25 09:44

从你给的代码看,最大的可能是出现了0/0的情况
自己慢慢检查一下,什么地方有问题吧

这个只能靠自己一步步地检查

zhailiangjun 发表于 2007-12-25 10:52

回复 #2 happy 的帖子

谢谢了,我好好检查一下
页: [1]
查看完整版本: 在运用ode4求解微分方程組時出現NAN是什么問題啊