使用 Runge-Kutta 4 求解罗斯勒吸引子答案

【问题标题】：Solving Rossler Attractor using Runge-Kutta 4使用 Runge-Kutta 4 求解罗斯勒吸引子
【发布时间】：2015-12-21 03:17:55
【问题描述】：

我正在尝试使用 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

【问题讨论】：

你犯的第一个错误是没有在缩进上花费足够的钱来使代码结构清晰，无论是对你自己还是对像我这样的随机陌生人。第二个错误可能与rossler的参数列表中的参数T有关；子程序根本不使用它。检查那个。然后摆脱无用的continue 语句——据我所知，这意味着所有这些语句。虽然 Fortran 对您在名称中使用的大小写不敏感（即 Yb 是与 yb 相同的变量），但我们人类发现对这些细节的注意力不集中会分散注意力。
除了上面的那些 cmets（虽然我不知道 t 因为我不知道方程式所以我不打算检查），我建议使用子例程的显式接口：编译器会比我更愉快地检查您的调用。
三十年来没有接触过fortran，我对subroutine rossler() 中的real […,] a,b,c 持怀疑态度——那不是声明局部变量a, b, c，让它们初始化在 0.0 时，可能隐藏读取的内容？
@greybeard 这肯定是答案——除了没有初始化（到 0.0 或其他任何东西），也没有隐藏任何东西。其他名为 a 等的变量不在同一范围内。
@greybeard 如果我没有在子例程 (rossler) 中包含 a、b、c 的声明，则会显示一条错误消息，要求指定它们的隐式类型。我该如何解决这个问题？

标签： algorithm fortran runge-kutta

【解决方案1】：

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

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 平面上的 2D 投影：

【讨论】：

是否有理由在rk4 中在deriv 的函数调用中使用矢量化操作但在两者之间循环？难道不能也表示为k1 = h*d等吗？
@LutzL 你是绝对正确的，但我故意留下显式循环，以便 OP 更容易进行比较（例如，不需要“CONTINUE”语句）。但是，是的，我希望 OP 也将它们重写为矢量化形式。
嗨，我从来不知道deriv。是标准的吗？（我检查了它是否适用于 gfortran 和 ifort），但看不到它是在哪里定义的
嗨@BaRud 这个deriv 是子程序rk4 的虚拟参数，所以它的名字实际上是任意的。它将绑定（或连接）到call rk4 ( rossler, t, n, y ) 的相应实际参数，以便deriv 在rk4 例程内等于rossler。没有名为deriv 的子例程（在其他文件或标准中）。

【解决方案2】：

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

明确并在问题上添加 cmets：a、b 和 c 代替该问题中的 omega；子程序rossler 作为函数fcn。

对该问题的回答说明了如何解决此问题。

【讨论】：