N_BOM

New Burst-Oscillation Mode in Paced One-Dimensional Excitable Systems

阅读

• 下载地址：https://chinazhang-my.sharepoint.com/:b:/g/personal/rep_rebirth_zh_zhangwenhao_icu/EeskR_6Y7XRItPh5hr9ZhYABzvQ5x8RsJeja8GYQgsm54g?e=OSbvS7

复现

figure 1

figure 1_A

• 链

  四阶龙格库塔法

• 四阶龙格库塔遍历步长下的链，返回链，循环，计算量太大,除首节点外有问题还需要修改 。

module BOM
    implicit none
    integer,parameter :: N=2,L=500,MaxT=10000
    integer,allocatable :: chain_matrix(:,:)
    real,parameter :: h=0.02,PI=3.1415926
contains

    subroutine fnf(x,fx,n,i,j,t)
        implicit none
        real :: x(n,L),fx(n,L)
        real :: epsilon=0.04,a=0.84,b=0.07,D=1.0,A1=2.5,f=4.5,effect
        integer :: N,i,j,t,c
        effect=0.0
        do c=1,L,1
            effect=effect+chain_matrix(i,c)*(x(1,j)-x(1,i))
        end do
        if(i==1) then
            fx(1,i)=1.0/epsilon*(x(1,i)*(1-x(1,i))*(x(1,i)-(x(2,i)+b)/a))+A1*sin(2*PI*f*t*h)
        else
            fx(1,i)=1.0/epsilon*(x(1,i)*(1-x(1,i))*(x(1,i)-(x(2,i)+b)/a))+D*effect
        end if
        if(x(1,i)<1.0/3.0) then
            fx(2,i)=-x(2,i)
        elseif(x(1,i)>=1.0/3.0.and.x(1,i)<=1.0) then
            fx(2,i)=1.0-6.75*x(1,i)*(x(1,i)-1)**2-x(2,i)
        else
            fx(2,i)=1-x(2,i)
        end if
        return
    end subroutine fnf

    subroutine rk4(N,h,x,i,j,t)
        implicit none
        integer :: N,i,j,t
        real :: h,xx(N,L),x(N,L),fx(N,L)
        real :: kx1(N,L),kx2(N,L),kx3(N,L),kx4(N,L)
        xx=x
        call fnf(xx,fx,N,i,j,t)
        kx1=h*fx
        xx=x+0.5*kx1
        call fnf(xx,fx,N,i,j,t)
        kx2=h*fx
        xx=x+0.5*kx2
        call fnf(xx,fx,N,i,j,t)
        kx3=h*fx
        xx=x+kx3
        call fnf(xx,fx,N,i,j,t)
        kx4=h*fx
        x=x+(kx1+2.0*kx2+2.0*kx3+kx4)/6.0
        return
    end subroutine rk4

    !链
    subroutine chain(node_number)
        implicit none
        integer :: node_number,a
        allocate(chain_matrix(node_number,node_number))
        !初始化网络矩阵
        chain_matrix=0
        !赋值节点关系
        do a=2,node_number,1
            chain_matrix(a,a-1)=1
        end do
    end subroutine chain

end module BOM


program main
    use BOM
    implicit none
    integer :: i,j,t,b
    real :: x(N,L)
    x=0.0
    call chain(L)
    open(20,file="data_1.txt")
    open(40,file="data_50.txt")
    do t=1,MaxT,1
        do i=1,L,1
            if(i==1) then
                    call rk4(N,h,x(:,i),i,1,t)
                    write(20,*) t*h,x(1,i)
                end if
            do j=1,L,1
                if(chain_matrix(i,j)==1) then
                    call rk4(N,h,x(:,i),i,j,t)
                    if(i==50) then
                        write(40,*) t*h,x(1,i)
                    end if
                end if
            end do
        end do
    end do
    deallocate(chain_matrix)
end program main

• 四阶龙格库塔遍历步长下的节点，返回节点，循环。

module BOM
    implicit none
    integer,parameter :: N=2,L=500,MaxT=10000
    integer,allocatable :: chain_matrix(:,:)
    real,parameter :: h=0.02,PI=3.1415926
    real :: preNode(L),x(N,L)
contains

    subroutine fnf(x,fx,n,i,j,t)
        implicit none
        real :: x(N),fx(N)
        real :: epsilon=0.04,a=0.84,b=0.07,D=1.0,A1=2.5,f=4.5,effect
        integer :: N,i,j,t,c
        effect=0.0
        do c=1,L,1
            effect=effect+chain_matrix(i,c)*(preNode(j)-x(1))
        end do
        if(i==1) then
            fx(1)=1.0/epsilon*(x(1)*(1-x(1))*(x(1)-(x(2)+b)/a))+A1*sin(2*PI*f*t*h)
        else
            fx(1)=1.0/epsilon*(x(1)*(1-x(1))*(x(1)-(x(2)+b)/a))+D*effect
        end if
        if(x(1)<1.0/3.0) then
            fx(2)=-x(2)
        elseif(x(1)>=1.0/3.0.and.x(1)<=1.0) then
            fx(2)=1.0-6.75*x(1)*(x(1)-1)**2-x(2)
        else
            fx(2)=1-x(2)
        end if
        return
    end subroutine fnf

    subroutine rk4(N,h,x,i,j,t)
        implicit none
        integer :: N,i,j,t,i1
        real :: h,xx(N),fx(N),x(N)
        real :: kx1(N),kx2(N),kx3(N),kx4(N)
        xx=x
        call fnf(xx,fx,N,i,j,t)
        kx1=h*fx
        xx=xx+0.5*kx1
        call fnf(xx,fx,N,i,j,t)
        kx2=h*fx
        xx=xx+0.5*kx2
        call fnf(xx,fx,N,i,j,t)
        kx3=h*fx
        xx=xx+kx3
        call fnf(xx,fx,N,i,j,t)
        kx4=h*fx
        x=x+(kx1+2.0*kx2+2.0*kx3+kx4)/6.0
        return
    end subroutine rk4

    !链
    subroutine chain(node_number)
        implicit none
        integer :: node_number,a
        allocate(chain_matrix(node_number,node_number))
        !初始化网络矩阵
        chain_matrix=0
        !赋值节点关系
        do a=2,node_number,1
            chain_matrix(a,a-1)=1
        end do
    end subroutine chain
end module BOM

program main
    use BOM
    implicit none
    integer :: i,j,t,b
    x=0.0
    call chain(L)
    open(20,file="data_1.txt")
    open(40,file="data_50.txt")
    do t=1,MaxT,1
        preNode=0.0
        do i=1,L,1
            if(i==1) then
                call rk4(N,h,x(:,i),i,1,t)
                preNode(i)=x(1,i)
                write(20,*) t*h,x(1,i)
            end if
            do j=1,L,1
                if(chain_matrix(i,j)==1) then
                    call rk4(N,h,x(:,i),i,j,t)
                    preNode(i)=x(1,i)
                    if(i==50) then
                       write(40,*) t*h,x(1,i)
                    end if
                end if
            end do
        end do
    end do
    deallocate(chain_matrix)
end program main

一阶欧拉法

• $f(x1)$=$f(x0)$+h*$f(x0)$

module BOM
    implicit none
    integer,parameter :: L=500,MaxT=10000
    integer,allocatable :: chain_matrix(:,:)
    real,parameter :: h=0.02,PI=3.1415926
contains

    subroutine fnf()
        implicit none
        real :: epsilon=0.04,a=0.84,b=0.07,D=1.0,A1=2.5,f=4.5,effect,u(L),v(L),fun_u(L),fun_v(L)
        integer :: i,j,t
        u=0.0
        v=0.0
        do t=1,MaxT,1
            do i=1,L,1
                effect=0
                do j=1,L
                    effect=effect+chain_matrix(i,j)*(u(j)-u(i))
                end do
                if(i==1) then
                    fun_u(i)=u(i)+h*(1.0/epsilon*u(i)*(1.0-u(i))*(u(i)-(v(i)+b)/a)+A1*sin(2.0*PI*f*t*h))
                else
                    fun_u(i)=u(i)+h*(1.0/epsilon*u(i)*(1.0-u(i))*(u(i)-(v(i)+b)/a)+D*effect)
                end if
                if(u(i)<1.0/3.0)then
                    fun_v(i)=v(i)+h*(-v(i))
                else if(u(i)>=1.0/3.0.and.u(i)<=1)then
                    fun_v(i)=v(i)+h*(1-6.75D0*u(i)*(u(i)-1)**2-v(i))
                else if(u(i)>1.0)then
                    fun_v(i)=v(i)+h*(1-v(i))
                end if
                if(i==1) then
                    write(20,*)t*h,fun_u(i)
                end if
                if(i==50) then
                    write(40,*)t*h,fun_v(i)
                end if
            end do
            u=fun_u
            v=fun_v
        end do
        close(20)
        close(40)
        return
    end subroutine fnf
    !链
    subroutine chain(node_number)
        implicit none
        integer :: node_number,a
        allocate(chain_matrix(node_number,node_number))
        !初始化网络矩阵
        chain_matrix=0
        !赋值节点关系
        do a=2,node_number,1
            chain_matrix(a,a-1)=1
        end do
    end subroutine chain

