有关混沌控制的一个程序

shunfuli

关于Lorenz系统的控制，驱动系统方程是：
dx1/dt=10*(x2-x1),
dx2/dt=28*x1-x2-x1*x3,
dx3/dt=-8*x3/3+x1*x2;
响应系统方程是：
dy1/dt=10*(y2-y1)
dy2/dt=28*y1-y2-y1*y3
dy3/dt=-8*y3/3+y1*y2+k1*e.
其中dk1/dt=-10*e^2; e=y3-x3。
我用fortran编这个程序，但是结果溢出（出现NaN），下面是我的程序(fortran编写的，四阶Longe-Kutta积分)，不知道错在哪里，请高手指教，谢谢！

program Lorenz
implicit none
real*8:: o1,o2,o3,o4,p1,p2,p3,p4,q1,q2,q3,q4
real*8:: r1,r2,r3,r4,s1,s2,s3,s4,t1,t2,t3,t4
real*8:: u1,u2,u3,u4
real*8:: x1,x2,x3,y1,y2,y3
real*8:: e1,e2,e3
real*8:: h,t,k1,u
real*8:: f1,f2,f3,g1,g2,g3,v
integer:: i,n
open(9,file="adaptive.txt")
h=0.01d0
n=10000
x1=0.1
x2=-0.2
x3=0.3
y1=1.0
y2=2.0
y3=3.0
k1=-1.0
do i=1,n
     t=h*i
     o1=h*f1(x1,x2,x3)
     p1=h*f2(x1,x2,x3)
     q1=h*f3(x1,x2,x3)
     r1=h*g1(y1,y2,y3)                                    
     s1=h*g2(y1,y2,y3)
     t1=h*g3(y1,y2,y3,k1,x3)
     u1=h*v(x3,y3)  

     o2=h*f1(x1+0.5*o1,x2+0.5*p1,x3+0.5*q1)
     p2=h*f2(x1+0.5*o1,x2+0.5*p1,x3+0.5*q1)
     q2=h*f3(x1+0.5*o1,x2+0.5*p1,x3+0.5*q1)
     r2=h*g1(y1+0.5*r1,y2+0.5*s1,y3+0.5*t1)  
     s2=h*g2(y1+0.5*r1,y2+0.5*s1,y3+0.5*t1)
     t2=h*g3(y1+0.5*r1,y2+0.5*s1,y3+0.5*t1,k1+0.5*u1,x3+0.5*q1)
     u2=h*v(x3+0.5*q1,y3+0.5*t1)

     o3=h*f1(x1+0.5*o2,x2+0.5*p2,x3+0.5*q2)
     p3=h*f2(x1+0.5*o2,x2+0.5*p2,x3+0.5*q2)
     q3=h*f3(x1+0.5*o2,x2+0.5*p2,x3+0.5*q2)
     r3=h*g1(y1+0.5*r2,y2+0.5*s2,y3+0.5*t2)  
     s3=h*g2(y1+0.5*r2,y2+0.5*s2,y3+0.5*t2)
     t3=h*g3(y1+0.5*r2,y2+0.5*s2,y3+0.5*t2,k1+0.5*u2,x3+0.5*q2)
     u3=h*v(x3+0.5*q2,y3+0.5*t2)

     o4=h*f1(x1+o3,x2+p3,x3+q3)
     p4=h*f2(x1+o3,x2+p3,x3+q3)
     q4=h*f3(x1+o3,x2+p3,x3+q3)
     r4=h*g1(y1+r3,y2+s3,y3+t3)  
     s4=h*g2(y1+r3,y2+s3,y3+t3)
     t4=h*g3(y1+r3,y2+s3,y3+t3,k1+u3,x3+q3)
     u4=h*v(x3+q3,y3+t3)

     x1=x1+(o1+2.0*o2+2.0*o3+o4)/6.0
     x2=x2+(p1+2.0*p2+2.0*p3+p4)/6.0
     x3=x3+(q1+2.0*q2+2.0*q3+q4)/6.0
     y1=y1+(r1+2.0*r2+2.0*r3+r4)/6.0
     y2=y2+(s1+2.0*s2+2.0*s3+s4)/6.0
     y3=y3+(t1+2.0*t2+2.0*t3+t4)/6.0
     k1=k1+(u1+2.0*u2+2.0*u3+u4)/6.0
     write(*,*) y1,y2,y3
  enddo
end program


function f1(x1,x2,x3)
  real*8:: f1,x1,x2,x3
  f1=10.0*(x2-x1)
end

function f2(x1,x2,x3)
  real*8:: f2,x1,x2,x3
  f2=28.0*x1-x2-x1*x3
end
function f3(x1,x2,x3)
  real*8:: f3,x1,x2,x3
  f3=-8.0*x3/3.0+x1*x2
end

function g1(y1,y2,y3)
  real*8:: g1,y1,y2,y3
  g1=10.0*(y2-y1)
end

function g2(y1,y2,y3)
  real*8:: g2,y1,y2,y3
  g2=28.0*y1-y2-y1*y3
end

function g3(y1,y2,y3,k1,x3)
  real*8:: g3,y1,y2,y3,x3,k1
  g3=-8.0*y3/3.0+y1*y2+k1*(y3-x3)
end

function v(x3,y3)
  real*8:: v,x3,y3
  v=-10.0*(y3-x3)**2
end

shunfuli

自己顶一下，我的疑问是控制项参与积分吗？

cailiang

哥们，现在我正在学混沌控制，请问你的是什么控制方法？OGY?我用fortran90,你的错误好像是write(*,*),应该是write(9,*)吧，你open的是open(9,file....),而且调用函数的时候也出现几处错误,unused dummy，我的QQ：254875089（主论坛）.交流一下啊。我把程序程序修改了一下，能够运行了。结果也有了，我不太清楚是不是对的。一并附上，可以讨论下。

shunfuli

谢谢cailiang ，你将我的程序改了几处，我觉得起关键作用的是积分步长h由原来的0.01改为0.001，将你的程序中的h改为0.01还是会溢出的。为什么积分步长能影响程序呢？

cailiang

这个我也试了，我也不太清楚，步长大一点就会溢出，关于这个我也没有找到原因，我想一定还是程序里的某处错误了。那个画出来的图有意义吗？我已加你的QQ，有什么不懂的，还须向你请教。

mechanic05

控制项肯定要参与积分。对混沌系统来讲，加入的控制就是相当于外加一个驱动，而驱动要参与系统运动的。建议用MATLAB试试，简单好用。

shunfuli

哦，谢谢mechanic05.