0

我正在尝试使用 RK-4 获得罗斯勒吸引子系统的解决方案，参数 a=0.2，b=0.2，c=6 和初始条件 x0=-5.6，y0=0，z0=0。我尝试使用 Fortran 求解，但结果即使在 1000 次迭代后也只显示初始条件。我犯了什么错误？

implicit none
external rossler
integer::i,j=0,n,nstep
real::a,b,c,y1(3),t0,dt,t1,t2,ya(3),yb(3),yd(3),t,x0,y0,z0,x(1000),y(1000),z(1000),k1(3),k2(3),k3(3),k4(3),h
print *, "enter the values of a,b,c"
read (*,*) a,b,c
print *, "enter the values of x0,y0,z0"
read (*,*) x0,y0,z0
n=3
t0=0.0
h=0.05
ya(1)=x0
ya(2)=y0
ya(3)=z0
nstep=1000
do i=1,nstep
t1=t0
t2=t0+h
call rk4(rossler,t1,t2,1,N,k1,k2,k3,k4,Ya,Y1,Yb)

x(i)=ya(1)
y(i)=ya(2)
z(i)=ya(3)
open (99,file="rossler.txt")
write(99,*) x(i),y(i),z(i)
end do

end program

subroutine rossler(T,Yd,YB,N)
implicit none
integer n
real t,yb(n),yd(n),a,b,c
yd(1)=-yb(2)-yb(3)
Yd(2)=yb(1)+a*yb(2)
Yd(3)=b+yb(3)*(yb(1)-c)
return
end

subroutine rk4(rossler,t1,t2,nstep,N,k1,k2,k3,k4,Ya,Y1,Yb)
implicit none
external rossler
integer nstep,n,i,j
REAL T1,T2,Ya(N),k1(n),k2(n),k3(n),k4(n),H,Y1(N),T,yb(n)
T=T1+(I-1)*H
CALL rossler(T,Yb,Ya,N)
DO J=1,N
k1(j)=YB(J)*H
end do
CONTINUE
CALL rossler(T+0.5*H,Yb,Ya+k1*0.5,N)
DO J=1,N
k2(j)=YB(J)*H
enddo
CONTINUE
CALL rossler(T+0.5*H,Yb,Ya+k2*0.5,N)
DO J=1,N
K3(J)=YB(J)*H
enddo
CONTINUE
CALL rossler(T+H,Yb,Ya+k3,N)
DO J=1,n
K4(J)=YB(J)*H
Y1(J)=Ya(J)+(k1(j)+k4(j)+2.0*(k2(j)+k3(j)))/6.0
enddo
CONTINUE
DO J=1,N
Ya(J)=Y1(j)
enddo
CONTINUE
enddo
RETURN
END
4

2 回答 2

2

尽管该问题似乎与另一个问题重复，但在这里我附上了一个经过最少修改的代码，以便 OP 可以将其与原始代码进行比较。基本的修改是我删除了所有未使用的变量，将a、b、c和h移至参数模块，并清理了不必要的语句（如CONTINUE）。没有引入 Fortran 的新功能（包括interface块 for rossler），因此希望可以直接看到代码是如何更改的。

module params
    real :: a, b, c, h
end module

program main
    use params, only: a, b, c, h
    implicit none
    external rossler
    integer :: i, n, nstep
    real :: t, y(3)

    a = 0.2
    b = 0.2
    c = 5.7
    n = 3
    t = 0.0
    h = 0.05
    y(1) = -5.6
    y(2) = 0.0
    y(3) = 0.0
    nstep = 7000

    open(99, file="rossler.txt")
    do i = 1,nstep
        call rk4 ( rossler, t, n, y )
        write(99,*) y(1), y(2), y(3)
    end do

end program

subroutine rossler ( t, dy, y, n )
    use params, only: a, b, c
    implicit none
    integer n
    real t, dy(n), y(n)
    dy(1) = -y(2) - y(3)
    dy(2) = y(1) + a * y(2)
    dy(3) = b + ( y(1) - c ) * y(3)
end

subroutine rk4 ( deriv, t, n, y )
    use params, only: h
    implicit none
    external deriv
    integer n, j
    real y(n), t, k1(n), k2(n), k3(n), k4(n), d(n)

    call deriv ( t, d, y, n )
    do j = 1,n
        k1(j) = d(j) * h
    enddo

    call deriv ( t+0.5*h, d, y+k1*0.5, n )
    DO j = 1,n
        k2(j) = d(j) * h
    enddo

    call deriv ( t+0.5*h, d, y+k2*0.5, n )
    do j = 1,n
        k3(j) = d(j) * h
    enddo

    call deriv ( t+h, d, y+k3, n )
    do j = 1,n
        k4(j) = d(j) * h
        y(j) = y(j) + ( k1(j) + k4(j) + 2.0 * (k2(j) + k3(j)) ) / 6.0
    enddo

    t = t + h
end

通过选择参数 asa = 0.2, b = 0.2, c = 5.7和nstep = 7000，修改后的代码给出了所谓的Rössler 吸引子，它非常漂亮，并且看起来与 Wiki 页面中显示的模式非常接近。因此，通过最小的修改，我相信 OP 也会得到类似的图片（看看模式如何根据参数变化可能会很有趣）。

轨迹在 xy 平面上的 2D 投影： 投影到 xy 平面

于 2015-12-20T18:45:51.553 回答
0

这里的问题与另一个问题完全相同，尽管我不能再投票关闭作为重复项。

明确并添加对问题的评论：a，b并从该问题c中代替omega；子程序rossler作为函数fcn。

该问题的答案解决了如何解决此问题。

于 2015-12-20T14:22:49.473 回答