algorithm - 使用Runge-Kutta 4解决Rossler吸引子

我试图用RK-4，参数a=0.2，b=0.2，c=6，初始条件x0=-5.6，y0=0，z0=0，得到Rossler吸引子系统的解。我试着用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

最佳答案

尽管这个问题似乎是another question的重复，但在这里我附加了一个经过最少修改的代码，以便op可以将其与原始代码进行比较。基本的修改是我删除了所有未使用的变量，将a、b、c和h移动到参数模块，并清除了不必要的语句（如CONTINUE）。Fortran没有引入新的特性（包括interface的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

通过选择参数a = 0.2, b = 0.2, c = 5.7和nstep = 7000，修改后的代码给出了所谓的Rössler attractor，它非常漂亮，并且与Wiki页面中显示的模式非常接近因此，只要稍加修改，我相信OP也会得到类似的图片（看看模式如何根据参数变化可能会很有趣）。
轨迹在xy平面上的二维投影：

关于algorithm - 使用Runge-Kutta 4解决Rossler吸引子，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/34381444/