end module BOM


program main
    use BOM
    implicit none
    open(20,file="data_1.txt")
    open(40,file="data_50.txt")
    call chain(L)
    call fnf
    deallocate(chain_matrix)
end program main

• $(A,f)$=$(2.5,4.5)$

figure 1_B

• $(A,f)$=$(1.0,4.5)$

figure 1_C

• 环

四阶龙格库塔法

module BOM
    implicit none
    integer,parameter :: N=3,L=500,MaxT=10000
    real,parameter :: h=0.02,PI=3.1415926
    integer,allocatable :: ring_matrix(:,:)
    real :: preNode(L)
contains

    subroutine fnf(x,fx,N,i,j)
        implicit none
        real :: x(N),fx(N)
        real :: epsilon=0.04,a=0.84,b=0.07,D=1.0,A1=2.5,f=4.5,delta,ans
        integer :: N,i,j
        integer :: c
        if(i==1) then
            delta=1.0
        else
            delta=0.0
        end if
        ans=0.0
        do c=1,L,1
            ans=ans+ring_matrix(i,c)*(preNode(j)-x(1))
        end do
        fx(1)=1.0/epsilon*(x(1)*(1.0-x(1))*(x(1)-(x(2)+b)/a))+delta*A1*sin(2.0*PI*f*x(3))+D*ans
        if(x(1)<1.0/3.0) then
            fx(2)=-x(2)
        elseif(x(1)>=1.0/3.0.and.x(1)<=1.0) then
            fx(2)=1.0-6.75*x(1)*(x(1)-1.0)**2-x(2)
        else
            fx(2)=1.0-x(2)
        end if
        fx(3)=1.0
        return
    end subroutine fnf

    subroutine rk4(N,x,i,j)
        implicit none
        integer :: i,j,N
        real :: xx(N),x(N),fx(N)
        real :: kx1(N),kx2(N),kx3(N),kx4(N)
        xx=x
        call fnf(xx,fx,N,i,j)
        kx1=h*fx
        xx=x+0.5*kx1
        call fnf(xx,fx,N,i,j)
        kx2=h*fx
        xx=x+0.5*kx2
        call fnf(xx,fx,N,i,j)
        kx3=h*fx
        xx=x+kx3
        call fnf(xx,fx,N,i,j)
        kx4=h*fx
        x=x+(kx1+2.0*kx2+2.0*kx3+kx4)/6.0
        return
    end subroutine rk4

    !链
    subroutine ring(node_number)
        implicit none
        integer :: node_number,a
        allocate(ring_matrix(node_number,node_number))
        !初始化网络矩阵
        ring_matrix=0
        !赋值节点关系
        do a=2,node_number,1
            ring_matrix(a,a-1)=1
        end do
        !首位相连
        ring_matrix(1,node_number)=1
    end subroutine ring

end module BOM


program main
    use BOM
    implicit none
    integer :: i,j,t
    real :: x(N,L)
    call ring(L)
    open(20,file="data.txt")
    preNode=0.0
    x=0.0
    do t=1,MaxT,1
        do i=1,L,1
            do j=1,L,1
                if(ring_matrix(i,j)==1) then
                    call rk4(N,x(:,i),i,j)
                    preNode(i)=x(1,i)
                    if(i==1) then
                        write(20,*) t*h,x(1,i),x(3,i)
                    end if
                end if
            end do
        end do
        preNode(L)=0.0
    end do
    deallocate(ring_matrix)
end program main

一阶欧拉法

module BOM
    implicit none
    integer,parameter :: L=500,MaxT=10000
    integer,allocatable :: ring_matrix(:,:)
    real,parameter :: h=0.02,PI=3.1415926
contains

    subroutine fnf()
        implicit none
        real :: epsilon=0.04,a=0.84,b=0.07,D=1.0,A1=2.5,f=4.5,effect,u(L),v(L),fun_u(L),fun_v(L),delta
        integer :: i,j,t
        u=0.0
        v=0.0
        do t=1,MaxT,1
            do i=1,L,1
                effect=0.0
                do j=1,L,1
                    effect=effect+ring_matrix(i,j)*(u(j)-u(i))
                end do
                if(i==1) then
                    delta=1.0
                else
                    delta=0.0
                end if
                fun_u(i)=u(i)+h*(1.0/epsilon*(u(i)*(1.0-u(i))*(u(i)-(v(i)+b)/a))+delta*A1*sin(2*PI*f*t*h)+D*effect)
                if(u(i)<1.0/3.0) then
                    fun_v(i)=v(i)+h*(-v(i))
                else if(u(i)>=1.0/3.0.and.u(i)<=1) then
                    fun_v(i)=v(i)+h*(1-6.75*u(i)*(u(i)-1)**2-v(i))
                else if(u(i)>1.0) then
                    fun_v(i)=v(i)+h*(1-v(i))
                end if
                if(i==1) then
                    write(20,*)t*h,fun_u(i)
                end if
                if(i==500) then
                    write(40,*)t*h,fun_v(i)
                end if
            end do
            fun_u(L)=0.0
            fun_v(L)=0.0
            u=fun_u
            v=fun_v
        end do
        close(20)
        close(40)
        return
    end subroutine fnf
    !环
    subroutine ring(node_number)
        implicit none
        integer :: node_number,a
        allocate(ring_matrix(node_number,node_number))
        !初始化网络矩阵
        ring_matrix=0
        !赋值节点关系
        do a=2,node_number,1
            ring_matrix(a,a-1)=1
        end do
        !首尾相连
        ring_matrix(1,node_number)=1
    end subroutine ring

end module BOM


program main
    use BOM
    implicit none
    open(20,file="data_1.txt")
    open(40,file="data_50.txt")
    call ring(L)
    call fnf
    deallocate(ring_matrix)
end program main

• $(A,f)$=$(2.5,4.5)$

figure 1_D

四阶龙格库塔法

module BOM
    implicit none
    integer,parameter :: N=3,L=500,MaxT=10000
    real,parameter :: h=0.02,PI=3.1415926
    integer,allocatable :: ring_matrix(:,:)
    real :: preNode(L)
contains

    subroutine fnf(x,fx,N,i,j)
        implicit none
        real :: x(N),fx(N)
        real :: epsilon=0.04,a=0.84,b=0.07,D=1.0,A1=2.5,f=4.5,delta,ans
        integer :: N,i,j
        integer :: c
        if(i==1) then
            delta=1.0
        else
            delta=0.0
        end if
        ans=0.0
        do c=1,L,1
            ans=ans+ring_matrix(i,c)*(preNode(j)-x(1))
        end do
        fx(1)=1.0/epsilon*(x(1)*(1.0-x(1))*(x(1)-(x(2)+b)/a))+delta*A1*sin(2.0*PI*f*x(3))+D*ans
        if(x(1)<1.0/3.0) then
            fx(2)=-x(2)
        elseif(x(1)>=1.0/3.0.and.x(1)<=1.0) then
            fx(2)=1.0-6.75*x(1)*(x(1)-1.0)**2-x(2)
        else
            fx(2)=1.0-x(2)
        end if
        fx(3)=1.0
        return
    end subroutine fnf

    subroutine rk4(N,x,i,j)
        implicit none
        integer :: i,j,N
        real :: xx(N),x(N),fx(N)
        real :: kx1(N),kx2(N),kx3(N),kx4(N)
        xx=x
        call fnf(xx,fx,N,i,j)
        kx1=h*fx
        xx=x+0.5*kx1
        call fnf(xx,fx,N,i,j)
        kx2=h*fx
        xx=x+0.5*kx2
        call fnf(xx,fx,N,i,j)
        kx3=h*fx
        xx=x+kx3
        call fnf(xx,fx,N,i,j)
        kx4=h*fx
        x=x+(kx1+2.0*kx2+2.0*kx3+kx4)/6.0
        return
    end subroutine rk4

    !环
    subroutine ring(node_number)
        implicit none
        integer :: node_number,a
        allocate(ring_matrix(node_number,node_number))
        !初始化网络矩阵
        ring_matrix=0
        !赋值节点关系
        do a=2,node_number,1
            ring_matrix(a,a-1)=1
        end do
        !首位相连
        ring_matrix(1,node_number)=1
    end subroutine ring

