Solving Rossler Attractor using Runge-Kutta 4

Question

I am trying to get a solution for the Rossler attractor system using RK-4, with parameters a=0.2, b=0.2, c=6 and initial conditions x0=-5.6, y0=0, z0=0. I tried solving using Fortran but the result is only displaying the initial conditions even after 1000 iterations. What mistakes am I making?

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

Answer 1

Although the question seems a duplicate of another question, here I am attaching a minimally modified code so that the OP can compare it with the original one. The essential modifications are that I have removed all the unused variables, moved a, b, c, and h to a parameter module, and cleaned up unnecessary statements (like CONTINUE). No newer features of Fortran introduced (including interface block for rossler), so it is hopefully straight-forward to see how the code has been changed.

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

By choosing the parameters as a = 0.2, b = 0.2, c = 5.7 and nstep = 7000, the modified code gave the so-called Rössler attractor, which is very beautiful and appears close in pattern to that displayed in the Wiki page. So with the minimal modifications, I believe the OP will also get a similar picture (it may be interesting to see how the pattern changes depending on parameters).

2D projection of the trajectory onto the xy plane:

Answer 2

The problem here is exactly the same as in another question, although I can no longer vote to close as a duplicate.

To make explicit and add the comments on the question: a, b and c take the place of omega from that question; the subroutine rossler as the function fcn.

An answer to that question addresses how this issue can be resolved.