|
楼主 |
发表于 2009-3-25 14:16
|
显示全部楼层
这是分岔图的程序:
real function f(v,w)
if (v.eq.0.0) then
f=-1.0
else if (v.lt.0.0) then
f=1.0 /(1.0+v*w)
else if (v.gt.0.0) then
f=-1.0/(1.0+v*w)
end if
end function f
real a1,a2,a3,a4,b1,b2,b3,b4,c1,c2,c3,c4,d1,d2,d3,d4
real x1,x10,x2,x20,y10,y20,t0,t,h,w,u1e,u2e,u
g1(t,x1,x2)=x2-x1-(x1-0.1*t)
g2(t,x1,x2)=x1-x2-(x2-0.1*t)
open(1,FILE='1分岔图.TXT',STATUS='old')
open(2,FILE='2分岔图.TXT',STATUS='old')
x10=0.00367950
y10=0.0
x20=0.00613586
y20=0.0
h=0.001
do w=0.0,1.5,0.005
do t0=0.0,500.0,0.001
a1=h*y10
b1=h*g1(t0,x10,x20)+h*f(y10,w)
c1=h*y20
d1=h*g2(t0,x10,x20)+h*f(y20,w)
a2=h*(y10+1.0*b1/2)
b2=h*g1(t0+1.0*h/2,x10+1.0*a1/2,x20+1.0*c1/2)+h*f(y10+1.0*b1/2,w)
c2=h*(y20+1.0*d1/2)
d2=h*g2(t0+1.0*h/2,x10+1.0*a1/2,x20+1.0*c1/2)+h*f(y20+1.0*d1/2,w)
a3=h*(y10+1.0*b2/2)
b3=h*g1(t0+1.0*h/2,x10+1.0*a2/2,x20+1.0*c2/2)+h*f(y10+1.0*b2/2,w)
c3=h*(y20+1.0*d2/2)
d3=h*g2(t0+1.0*h/2,x10+1.0*a2/2,x20+1.0*c2/2)+h*f(y20+1.0*d2/2,w)
a4=h*(y10+b3)
b4=h*g1(t0+h,x10+a3,x20+c3)+h*f(y10+b3,w)
c4=h*(y20+d3)
d4=h*g2(t0+h,x10+a3,x20+c3)+h*f(y20+d3,w)
x10=x10+1.0*(a1+2.0*a2+2.0*a3+a4)/6
y10=y10+1.0*(b1+2.0*b2+2.0*b3+b4)/6
x20=x20+1.0*(c1+2.0*c2+2.0*c3+c4)/6
y20=y20+1.0*(d1+2.0*d2+2.0*d3+d4)/6
u1e=0.1*t0-1.0/(1.0+0.1*w)
u2e=0.1*t0-1.0/(1.0+0.1*w)
u=1.0*(x10-u1e+x20-u2e)/2.0
if (u==0.0.and.t0>50.0) then
write(1,*)w,' ',y10
end if
!if (x20==u2e.and.t0>3000.0) then
!write(2,*)w,' ',y20
!end if
end do
end do
end |
|