end module BOM


program main
    use BOM
    implicit none
    integer :: i,j,t
    real :: x(N,L)
    call ring(L)
    open(20,file="data.txt")
    preNode=0.0
    x=0.0
    do t=1,MaxT,1
        do i=1,L,1
            do j=1,L,1
                if(ring_matrix(i,j)==1) then
                    call rk4(N,x(:,i),i,j)
                    preNode(i)=x(1,i)
                    if(i==1) then
                        write(20,*) t*h,x(1,i),x(3,i)
                    end if
                end if
            end do
        end do
    end do
    deallocate(ring_matrix)
end program main

• 思路，既然每个节点遍历相同的时间，可以不使用步长下遍历链或环，直接对每个 节点遍历相应时间，但这种方法可以很快看到节点趋势，时间数据庞大时可以使用，与实际走向相符，但这种方法适用较少传递节点低损失传播。

module BOM
    implicit none
    integer,parameter :: n=3,L=500
    integer,allocatable :: ring_matrix(:,:)
    integer :: i,j
    real :: preNode(10000)
    contains
    !环
    subroutine ring(L)
        implicit none
        integer :: L,i
        allocate(ring_matrix(L,L))
        !初始化网络矩阵
        ring_matrix=0
        !赋值节点关系
        do i=1,L-1,1
            ring_matrix(i,i+1)=1
        end do !首位相连
        ring_matrix(L,1)=1
        deallocate(ring_matrix)
    end subroutine ring
   subroutine fnf(x,fx,n,i,j)
       implicit none
       real :: x(n),fx(n),epsilon=0.04,a=0.84,b=0.07,D=1.0,A1=2.5,f=4.5,delta
       integer :: n,i,j
       real,parameter :: PI=3.1415926
       if(mod(i,500)==1) then
           delta=1.0
       else
           delta=0.0
       end if
       fx(1)=1.0/epsilon*(x(1)*(1-x(1))*(x(1)-(x(2)+b)/a))+delta*A1*sin(2*PI*f*x(3))+D*(preNode(j)-x(1))
       if(x(1)<1.0/3.0) then
           fx(2)=-x(2)
       elseif(x(1)>=1.0/3.0.and.x(1)<=1.0) then
           fx(2)=1.0-6.75*x(1)*(x(1)-1)**2-x(2)
       else
           fx(2)=1-x(2)
       end if
       fx(3)=1
       return
   end subroutine fnf

    subroutine rk4(n,h,x,i,j)
        implicit none
        integer :: n,i,j
        real :: h,xx(n),x(n),fx(n)
        real :: kx1(n),kx2(n),kx3(n),kx4(n)
        xx=x
        call fnf(xx,fx,n,i,j)
        kx1=h*fx
        xx=x+0.5*kx1
        call fnf(xx,fx,n,i,j)
        kx2=h*fx
        xx=x+0.5*kx2
        call fnf(xx,fx,n,i,j)
        kx3=h*fx
        xx=x+kx3
        call fnf(xx,fx,n,i,j)
        kx4=h*fx
        x=x+(kx1+2.0*kx2+2.0*kx3+kx4)/6.0
        return
    end subroutine rk4

end module BOM
program main
    use BOM
    implicit none
    real :: x(n),h=0.02
    call ring(L)
    x(1)=0.0
    x(2)=0.0
    preNode=0.0
    open(20,file="data.txt")
    do i=1,L+1,1
        x(3)=0.0
        do j=1,10000,1
            call rk4(n,h,x,i,j)
            preNode(j)=x(1)
            if(i==501) then
                write(20,*) x
            end if
        end do
    end do
    return
end program main

一阶欧拉法

module BOM
    implicit none
    integer,parameter :: L=500,MaxT=10000
    integer,allocatable :: ring_matrix(:,:)
    real,parameter :: h=0.02,PI=3.1415926
contains

    subroutine fnf()
        implicit none
        real :: epsilon=0.04,a=0.84,b=0.07,D=1.0,A1=2.5,f=4.5,effect,u(L),v(L),fun_u(L),fun_v(L),delta
        integer :: i,j,t
        u=0.0
        v=0.0
        do t=1,MaxT,1
            do i=1,L,1
                effect=0.0
                do j=1,L,1
                    effect=effect+ring_matrix(i,j)*(u(j)-u(i))
                end do
                if(i==1) then
                    delta=1.0
                else
                    delta=0.0
                end if
                fun_u(i)=u(i)+h*(1.0/epsilon*(u(i)*(1.0-u(i))*(u(i)-(v(i)+b)/a))+delta*A1*sin(2*PI*f*t*h)+D*effect)
                if(u(i)<1.0/3.0) then
                    fun_v(i)=v(i)+h*(-v(i))
                else if(u(i)>=1.0/3.0.and.u(i)<=1) then
                    fun_v(i)=v(i)+h*(1-6.75*u(i)*(u(i)-1)**2-v(i))
                else if(u(i)>1.0) then
                    fun_v(i)=v(i)+h*(1-v(i))
                end if
                if(i==1) then
                    write(20,*)t*h,fun_u(i)
                end if
                if(i==500) then
                    write(40,*)t*h,fun_v(i)
                end if
            end do
            u=fun_u
            v=fun_v
        end do
        close(20)
        close(40)
        return
    end subroutine fnf
    !环
    subroutine ring(node_number)
        implicit none
        integer :: node_number,a
        allocate(ring_matrix(node_number,node_number))
        !初始化网络矩阵
        ring_matrix=0
        !赋值节点关系
        do a=2,node_number,1
            ring_matrix(a,a-1)=1
        end do
        !首尾相连
        ring_matrix(1,node_number)=1
    end subroutine ring

end module BOM


program main
    use BOM
    implicit none
    open(20,file="data_1.txt")
    open(40,file="data_50.txt")
    call ring(L)
    call fnf
    deallocate(ring_matrix)
end program main

figure 2

• 查看大致趋势

module BOM
    implicit none
    integer,parameter :: n=3,L=500
    integer :: i,j
    real :: preNode(35000)
    contains

    subroutine fnf(x,fx,n,i,j)
        implicit none
        real :: x(n),fx(n),epsilon=0.04,a=0.84,b=0.07,D=1.0,A1=2.5,f=4.5,delta
        integer :: n,i,j
        real,parameter :: PI=3.1415926
        if(i==1.and.x(3)<=50.0) then
            fx(1)=1.0/epsilon*(x(1)*(1-x(1))*(x(1)-(x(2)+b)/a))+A1*sin(2*PI*f*x(3))
        else
            if(mod(i,500)==1) then
                delta=1.0
            else
                delta=0.0
            end if
            fx(1)=1.0/epsilon*(x(1)*(1-x(1))*(x(1)-(x(2)+b)/a))+delta*A1*sin(2*PI*f*x(3))+D*(preNode(j)-x(1))
        end if
        if(x(1)<1.0/3.0) then
            fx(2)=-x(2)
        elseif(x(1)>=1.0/3.0.and.x(1)<=1.0) then
            fx(2)=1.0-6.75*x(1)*(x(1)-1)**2-x(2)
        else
            fx(2)=1-x(2)
        end if
        fx(3)=1
        return
    end subroutine fnf

    subroutine rk4(n,h,x,i,j)
        implicit none
        integer :: n,i,j
        real :: h,xx(n),x(n),fx(n)
        real :: kx1(n),kx2(n),kx3(n),kx4(n)
        xx=x
        call fnf(xx,fx,n,i,j)
        kx1=h*fx
        xx=x+0.5*kx1
        call fnf(xx,fx,n,i,j)
        kx2=h*fx
        xx=x+0.5*kx2
        call fnf(xx,fx,n,i,j)
        kx3=h*fx
        xx=x+kx3
        call fnf(xx,fx,n,i,j)
        kx4=h*fx
        x=x+(kx1+2.0*kx2+2.0*kx3+kx4)/6.0
        return
    end subroutine rk4

end module BOM


program main
    use BOM
    implicit none
    real :: x(n),h=0.02
    x(1)=0.0
    x(2)=0.0
    preNode=0.0
    open(20,file="data.txt")
    do i=1,3*L+1,1
        x(3)=0.0
        do j=1,35000,1
            call rk4(n,h,x,i,j)
            preNode(j)=x(1)
            if(i==1) then
                write(20,*) x
            end if
        end do
    end do
