有关混沌控制的一个程序
关于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 自己顶一下，我的疑问是控制项参与积分吗？ 哥们，现在我正在学混沌控制，请问你的是什么控制方法？OGY?我用fortran90,你的错误好像是write(*,*),应该是write(9,*)吧，你open的是open(9,file....),而且调用函数的时候也出现几处错误,unused dummy，我的QQ：254875089（主论坛）.交流一下啊。我把程序程序修改了一下，能够运行了。结果也有了，我不太清楚是不是对的。一并附上，可以讨论下。 谢谢cailiang ，你将我的程序改了几处，我觉得起关键作用的是积分步长h由原来的0.01改为0.001，将你的程序中的h改为0.01还是会溢出的。为什么积分步长能影响程序呢？ 这个我也试了，我也不太清楚，步长大一点就会溢出，关于这个我也没有找到原因，我想一定还是程序里的某处错误了。那个画出来的图有意义吗？我已加你的QQ，有什么不懂的，还须向你请教。
回复 沙发 shunfuli 的帖子
控制项肯定要参与积分。对混沌系统来讲，加入的控制就是相当于外加一个驱动，而驱动要参与系统运动的。建议用MATLAB试试，简单好用。 哦，谢谢mechanic05.
页:
[1]