我正在運用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))