end program main

• 图片左轴有误，应为500节点

figure 2_A

四阶龙格库塔法

module BOM
    implicit none
    integer,parameter :: N=3,L=500,MaxT=10000
    real,parameter :: h=0.02,PI=3.1415926
    integer,allocatable :: ring_matrix(:,:)
    real :: preNode(L)
contains

    subroutine fnf(x,fx,N,i,j)
        implicit none
        real :: x(N),fx(N)
        real :: epsilon=0.04,a=0.84,b=0.07,D=1.0,A1=2.5,f=4.5,delta,ans
        integer :: N,i,j
        integer :: c
        if(i==1) then
            delta=1.0
        else
            delta=0.0
        end if
        ans=0.0
        do c=1,L,1
            ans=ans+ring_matrix(i,c)*(preNode(j)-x(1))
        end do
        if(x(3)<=50.and.i==1) then
            fx(1)=1.0/epsilon*(x(1)*(1.0-x(1))*(x(1)-(x(2)+b)/a))+A1*sin(2.0*PI*f*x(3))
        else
            fx(1)=1.0/epsilon*(x(1)*(1.0-x(1))*(x(1)-(x(2)+b)/a))+delta*A1*sin(2.0*PI*f*x(3))+D*ans
        end if
        if(x(1)<1.0/3.0) then
            fx(2)=-x(2)
        elseif(x(1)>=1.0/3.0.and.x(1)<=1.0) then
            fx(2)=1.0-6.75*x(1)*(x(1)-1.0)**2-x(2)
        else
            fx(2)=1.0-x(2)
        end if
        fx(3)=1.0
        return
    end subroutine fnf

    subroutine rk4(N,x,i,j)
        implicit none
        integer :: i,j,N
        real :: xx(N),x(N),fx(N)
        real :: kx1(N),kx2(N),kx3(N),kx4(N)
        xx=x
        call fnf(xx,fx,N,i,j)
        kx1=h*fx
        xx=x+0.5*kx1
        call fnf(xx,fx,N,i,j)
        kx2=h*fx
        xx=x+0.5*kx2
        call fnf(xx,fx,N,i,j)
        kx3=h*fx
        xx=x+kx3
        call fnf(xx,fx,N,i,j)
        kx4=h*fx
        x=x+(kx1+2.0*kx2+2.0*kx3+kx4)/6.0
        return
    end subroutine rk4

    !链
    subroutine ring(node_number)
        implicit none
        integer :: node_number,a
        allocate(ring_matrix(node_number,node_number))
        !初始化网络矩阵
        ring_matrix=0
        !赋值节点关系
        do a=2,node_number,1
            ring_matrix(a,a-1)=1
        end do
        !首位相连
        ring_matrix(1,node_number)=1
    end subroutine ring

end module BOM


program main
    use BOM
    implicit none
    integer :: i,j,t
    real :: x(N,L)
    call ring(L)
    open(20,file="data_1.txt")
    preNode=0.0
    x=0.0
    do t=1,MaxT,1
        do i=1,L,1
            do j=1,L,1
                if(ring_matrix(i,j)==1) then
                    call rk4(N,x(:,i),i,j)
                    preNode(i)=x(1,i)
                    if(i==1) then
                        write(20,*) t*h,x(1,i)
                    end if
                end if
            end do
        end do
    end do
    deallocate(ring_matrix)
end program main

一阶欧拉法

module BOM
    implicit none
    integer,parameter :: L=500,MaxT=35000
    integer,allocatable :: ring_matrix(:,:)
    real,parameter :: h=0.02,PI=3.1415926
contains

    subroutine fnf()
        implicit none
        real :: epsilon=0.04,a=0.84,b=0.07,D=1.0,A1=2.5,f=4.5,effect,u(L),v(L),fun_u(L),fun_v(L),delta,stepT
        integer :: i,j,t
        u=0.0
        v=0.0
        do t=1,MaxT,1
            do i=1,L,1
                stepT=h*t
                effect=0.0
                do j=1,L,1
                    effect=effect+ring_matrix(i,j)*(u(j)-u(i))
                end do
                if(i==1) then
                    delta=1.0
                else
                    delta=0.0
                end if
                if(i==1.and.t<=2500) then
                    fun_u(i)=u(i)+h*(1.0/epsilon*(u(i)*(1.0-u(i))*(u(i)-(v(i)+b)/a))+A1*sin(2.0*PI*f*stepT))
                else
                    fun_u(i)=u(i)+h*(1.0/epsilon*(u(i)*(1.0-u(i))*(u(i)-(v(i)+b)/a))+delta*A1*sin(2.0*PI*f*stepT)+D*effect)
                end if
                if(u(i)<1.0/3.0) then
                    fun_v(i)=v(i)+h*(-v(i))
                else if(u(i)>=1.0/3.0.and.u(i)<=1.0) then
                    fun_v(i)=v(i)+h*(1.0-6.75*u(i)*(u(i)-1.0)**2-v(i))
                else if(u(i)>1.0) then
                    fun_v(i)=v(i)+h*(1.0-v(i))
                end if
                if(i==1) then
                    write(20,*)t*h,fun_u(i)
                end if
                if(i==500) then
                    write(40,*)t*h,fun_u(i)
                end if
                if(mod(t,100)==0)then
                    !write(20,*) i，t*h,u(i)
                end if  
            end do
            u=fun_u
            v=fun_v
        end do
        close(20)
        close(40)
        return
    end subroutine fnf
    !环
    subroutine ring(node_number)
        implicit none
        integer :: node_number,a
        allocate(ring_matrix(node_number,node_number))
        !初始化网络矩阵
        ring_matrix=0
        !赋值节点关系
        do a=2,node_number,1
            ring_matrix(a,a-1)=1
        end do
        !首尾相连
        ring_matrix(1,node_number)=1
    end subroutine ring

end module BOM


program main
    use BOM
    implicit none
    open(20,file="data_1.txt")
    open(40,file="data_500.txt")
    call ring(L)
    call fnf
    deallocate(ring_matrix)
end program main

figure 2_B

module BOM
    implicit none
    integer,parameter :: L=500,MaxT=35000
    integer,allocatable :: ring_matrix(:,:)
    real,parameter :: h=0.02,PI=3.1415926
contains

    subroutine fnf()
        implicit none
        real :: epsilon=0.04,a=0.84,b=0.07,D=1.0,A1=2.5,f=4.5,effect,u(L),v(L),fun_u(L),fun_v(L),delta,stepT
        integer :: i,j,t
        u=0.0
        v=0.0
        do t=1,MaxT,1
            do i=1,L,1
                stepT=h*t
                effect=0.0
                do j=1,L,1
                    effect=effect+ring_matrix(i,j)*(u(j)-u(i))
                end do
                if(i==1) then
                    delta=1.0
                else
                    delta=0.0
                end if
                if(i==1.and.t<=2500) then
                    fun_u(i)=u(i)+h*(1.0/epsilon*(u(i)*(1.0-u(i))*(u(i)-(v(i)+b)/a))+A1*sin(2.0*PI*f*stepT))
                else
                    fun_u(i)=u(i)+h*(1.0/epsilon*(u(i)*(1.0-u(i))*(u(i)-(v(i)+b)/a))+delta*A1*sin(2.0*PI*f*stepT)+D*effect)
                end if
                if(u(i)<1.0/3.0) then
                    fun_v(i)=v(i)+h*(-v(i))
                else if(u(i)>=1.0/3.0.and.u(i)<=1.0) then
                    fun_v(i)=v(i)+h*(1.0-6.75*u(i)*(u(i)-1.0)**2-v(i))
                else if(u(i)>1.0) then
                    fun_v(i)=v(i)+h*(1.0-v(i))
                end if
                if(mod(t,100)==0)then
                    write(20,*) i,t*h,u(i)
                end if  
            end do
            u=fun_u
            v=fun_v
        end do
        close(20)
        close(40)
        return
    end subroutine fnf
    !环
    subroutine ring(node_number)
        implicit none
        integer :: node_number,a
        allocate(ring_matrix(node_number,node_number))
        !初始化网络矩阵
        ring_matrix=0
        !赋值节点关系
        do a=2,node_number,1
            ring_matrix(a,a-1)=1
        end do
        !首尾相连
        ring_matrix(1,node_number)=1
    end subroutine ring

end module BOM


program main
    use BOM
    implicit none
    open(20,file="data_1.txt")
    open(40,file="data_500.txt")
    call ring(L)
    call fnf
    deallocate(ring_matrix)
end program main

figure 2_C

figure 2_D

figure 2_E

figure 2_F

figure 2_G

figure 2_H

figure 2_I

figure 5

figure 5_A

module BOM
    implicit none
    integer,parameter :: L=500,MaxT=8750 !MaxT一个周期时间，T_Length
    integer,allocatable :: ring_matrix(:,:)
    real,parameter :: h=0.02,PI=3.1415926
contains

    subroutine fnf(T_Link)
        implicit none
        real :: epsilon=0.04,a=0.84,b=0.07,D=1.0,A1=2.5,f=4.5,effect,u(L),v(L),fun_u(L),fun_v(L),delta,stepT,T_Link,T_Interval=0.0,&
                U_change(1:MaxT)
        integer :: i,j,t,N_pulse=0,count_number=0
        u=0.0
        v=0.0
        do t=1,MaxT,1
            do i=1,L,1
                stepT=h*t
                effect=0.0
                do j=1,L,1
                    effect=effect+ring_matrix(i,j)*(u(j)-u(i))
                end do
                if(i==1) then
                    delta=1.0
                else
                    delta=0.0
                end if
                if(i==1.and.stepT<=T_Link) then
                    fun_u(i)=u(i)+h*(1.0/epsilon*(u(i)*(1.0-u(i))*(u(i)-(v(i)+b)/a))+A1*sin(2.0*PI*f*stepT))
                else
                    fun_u(i)=u(i)+h*(1.0/epsilon*(u(i)*(1.0-u(i))*(u(i)-(v(i)+b)/a))+delta*A1*sin(2.0*PI*f*stepT)+D*effect)
                end if
                if(u(i)<1.0/3.0) then
                    fun_v(i)=v(i)+h*(-v(i))
                else if(u(i)>=1.0/3.0.and.u(i)<=1.0) then
                    fun_v(i)=v(i)+h*(1.0-6.75*u(i)*(u(i)-1.0)**2-v(i))
                else if(u(i)>1.0) then
                    fun_v(i)=v(i)+h*(1.0-v(i))
                end if
                if(i==1) then
                    U_change(t)=fun_u(i)
                    if(U_change(t)>0.5) then
                        count_number=count_number+1
                    else if(count_number>30.and.U_change(t)<0.2) then
                        N_pulse=N_pulse+1
                        count_number=0
                    end if
                end if
            end do
            u=fun_u
            v=fun_v
        end do
        write(20,*) T_Link,N_pulse
        N_pulse=0
        return
    end subroutine fnf
    !环
    subroutine ring(node_number)
        implicit none
        integer :: node_number,a
        allocate(ring_matrix(node_number,node_number))
        !初始化网络矩阵
        ring_matrix=0
        !赋值节点关系
        do a=2,node_number,1
            ring_matrix(a,a-1)=1
        end do
        !首尾相连
        ring_matrix(1,node_number)=1
    end subroutine ring

end module BOM


program main
    use BOM
    implicit none
    real :: T_Link=0.0
    integer :: i
    open(20,file="T_LinkN_Pulse.txt")
    open(40,file="T_LinkT_Interval.txt")
    call ring(L)
    do i=1,8751,1 !T_Link_max=175s
        T_Link=i*h-h
        call fnf(T_Link)
    end do
    close(20)
    close(40)
    deallocate(ring_matrix)
end program main

figure 5_B

module BOM
    implicit none
    integer,parameter :: L=500,MaxT=9000 !MaxT一个周期时间，T_Length
    integer,allocatable :: ring_matrix(:,:)
    real,parameter :: h=0.02,PI=3.1415926
contains

    subroutine fnf(T_Link)
        implicit none
        real :: epsilon=0.04,a=0.84,b=0.07,D=1.0,A1=2.5,f=4.5,effect,u(L),v(L),fun_u(L),fun_v(L),delta,stepT,T_Link,T_Interval=0.0,&
                U_change(1:MaxT)
        integer :: i,j,t,N_pulse=0,count_number=0
        u=0.0
        v=0.0
        do t=1,MaxT,1
            do i=1,L,1
                stepT=h*t
                effect=0.0
                do j=1,L,1
                    effect=effect+ring_matrix(i,j)*(u(j)-u(i))
                end do
                if(i==1) then
                    delta=1.0
                else
                    delta=0.0
                end if
                if(i==1.and.stepT<=T_Link) then
                    fun_u(i)=u(i)+h*(1.0/epsilon*(u(i)*(1.0-u(i))*(u(i)-(v(i)+b)/a))+A1*sin(2.0*PI*f*stepT))
                else
                    fun_u(i)=u(i)+h*(1.0/epsilon*(u(i)*(1.0-u(i))*(u(i)-(v(i)+b)/a))+delta*A1*sin(2.0*PI*f*stepT)+D*effect)
                end if
                if(u(i)<1.0/3.0) then
                    fun_v(i)=v(i)+h*(-v(i))
                else if(u(i)>=1.0/3.0.and.u(i)<=1.0) then
                    fun_v(i)=v(i)+h*(1.0-6.75*u(i)*(u(i)-1.0)**2-v(i))
                else if(u(i)>1.0) then
                    fun_v(i)=v(i)+h*(1.0-v(i))
                end if
                if(i==2) then
                    U_change(t)=fun_u(i)
                    if(U_change(t)<0.5) then
                        count_number=count_number+1
                    else if(count_number>270.and.U_change(t)>0.5) then
                        T_Interval=count_number*h
                        count_number=0
                    else
                        count_number=0
                    end if
                end if
            end do
            u=fun_u
            v=fun_v
        end do
        write(*,*) T_Link,T_Interval
        write(40,*) T_Link,T_Interval
        N_pulse=0
        return
    end subroutine fnf
    !环
    subroutine ring(node_number)
        implicit none
        integer :: node_number,a
        allocate(ring_matrix(node_number,node_number))
        !初始化网络矩阵
        ring_matrix=0
        !赋值节点关系
        do a=2,node_number,1
            ring_matrix(a,a-1)=1
        end do
        !首尾相连
        ring_matrix(1,node_number)=1
    end subroutine ring

end module BOM


program main
    use BOM
    implicit none
    real :: T_Link=0.0
    integer :: i
    open(20,file="T_LinkN_Pulse.txt")
    open(40,file="T_LinkT_Interval.txt")
    call ring(L)
    do i=1,8751,1 !T_Link_max=175s
        T_Link=i*h-h
        call fnf(T_Link)
    end do
    close(20)
    close(40)
    deallocate(ring_matrix)
end program main

figure 6

• 这里A图的环不加驱动，给定初始值，传递振荡，一下代码是加了驱动，还没修改。

figure 6_A

module BOM
    implicit none
    integer,parameter :: MaxT=15000
    integer,allocatable :: ring_matrix(:,:)
    real,parameter :: h=0.02,PI=3.1415926
contains

    subroutine fnf(L)
        implicit none
        real :: epsilon=0.04,a=0.84,b=0.07,D=1.0,A1=2.5,f=4.5,effect,u(L),v(L),fun_u(L),fun_v(L),delta,stepT,&
                U_change(1:MaxT),T_chain=0.0,T_peak(2),T_length=0.0
        integer :: i,j,t,L,N_pulse=0,count_number=0
        u=0.0
        v=0.0
        do t=1,MaxT,1
            do i=1,L,1
                stepT=h*t
                effect=0.0
                do j=1,L,1
                    effect=effect+ring_matrix(i,j)*(u(j)-u(i))
                end do
                if(i==1) then
                    delta=1.0
                else
                    delta=0.0
                end if
                if(i==1.and.stepT<=10) then
                    fun_u(i)=u(i)+h*(1.0/epsilon*(u(i)*(1.0-u(i))*(u(i)-(v(i)+b)/a))+A1*sin(2.0*PI*f*stepT))
                else
                    fun_u(i)=u(i)+h*(1.0/epsilon*(u(i)*(1.0-u(i))*(u(i)-(v(i)+b)/a))+delta*A1*sin(2.0*PI*f*stepT)+D*effect)
                end if
                if(u(i)<1.0/3.0) then
                    fun_v(i)=v(i)+h*(-v(i))
                else if(u(i)>=1.0/3.0.and.u(i)<=1.0) then
                    fun_v(i)=v(i)+h*(1.0-6.75*u(i)*(u(i)-1.0)**2-v(i))
                else if(u(i)>1.0) then
                    fun_v(i)=v(i)+h*(1.0-v(i))
                end if
                if(i==1) then
                    U_change(t)=fun_u(i)
                    if(U_change(t)>0.5) then
                        count_number=count_number+1
                        if(N_pulse==2) then
                        T_length=stepT
                        end if
                    else if(count_number>30.and.U_change(t)<0.2) then
                        N_pulse=N_pulse+1
                        count_number=0
                        if(N_pulse==1) then
                            T_peak(1)=stepT
                        end if
                        if(N_pulse==2) then
                            T_peak(2)=stepT
                        end if
                    end if
                end if
            end do
            u=fun_u
            v=fun_v
        end do
        T_chain=T_peak(2)-T_peak(1)
        write(20,*) L,T_chain,T_length
        N_pulse=0
        return
    end subroutine fnf
    !环
    subroutine ring(node_number)
        implicit none
        integer :: node_number,a
        allocate(ring_matrix(node_number,node_number))
        !初始化网络矩阵
        ring_matrix=0
        !赋值节点关系
        do a=2,node_number,1
            ring_matrix(a,a-1)=1
        end do
        !首尾相连
        ring_matrix(1,node_number)=1
    end subroutine ring

end module BOM


program main
    use BOM
    implicit none
    integer :: i,L
    open(20,file="L_Tchain_Tlength.txt")
    do i=300,700,1  !间隔应该为10
        call ring(i)
        call fnf(i)
        deallocate(ring_matrix)
    end do
    close(20)
end program main

figure 6_B

module BOM
    implicit none
    integer,allocatable :: ring_matrix(:,:)
    real,parameter :: h=0.02,PI=3.1415926
contains

    subroutine fnf(L,T_link,MaxT)
        implicit none
        integer :: i,j,t,L,N_pulse=0,count_number=0,MaxT
        real :: epsilon=0.04,a=0.84,b=0.07,D=1.0,A1=2.5,f=4.5,effect,u(L),v(L),fun_u(L),fun_v(L),delta,stepT,&
                U_change(1:MaxT),T_chain=0.0,T_peak(2),T_length=0.0,T_link
        u=0.0
        v=0.0
        do t=1,MaxT,1
            do i=1,L,1
                stepT=h*t
                effect=0.0
                do j=1,L,1
                    effect=effect+ring_matrix(i,j)*(u(j)-u(i))
                end do
                if(i==1) then
                    delta=1.0
                else
                    delta=0.0
                end if
                if(i==1.and.stepT<=T_link) then
                    fun_u(i)=u(i)+h*(1.0/epsilon*(u(i)*(1.0-u(i))*(u(i)-(v(i)+b)/a))+A1*sin(2.0*PI*f*stepT))
                else
                    fun_u(i)=u(i)+h*(1.0/epsilon*(u(i)*(1.0-u(i))*(u(i)-(v(i)+b)/a))+delta*A1*sin(2.0*PI*f*stepT)+D*effect)
                end if
                if(u(i)<1.0/3.0) then
                    fun_v(i)=v(i)+h*(-v(i))
                else if(u(i)>=1.0/3.0.and.u(i)<=1.0) then
                    fun_v(i)=v(i)+h*(1.0-6.75*u(i)*(u(i)-1.0)**2-v(i))
                else if(u(i)>1.0) then
                    fun_v(i)=v(i)+h*(1.0-v(i))
                end if
                if(i==1) then
                    U_change(t)=fun_u(i)
                    if(U_change(t)>0.5) then
                        count_number=count_number+1
                    else if(count_number>=97.and.U_change(t)<0.2) then
                        N_pulse=N_pulse+1
                        count_number=0
                    end if
                end if
            end do
            u=fun_u
            v=fun_v
        end do
        write(20,*) L,N_pulse,T_link,MaxT
        N_pulse=0
        return
    end subroutine fnf
    !环
    subroutine ring(node_number)
        implicit none
        integer :: node_number,a
        allocate(ring_matrix(node_number,node_number))
        !初始化网络矩阵
        ring_matrix=0
        !赋值节点关系
        do a=2,node_number,1
            ring_matrix(a,a-1)=1
        end do
        !首尾相连
        ring_matrix(1,node_number)=1
    end subroutine ring

end module BOM


program main
    use BOM
    implicit none
    integer :: i,j=1,L
    real(kind=4) :: T_Link_Array(41)=(/106.879997,110.399994,113.899994,117.399994,120.860001,124.360001,127.879997,131.300003,&
            134.819992,138.339996,141.779999,145.279999,148.800003,152.259995,155.73999,159.259995,162.759995,166.199997,&
            169.699997,173.220001,176.660004,180.159988,183.679993,187.139999,190.619995,194.139999,197.639999,201.080002,&
            204.599991,208.099991,211.539993,215.059998,218.559998,222.019989,225.5,229.019989,232.519989,235.959991,&
            239.479996,243.0,246.419998/)
    open(20,file="L_N_pulse_max.txt")
    do i=300,700,10
        call ring(i)
        call fnf(i,T_Link_Array(j),int(T_Link_Array(j)/0.02))
        deallocate(ring_matrix)
        j=j+1
    end do
    close(20)
end program main

figure 6_C

module BOM
    implicit none
    integer,allocatable :: ring_matrix(:,:)
    real,parameter :: h=0.02,PI=3.1415926
contains

    subroutine fnf(L,MaxT)
        implicit none
        integer :: i,j,t,L,N_pulse=0,count_number=0,MaxT
        real :: epsilon=0.04,a=0.84,b=0.07,D=1.0,A1=2.5,f=4.5,effect,u(L),v(L),fun_u(L),fun_v(L),delta,stepT,&
                U_change(1:MaxT),T_spike=0.0,T_interval=0.0
        u=0.0
        v=0.0
        do t=1,MaxT,1
            do i=1,L,1
                stepT=h*t
                effect=0.0
                do j=1,L,1
                    effect=effect+ring_matrix(i,j)*(u(j)-u(i))
                end do
                if(i==1) then
                    delta=1.0
                else
                    delta=0.0
                end if
                if(i==1.and.stepT<100) then
                    fun_u(i)=u(i)+h*(1.0/epsilon*(u(i)*(1.0-u(i))*(u(i)-(v(i)+b)/a))+A1*sin(2.0*PI*f*stepT))
                else
                    fun_u(i)=u(i)+h*(1.0/epsilon*(u(i)*(1.0-u(i))*(u(i)-(v(i)+b)/a))+delta*A1*sin(2.0*PI*f*stepT)+D*effect)
                end if
                if(u(i)<1.0/3.0) then
                    fun_v(i)=v(i)+h*(-v(i))
                else if(u(i)>=1.0/3.0.and.u(i)<=1.0) then
                    fun_v(i)=v(i)+h*(1.0-6.75*u(i)*(u(i)-1.0)**2-v(i))
                else if(u(i)>1.0) then
                    fun_v(i)=v(i)+h*(1.0-v(i))
                end if
                if(i==2) then
                    U_change(t)=fun_u(i)
                    if(U_change(t)>0.5) then
                        count_number=count_number+1
                    else if(count_number>70.and.U_change(t)<0.2) then
                        N_pulse=N_pulse+1
                        if(N_pulse==14) then
                            T_spike=stepT/13.0
                            T_interval=MaxT*0.02-stepT
                        end if
                        count_number=0
                    end if
                end if
            end do
            u=fun_u
            v=fun_v
        end do
        write(20,*) L,T_spike
        write(40,*) L,T_interval
        N_pulse=0
        return
    end subroutine fnf
    !环
    subroutine ring(node_number)
        implicit none
        integer :: node_number,a
        allocate(ring_matrix(node_number,node_number))
        !初始化网络矩阵
        ring_matrix=0
        !赋值节点关系
        do a=2,node_number,1
            ring_matrix(a,a-1)=1
        end do
        !首尾相连
        ring_matrix(1,node_number)=1
    end subroutine ring

end module BOM


program main
    use BOM
    implicit none
    integer :: i,j=1,L
    real(kind=4) :: T_Link_Array(41)=(/106.879997,110.399994,113.899994,117.399994,120.860001,124.360001,127.879997,131.300003,&
            134.819992,138.339996,141.779999,145.279999,148.800003,152.259995,155.73999,159.259995,162.759995,166.199997,&
            169.699997,173.220001,176.660004,180.159988,183.679993,187.139999,190.619995,194.139999,197.639999,201.080002,&
            204.599991,208.099991,211.539993,215.059998,218.559998,222.019989,225.5,229.019989,232.519989,235.959991,&
            239.479996,243.0,246.419998/)
    open(20,file="L_T_Spike.txt")
    open(40,file="L_T_Interval.txt")
    do i=300,700,10
        call ring(i)
        call fnf(i,int(T_Link_Array(j)/0.02))
        deallocate(ring_matrix)
        j=j+1
    end do
    close(20)
end program main

figure 6_D

module BOM
    implicit none
    integer,allocatable :: ring_matrix(:,:)
    real,parameter :: h=0.02,PI=3.1415926
contains

    subroutine fnf(L,MaxT)
        implicit none
        integer :: i,j,t,L,N_pulse=0,count_number=0,MaxT
        real :: epsilon=0.04,a=0.84,b=0.07,D=1.0,A1=2.5,f=4.5,effect,u(L),v(L),fun_u(L),fun_v(L),delta,stepT,T_limit
        u=0.0
        v=0.0
        do t=1,MaxT,1
            do i=1,L,1
                stepT=h*t
                effect=0.0
                do j=1,L,1
                    effect=effect+ring_matrix(i,j)*(u(j)-u(i))
                end do
                if(i==1) then
                    delta=1.0
                else
                    delta=0.0
                end if
                if(i==1.and.stepT<100) then
                    fun_u(i)=u(i)+h*(1.0/epsilon*(u(i)*(1.0-u(i))*(u(i)-(v(i)+b)/a))+A1*sin(2.0*PI*f*stepT))
                else
                    fun_u(i)=u(i)+h*(1.0/epsilon*(u(i)*(1.0-u(i))*(u(i)-(v(i)+b)/a))+delta*A1*sin(2.0*PI*f*stepT)+D*effect)
                end if
                if(u(i)<1.0/3.0) then
                    fun_v(i)=v(i)+h*(-v(i))
                else if(u(i)>=1.0/3.0.and.u(i)<=1.0) then
                    fun_v(i)=v(i)+h*(1.0-6.75*u(i)*(u(i)-1.0)**2-v(i))
                else if(u(i)>1.0) then
                    fun_v(i)=v(i)+h*(1.0-v(i))
                end if
            end do
            u=fun_u
            v=fun_v
        end do
        T_limit=MaxT*h/14.0
        write(20,*) L,T_limit
        N_pulse=0
        return
    end subroutine fnf
    !环
    subroutine ring(node_number)
        implicit none
        integer :: node_number,a
        allocate(ring_matrix(node_number,node_number))
        !初始化网络矩阵
        ring_matrix=0
        !赋值节点关系
        do a=2,node_number,1
            ring_matrix(a,a-1)=1
        end do
        !首尾相连
        ring_matrix(1,node_number)=1
    end subroutine ring

end module BOM


program main
    use BOM
    implicit none
    integer :: i,j=1,L
    real(kind=4) :: T_Link_Array(41)=(/106.879997,110.399994,113.899994,117.399994,120.860001,124.360001,127.879997,131.300003,&
            134.819992,138.339996,141.779999,145.279999,148.800003,152.259995,155.73999,159.259995,162.759995,166.199997,&
            169.699997,173.220001,176.660004,180.159988,183.679993,187.139999,190.619995,194.139999,197.639999,201.080002,&
            204.599991,208.099991,211.539993,215.059998,218.559998,222.019989,225.5,229.019989,232.519989,235.959991,&
            239.479996,243.0,246.419998/)
    open(20,file="L_T_Limit.txt")
    do i=300,700,10
        call ring(i)
        call fnf(i,int(T_Link_Array(j)/0.02))
        deallocate(ring_matrix)
        j=j+1
    end do
    close(20)
end program main

• 红线和黑点同为数值结果

figure 7

figure 7_A

module BOM
    implicit none
    integer,parameter :: MaxT=10000,L=500
    integer,allocatable :: ring_matrix(:,:)
    real,parameter :: h=0.02,PI=3.1415926
contains

    subroutine fnf(a)
        implicit none
        real :: epsilon=0.04,a,b=0.07,D=1.0,A1=2.5,f=4.5,effect,u(L),v(L),fun_u(L),fun_v(L),delta,stepT,&
                U_change(1:MaxT),T_chain=0.0,T_peak(2),T_length=0.0
        integer :: i,j,t,L,N_pulse=0,count_number=0
        logical :: status
        u=0.0
        v=0.0
        status=.true.
        do t=1,MaxT,1
            do i=1,L,1
                stepT=h*t
                effect=0.0
                do j=1,L,1
                    effect=effect+ring_matrix(i,j)*(u(j)-u(i))
                end do
                if(i==1) then
                    delta=1.0
                else
                    delta=0.0
                end if
                if(i==1.and.status) then
                    fun_u(i)=u(i)+h*(1.0/epsilon*(u(i)*(1.0-u(i))*(u(i)-(v(i)+b)/a))+A1*sin(2.0*PI*f*stepT))
                else
                    fun_u(i)=u(i)+h*(1.0/epsilon*(u(i)*(1.0-u(i))*(u(i)-(v(i)+b)/a))+delta*A1*sin(2.0*PI*f*stepT)+D*effect)
                end if
                if(u(i)<1.0/3.0) then
                    fun_v(i)=v(i)+h*(-v(i))
                else if(u(i)>=1.0/3.0.and.u(i)<=1.0) then
                    fun_v(i)=v(i)+h*(1.0-6.75*u(i)*(u(i)-1.0)**2-v(i))
                else if(u(i)>1.0) then
                    fun_v(i)=v(i)+h*(1.0-v(i))
                end if
                if(i==1) then
                    U_change(t)=fun_u(i)
                    if(U_change(t)>0.5) then
                        count_number=count_number+1
                        if(N_pulse==2) then
                            T_length=stepT
                        end if
                    else if(count_number>=30.and.U_change(t)<0.2) then
                        N_pulse=N_pulse+1
                        count_number=0
                        if(N_pulse==1) then
                            T_peak(1)=stepT
                        end if
                        if(N_pulse==2) then
                            T_peak(2)=stepT
                            status=.false.
                        end if
                    end if
                end if
            end do
            u=fun_u
            v=fun_v
        end do
        T_chain=T_peak(2)-T_peak(1)
        write(20,*) a,T_chain,T_length,N_pulse
        N_pulse=0
        return
    end subroutine fnf
    !环
    subroutine ring(node_number)
        implicit none
        integer :: node_number,a
        allocate(ring_matrix(node_number,node_number))
        !初始化网络矩阵
        ring_matrix=0
        !赋值节点关系
        do a=2,node_number,1
            ring_matrix(a,a-1)=1
        end do
        !首尾相连
        ring_matrix(1,node_number)=1
    end subroutine ring

end module BOM


program main
    use BOM
    implicit none
    integer :: i
    real :: a
    open(20,file="a_Tchain_Tlength.txt")
    open(40,file="data.txt")
    call ring(L)
    do i=1,26,1
        a=0.59+i*0.01
        call fnf(a)
    end do
    close(20)
    deallocate(ring_matrix)
end program main

figure 7_B

module BOM
    implicit none
    integer,parameter :: L=500
    integer,allocatable :: ring_matrix(:,:)
    real,parameter :: h=0.02,PI=3.1415926
contains

    subroutine fnf(a,T_link,MaxT)
        implicit none
        integer :: i,j,t,N_pulse=0,count_number=0,MaxT
        real :: epsilon=0.04,a,b=0.07,D=1.0,A1=2.5,f=4.5,effect,u(L),v(L),fun_u(L),fun_v(L),delta,stepT,&
                U_change(1:MaxT),T_link
        u=0.0
        v=0.0
        do t=1,MaxT,1
            do i=1,L,1
                stepT=h*t
                effect=0.0
                do j=1,L,1
                    effect=effect+ring_matrix(i,j)*(u(j)-u(i))
                end do
                if(i==1) then
                    delta=1.0
                else
                    delta=0.0
                end if
                if(i==1.and.stepT<=T_link) then
                    fun_u(i)=u(i)+h*(1.0/epsilon*(u(i)*(1.0-u(i))*(u(i)-(v(i)+b)/a))+A1*sin(2.0*PI*f*stepT))
                else
                    fun_u(i)=u(i)+h*(1.0/epsilon*(u(i)*(1.0-u(i))*(u(i)-(v(i)+b)/a))+delta*A1*sin(2.0*PI*f*stepT)+D*effect)
                end if
                if(u(i)<1.0/3.0) then
                    fun_v(i)=v(i)+h*(-v(i))
                else if(u(i)>=1.0/3.0.and.u(i)<=1.0) then
                    fun_v(i)=v(i)+h*(1.0-6.75*u(i)*(u(i)-1.0)**2-v(i))
                else if(u(i)>1.0) then
                    fun_v(i)=v(i)+h*(1.0-v(i))
                end if
                if(i==1) then
                    U_change(t)=fun_u(i)
                    if(U_change(t)>0.5) then
                        count_number=count_number+1
                    else if(count_number>30.and.U_change(t)<0.2) then
                        N_pulse=N_pulse+1
                        count_number=0
                    end if
                end if
            end do
            u=fun_u
            v=fun_v
        end do
        write(20,*) a,N_pulse
        N_pulse=0
        return
    end subroutine fnf
    !环
    subroutine ring(node_number)
        implicit none
        integer :: node_number,a
        allocate(ring_matrix(node_number,node_number))
        !初始化网络矩阵
        ring_matrix=0
        !赋值节点关系
        do a=2,node_number,1
            ring_matrix(a,a-1)=1
        end do
        !首尾相连
        ring_matrix(1,node_number)=1
    end subroutine ring

end module BOM


program main
    use BOM
    implicit none
    integer :: i
    real :: a
    real(kind=4) :: T_Link_Array(26)=(/195.659988,194.339996,193.080002,191.940002,190.839996,189.860001,188.839996,&
            187.919998,187.019989,186.159988,185.319992,184.559998,183.779999,183.059998,182.37999,181.720001,&
            181.059998,180.419998,179.87999,179.279999,178.73999,178.179993,177.679993,177.179993,176.660004,176.199997/)
    open(20,file="a_Tchain_Tlength.txt")
    call ring(L)
    do i=1,26,1
        a=0.59+i*0.01
        call fnf(a,T_Link_Array(i),int(T_Link_Array(i)/0.02))
    end do
    close(20)
    deallocate(ring_matrix)
end program main

figure 7_C

module BOM
    implicit none
    integer,parameter :: L=500
    integer,allocatable :: ring_matrix(:,:)
    real,parameter :: h=0.02,PI=3.1415926
contains

    subroutine fnf(a,MaxT,number)
        implicit none
        integer :: i,j,t,N_pulse=0,count_number=0,MaxT,number
        real :: epsilon=0.04,a,b=0.07,D=1.0,A1=2.5,f=4.5,effect,u(L),v(L),fun_u(L),fun_v(L),delta,stepT,&
                U_change(1:MaxT),T_spike=0.0,T_interval=0.0
        logical :: status
        status=.true.
        u=0.0
        v=0.0
        do t=1,MaxT,1
            do i=1,L,1
                stepT=h*t
                effect=0.0
                do j=1,L,1
                    effect=effect+ring_matrix(i,j)*(u(j)-u(i))
                end do
                if(i==1) then
                    delta=1.0
                else
                    delta=0.0
                end if
                if(i==1.and.status) then
                    fun_u(i)=u(i)+h*(1.0/epsilon*(u(i)*(1.0-u(i))*(u(i)-(v(i)+b)/a))+A1*sin(2.0*PI*f*stepT))
                else
                    fun_u(i)=u(i)+h*(1.0/epsilon*(u(i)*(1.0-u(i))*(u(i)-(v(i)+b)/a))+delta*A1*sin(2.0*PI*f*stepT)+D*effect)
                end if
                if(u(i)<1.0/3.0) then
                    fun_v(i)=v(i)+h*(-v(i))
                else if(u(i)>=1.0/3.0.and.u(i)<=1.0) then
                    fun_v(i)=v(i)+h*(1.0-6.75*u(i)*(u(i)-1.0)**2-v(i))
                else if(u(i)>1.0) then
                    fun_v(i)=v(i)+h*(1.0-v(i))
                end if
                if(i==1) then
                    U_change(t)=fun_u(i)
                    if(U_change(t)>0.5) then
                        count_number=count_number+1
                    else if(count_number>30.and.U_change(t)<0.2) then
                        N_pulse=N_pulse+1
                        if(N_pulse==14) then
                            T_spike=stepT/13.0
                            T_interval=MaxT*0.02-stepT
                            status=.false.
                        end if
                        count_number=0
                    end if
                end if
!
            end do
            u=fun_u
            v=fun_v
        end do
        write(20,*) a,T_spike,T_interval,N_pulse
        N_pulse=0
        return
    end subroutine fnf
    !环
    subroutine ring(node_number)
        implicit none
        integer :: node_number,a
        allocate(ring_matrix(node_number,node_number))
        !初始化网络矩阵
        ring_matrix=0
        !赋值节点关系
        do a=2,node_number,1
            ring_matrix(a,a-1)=1
        end do
        !首尾相连
        ring_matrix(1,node_number)=1
    end subroutine ring

end module BOM


program main
    use BOM
    implicit none
    integer :: i
    real :: a
    real(kind=4) :: T_Link_Array(26)=(/195.659988,194.339996,193.080002,191.940002,190.839996,189.860001,188.839996,&
            187.919998,187.019989,186.159988,185.319992,184.559998,183.779999,183.059998,182.37999,181.720001,&
            181.059998,180.419998,179.87999,179.279999,178.73999,178.179993,177.679993,177.179993,176.660004,176.199997/)
    open(20,file="a_Tspike_Tinterval.txt")
    open(40,file="data.txt")
    call ring(L)
    do i=1,26,1
        a=0.59+i*0.01
        call fnf(a,int(T_Link_Array(i)/0.02),i)
    end do
    close(20)
    deallocate(ring_matrix)
end program main

figure 7_D

module BOM
    implicit none
    integer,parameter :: L=500
    integer,allocatable :: ring_matrix(:,:)
    real,parameter :: h=0.02,PI=3.1415926
contains

    subroutine fnf(a,MaxT,number)
        implicit none
        integer :: i,j,t,N_pulse=0,count_number=0,MaxT,number
        real :: epsilon=0.04,a,b=0.07,D=1.0,A1=2.5,f=4.5,effect,u(L),v(L),fun_u(L),fun_v(L),delta,stepT,&
                U_change(1:MaxT),T_limit
        u=0.0
        v=0.0
        do t=1,MaxT,1
            do i=1,L,1
                stepT=h*t
                effect=0.0
                do j=1,L,1
                    effect=effect+ring_matrix(i,j)*(u(j)-u(i))
                end do
                if(i==1) then
                    delta=1.0
                else
                    delta=0.0
                end if
                if(i==1.and.stepT<=141.5-1.5*number) then
                    fun_u(i)=u(i)+h*(1.0/epsilon*(u(i)*(1.0-u(i))*(u(i)-(v(i)+b)/a))+A1*sin(2.0*PI*f*stepT))
                else
                    fun_u(i)=u(i)+h*(1.0/epsilon*(u(i)*(1.0-u(i))*(u(i)-(v(i)+b)/a))+delta*A1*sin(2.0*PI*f*stepT)+D*effect)
                end if
                if(u(i)<1.0/3.0) then
                    fun_v(i)=v(i)+h*(-v(i))
                else if(u(i)>=1.0/3.0.and.u(i)<=1.0) then
                    fun_v(i)=v(i)+h*(1.0-6.75*u(i)*(u(i)-1.0)**2-v(i))
                else if(u(i)>1.0) then
                    fun_v(i)=v(i)+h*(1.0-v(i))
                end if
                if(i==2) then
                    U_change(t)=fun_u(i)
                    if(U_change(t)>0.5) then
                        count_number=count_number+1
                    else if(count_number>30.and.U_change(t)<0.2) then
                        N_pulse=N_pulse+1
                        if(N_pulse==14) then
                            T_limit=MaxT*0.02/14.0
                        end if
                        count_number=0
                    end if
                end if
            end do
            u=fun_u
            v=fun_v
        end do
        write(20,*) a,T_limit
        N_pulse=0
        return
    end subroutine fnf
    !环
    subroutine ring(node_number)
        implicit none
        integer :: node_number,a
        allocate(ring_matrix(node_number,node_number))
        !初始化网络矩阵
        ring_matrix=0
        !赋值节点关系
        do a=2,node_number,1
            ring_matrix(a,a-1)=1
        end do
        !首尾相连
        ring_matrix(1,node_number)=1
    end subroutine ring

end module BOM


program main
    use BOM
    implicit none
    integer :: i
    real :: a
    real(kind=4) :: T_Link_Array(26)=(/195.659988,194.339996,193.080002,191.940002,190.839996,189.860001,188.839996,&
            187.919998,187.019989,186.159988,185.319992,184.559998,183.779999,183.059998,182.37999,181.720001,&
            181.059998,180.419998,179.87999,179.279999,178.73999,178.179993,177.679993,177.179993,176.660004,176.199997/)
    open(20,file="a_Tspike_Tinterval.txt")
    call ring(L)
    do i=1,26,1
        a=0.59+i*0.01
        call fnf(a,int(T_Link_Array(i)/0.02),i)
    end do
    close(20)
    deallocate(ring_matrix)
end program main