非局域相互作用的双稳Fitzhugh-Nagumo振子环中的新型奇异态结构

New type of chimera structures in a ring of bistable FitzHugh–Nagumo oscillators with nonlocal interaction

阅读

下载地址：https://chinazhang-my.sharepoint.com/:b:/g/personal/rep_rebirth_zh_zhangwenhao_icu/ETsX-FYc0ANPi7tCgLHDDSsBO-rPXZxqfdNDTcRRa8lAhw?e=iixX5y

• 研究了双稳区非局域耦合Fitzhugh-Nagumo振子环的时空动力学。在无噪声网络中发现了一种新类型的奇异态图案，当孤立的元件不振荡时。探讨了当网络单元之间的耦合范围和耦合强度不同时，这些结构的存在范围。
• 本文研究了非局域相互作用对双稳区FHN振子环的时空动力学的影响。我们注意到，环的每个孤立(未耦合)元素都有两个稳定的平衡点，并且不振荡。我们的数值模拟表明，对于非局域耦合，不仅可以产生行波。我们还发现了一种特殊类型的奇异态模式，这与以前的研究不同。虽然环单元不是自维持振荡器，但观察到的奇异态结构可以在没有随机或其他外力的情况下存在。

复现

figure 1

• FHN oscillator model具有三个解，两个稳定焦点，一个鞍点。

module FHN
    implicit none
    real :: h=0.002,x,y
contains

    subroutine fnf()
        implicit none
        real :: func_x,func_y,epsilon=0.2,gamma=0.7,beta=0.0001
        func_x=x+h*((x-x**3.0/3.0-y)/epsilon)
        func_y=y+h*(gamma*x-y-beta)
        x=func_x
        y=func_y
        write(20,*) x,y
        return
    end subroutine fnf

end module FHN


program main
    use FHN
    implicit none
    integer :: i,a,b
    open(20,file="data_x_y.txt")
    do a=-10,-10,1
        do b=-5,15,1
            x=a*0.1-1.0
            y=b*0.1-0.5
            do i=1,4000,1
                call fnf()
            end do
        end do
    end do
    do a=-10,30,1
        do b=-5,-5,1
            x=a*0.1-1.0
            y=b*0.1-0.5
            do i=1,4000,1
                call fnf()
            end do
        end do
    end do
    do a=30,30,1
        do b=-5,15,1
            x=a*0.1-1.0
            y=b*0.1-0.5
            do i=1,4000,1
                call fnf()
            end do
        end do
    end do
    do a=-10,30,1
        do b=15,15,1
            x=a*0.1-1.0
            y=b*0.1-0.5
            do i=1,4000,1
                call fnf()
            end do
        end do
    end do
    do a=-10,10,1
        x=a*0.1
        y=x-(x**3)/3.0-0.08
        do i=1,4000,1
            call fnf()
        end do
    end do
    do a=-10,10,1
        x=a*0.1
        y=x-(x**3)/3.0+0.08
        do i=1,4000,1
            call fnf()
        end do
    end do
    close(20)
end program main

module FHN
    implicit none
    real :: h=0.002,x,y
contains

    subroutine fnf()
        implicit none
        real :: func_x,func_y,epsilon=0.2,gamma=0.7,beta=0.0001
        func_x=x+h*((x-x**3.0/3.0-y)/epsilon)
        func_y=y+h*(gamma*x-y-beta)
        x=func_x
        y=func_y
        return
    end subroutine fnf

end module FHN


program main
    use FHN
    implicit none
    integer :: i,a,b,status=0
    real :: x1,y1
    open(20,file="data_x_y_status.txt")
    do a=-200,200,1
        do b=-100,100,1
            x=a*0.01
            y=b*0.01
            x1=x
            y1=y
            do i=1,10000,1
                call fnf()
            end do
            if(abs(x-0.9488)<=0.01.and.abs(y-0.6641)<=0.01) then
                status=1
            elseif(abs(x+0.9485)<=0.01.and.abs(y+0.6641)<=0.01) then
                status=2
            end if
            write(20,*) x1,y1,status
        end do
    end do
    close(20)
end program main

figure 2

• 疑问（待解决）：在不知道figure 3的情况下，怎么画figure 2，知道figure 3，未解决如何框住figure 3的A,B,C,D四种情况，特别是A和D。

module FHN
    implicit none
    real,parameter :: h=0.001
    integer,parameter :: N=300,MaxT=10000
    integer,allocatable :: neighbour_matrix(:,:)
    real :: x(N),y(N)
contains

    subroutine x0_y0()
        implicit none
        integer :: i
        real :: x0,y0
        do i=1,N,1
            read(10,*) x(i),y(i)
        end do
    end subroutine x0_y0

    subroutine fnf(p,sigma)
        implicit none
        integer :: i,j,t,p
        real :: func_x(N),func_y(N),epsilon=0.2,gamma=0.7,beta=0.0001,sigma,effect
        do t=1,MaxT,1
            do i=1,N,1
                effect=0.0
                do j=1,N
                    effect=effect+neighbour_matrix(i,j)*(x(j)-x(i))
                end do
                func_x(i)=x(i)+h*((x(i)-(x(i)**3)/3.0-y(i)+epsilon*(sigma/(2.0*p))*effect)/epsilon)
                func_y(i)=y(i)+h*(gamma*x(i)-y(i)-beta)
                if(i==100.and.t==10000) then
                    write(20,*) p/real(N),sigma,x(i)
                    write(*,*) p/real(N),sigma,x(i)
                end if
                x=func_x
                y=func_y
            end do
        end do
        return
    end subroutine fnf

    !建立网状结构
    subroutine neighbour(K)
        implicit none
        integer :: K !单节点K个邻居
        integer :: i,j,number
        allocate(neighbour_matrix(N,N))
        !初始化网络矩阵
        neighbour_matrix=0
        do i=1,N-1,1
            do j=i+1,N,1
                if(abs(i-j)<=K/2.or.abs(i-j)>=N-K/2) then
                    neighbour_matrix(i,j)=1
                    neighbour_matrix(j,i)=1
                end if
            end do
        end do
    end subroutine neighbour

end module FHN


program main
    use FHN
    implicit none
    integer :: a,b,p
    real :: sigma
    open(10,file="data_x0_y0.txt")
    open(20,file="data_r_sigma_x.txt")
    call x0_y0()
    do a=1,30,1
        do b=1,100,1
            p=nint(N*(a*0.01))
            sigma=b*0.01
            call neighbour(2*p)
            call fnf(p,sigma)
            deallocate(neighbour_matrix)
        end do
    end do
    close(10)
    close(20)
end program main

figure 3

module FHN
    implicit none
    real,parameter :: h=0.001D0
    integer,parameter :: N=300,MaxT=300000
    integer,allocatable :: neighbour_matrix(:,:)
    real :: x(N),y(N)
contains

    subroutine x0_y0()
        implicit none
        integer :: i
        real :: x0,y0
        do i=1,N,1
            read(10,*) x(i),y(i)
        end do
        !do i=1,N,1
        !    call random_number(x(i))
        !    call random_number(y(i))
        !    x(i)=2.0*x(i)-1.0
        !    y(i)=2.0*y(i)-1.0
        !    write(10,*) x(i),y(i)
        !end do
    end subroutine x0_y0

    subroutine fnf(p,sigma)
        implicit none
        integer :: i,j,t,p
        real :: func_x(N),func_y(N),epsilon=0.2D0,gamma=0.7D0,beta=0.0001D0,sigma,effect
        do t=1,MaxT,1
            do i=1,N,1
                effect=0.0D0
                do j=1,N,1
                    effect=effect+neighbour_matrix(i,j)*(x(j)-x(i))
                end do
                func_x(i)=x(i)+h*((x(i)-(x(i)**3)/3.0-y(i)+epsilon*(sigma/(2.0*p))*effect)/epsilon)
                func_y(i)=y(i)+h*(gamma*x(i)-y(i)-beta)
                if(mod(t,50)==0.and.t>100000) then
                    write(20,*) i,t*h-100,x(i)
                    write(*,*) i,t*h-100,x(i)
                end if       
            end do
            x=func_x
            y=func_y
        end do
        return
    end subroutine fnf

    !建立网状结构
    subroutine neighbour(K)
        implicit none
        integer :: K !单节点K个邻居
        integer :: i,j,number
        allocate(neighbour_matrix(N,N))
        !初始化网络矩阵
        neighbour_matrix=0
        do i=1,N-1,1
            do j=i+1,N,1
                if(abs(i-j)<=K/2.or.abs(i-j)>=N-K/2) then
                    neighbour_matrix(i,j)=1
                    neighbour_matrix(j,i)=1
                end if
            end do
        end do
    end subroutine neighbour

end module FHN


program main
    use FHN
    implicit none
    integer :: p
    real :: sigma
    open(10,file="data_x0_y0.txt")
    open(20,file="data_i_t_x.txt")
    call random_seed()
    call x0_y0()
    p=25
    sigma=0.38D0
    call neighbour(2*p)
    call fnf(p,sigma)
    deallocate(neighbour_matrix)
    close(10)
    close(20)
end program main

figure 3_A

figure 3_B

figure 3_C

figure 3_D

figure 4

figure 4_A

module FHN
    implicit none
    real,parameter :: h=0.001D0
    integer,parameter :: N=300,MaxT=200000
    integer,allocatable :: neighbour_matrix(:,:)
    real :: x(N),y(N)
contains

    subroutine x0_y0()
        implicit none
        integer :: i
        real :: x0,y0
        do i=1,N,1
            read(10,*) x(i),y(i)
        end do
    end subroutine x0_y0

    subroutine fnf(p,sigma)
        implicit none
        integer :: i,j,t,p
        real :: func_x(N),func_y(N),epsilon=0.2D0,gamma=0.7D0,beta=0.0001D0,sigma,effect
        do t=1,MaxT,1
            do i=1,N,1
                effect=0.0D0
                do j=1,N,1
                    effect=effect+neighbour_matrix(i,j)*(x(j)-x(i))
                end do
                func_x(i)=x(i)+h*((x(i)-(x(i)**3)/3.0-y(i)+epsilon*(sigma/(2.0*p))*effect)/epsilon)
                func_y(i)=y(i)+h*(gamma*x(i)-y(i)-beta)
                if(t==200000) then
                    write(20,*) i,x(i)
                    write(*,*) i,x(i)
                end if       
            end do
            x=func_x
            y=func_y
        end do
        return
    end subroutine fnf

    !建立网状结构
    subroutine neighbour(K)
        implicit none
        integer :: K !单节点K个邻居
        integer :: i,j,number
        allocate(neighbour_matrix(N,N))
        !初始化网络矩阵
        neighbour_matrix=0
        do i=1,N-1,1
            do j=i+1,N,1
                if(abs(i-j)<=K/2.or.abs(i-j)>=N-K/2) then
                    neighbour_matrix(i,j)=1
                    neighbour_matrix(j,i)=1
                end if
            end do
        end do
    end subroutine neighbour

end module FHN


program main
    use FHN
    implicit none
    integer :: p
    real :: sigma
    open(10,file="data_x0_y0.txt")
    open(20,file="data_i_x.txt")
    call random_seed()
    call x0_y0()
    p=25
    sigma=0.34D0
    call neighbour(2*p)
    call fnf(p,sigma)
    deallocate(neighbour_matrix)
    close(10)
    close(20)
end program main

figure 4_B

module FHN
    implicit none
    real,parameter :: h=0.001D0
    integer,parameter :: N=300,MaxT=450000
    integer,allocatable :: neighbour_matrix(:,:)
    real :: x(N),y(N)
contains

    subroutine x0_y0()
        implicit none
        integer :: i
        real :: x0,y0
        do i=1,N,1
            read(10,*) x(i),y(i)
        end do
    end subroutine x0_y0

    subroutine fnf(p,sigma)
        implicit none
        integer :: i,j,t,p,a
        real :: func_x(N),func_y(N),epsilon=0.2D0,gamma=0.7D0,beta=0.0001D0,sigma,effect
        do t=1,MaxT,1
            do i=1,N,1
                effect=0.0D0
                do j=1,N,1
                    effect=effect+neighbour_matrix(i,j)*(x(j)-x(i))
                end do
                func_x(i)=x(i)+h*((x(i)-(x(i)**3)/3.0-y(i)+epsilon*(sigma/(2.0*p))*effect)/epsilon)
                func_y(i)=y(i)+h*(gamma*x(i)-y(i)-beta)
                do a=0,150000,5000
                    if(t==200000+a) then
                       write(20,*) i,x(i)
                       write(*,*) i,x(i)
                    end if
                end do         
            end do
            x=func_x
            y=func_y
        end do
        return
    end subroutine fnf

    !建立网状结构
    subroutine neighbour(K)
        implicit none
        integer :: K !单节点K个邻居
        integer :: i,j,number
        allocate(neighbour_matrix(N,N))
        !初始化网络矩阵
        neighbour_matrix=0
        do i=1,N-1,1
            do j=i+1,N,1
                if(abs(i-j)<=K/2.or.abs(i-j)>=N-K/2) then
                    neighbour_matrix(i,j)=1
                    neighbour_matrix(j,i)=1
                end if
            end do
        end do
    end subroutine neighbour

end module FHN


program main
    use FHN
    implicit none
    integer :: p
    real :: sigma
    open(10,file="data_x0_y0.txt")
    open(20,file="data_i_x.txt")
    call random_seed()
    call x0_y0()
    p=25
    sigma=0.34D0
    call neighbour(2*p)
    call fnf(p,sigma)
    deallocate(neighbour_matrix)
    close(10)
    close(20)
end program main

figure 4_C

module FHN
    implicit none
    real,parameter :: h=0.001D0
    integer,parameter :: N=300,MaxT=450000
    integer,allocatable :: neighbour_matrix(:,:)
    real :: x(N),y(N)
contains

    subroutine x0_y0()
        implicit none
        integer :: i
        real :: x0,y0
        do i=1,N,1
            read(10,*) x(i),y(i)
        end do
    end subroutine x0_y0

    subroutine fnf(p,sigma)
        implicit none
        integer :: i,j,t,p,a
        real :: func_x(N),func_y(N),epsilon=0.2D0,gamma=0.7D0,beta=0.0001D0,sigma,effect
        do t=1,MaxT,1
            do i=1,N,1
                effect=0.0D0
                do j=1,N,1
                    effect=effect+neighbour_matrix(i,j)*(x(j)-x(i))
                end do
                func_x(i)=x(i)+h*((x(i)-(x(i)**3)/3.0-y(i)+epsilon*(sigma/(2.0*p))*effect)/epsilon)
                func_y(i)=y(i)+h*(gamma*x(i)-y(i)-beta)
                if(i==180.and.t>150000) then
                    write(20,*) t*h,x(i)
                    write(*,*) t*h,x(i)
                end if       
            end do
            x=func_x
            y=func_y
        end do
        return
    end subroutine fnf

    !建立网状结构
    subroutine neighbour(K)
        implicit none
        integer :: K !单节点K个邻居
        integer :: i,j,number
        allocate(neighbour_matrix(N,N))
        !初始化网络矩阵
        neighbour_matrix=0
        do i=1,N-1,1
            do j=i+1,N,1
                if(abs(i-j)<=K/2.or.abs(i-j)>=N-K/2) then
                    neighbour_matrix(i,j)=1
                    neighbour_matrix(j,i)=1
                end if
            end do
        end do
    end subroutine neighbour

end module FHN


program main
    use FHN
    implicit none
    integer :: p
    real :: sigma
    open(10,file="data_x0_y0.txt")
    open(20,file="data_i_x.txt")
    call random_seed()
    call x0_y0()
    p=25
    sigma=0.34D0
    call neighbour(2*p)
    call fnf(p,sigma)
    deallocate(neighbour_matrix)
    close(10)
    close(20)
end program main

figure 4_D

module FHN
    implicit none
    real,parameter :: h=0.001D0
    integer,parameter :: N=300,MaxT=450000
    integer,allocatable :: neighbour_matrix(:,:)
    real :: x(N),y(N)
contains

    subroutine x0_y0()
        implicit none
        integer :: i
        real :: x0,y0
        do i=1,N,1
            read(10,*) x(i),y(i)
        end do
    end subroutine x0_y0

    subroutine fnf(p,sigma)
        implicit none
        integer :: i,j,t,p,a
        real :: func_x(N),func_y(N),epsilon=0.2D0,gamma=0.7D0,beta=0.0001D0,sigma,effect
        do t=1,MaxT,1
            do i=1,N,1
                effect=0.0D0
                do j=1,N,1
                    effect=effect+neighbour_matrix(i,j)*(x(j)-x(i))
                end do
                func_x(i)=x(i)+h*((x(i)-(x(i)**3)/3.0-y(i)+epsilon*(sigma/(2.0*p))*effect)/epsilon)
                func_y(i)=y(i)+h*(gamma*x(i)-y(i)-beta)
                do a=300000,300000,6000
                    if(i==180.and.t>150000.and.mod(t,50)==0) then
                        write(20,*) x(i),y(i)
                        write(*,*) t*h,x(i),y(i)
                    end if
                end do          
            end do
            x=func_x
            y=func_y
        end do
        return
    end subroutine fnf

    !建立网状结构
    subroutine neighbour(K)
        implicit none
        integer :: K !单节点K个邻居
        integer :: i,j,number
        allocate(neighbour_matrix(N,N))
        !初始化网络矩阵
        neighbour_matrix=0
        do i=1,N-1,1
            do j=i+1,N,1
                if(abs(i-j)<=K/2.or.abs(i-j)>=N-K/2) then
                    neighbour_matrix(i,j)=1
                    neighbour_matrix(j,i)=1
                end if
            end do
        end do
    end subroutine neighbour

end module FHN


program main
    use FHN
    implicit none
    integer :: p
    real :: sigma
    open(10,file="data_x0_y0.txt")
    open(20,file="data_i_x.txt")
    call random_seed()
    call x0_y0()
    p=25
    sigma=0.34D0
    call neighbour(2*p)
    call fnf(p,sigma)
    deallocate(neighbour_matrix)
    close(10)
    close(20)
end program main

figure 5

figure 5_A

module FHN
    implicit none
    real,parameter :: h=0.001D0
    integer,parameter :: N=300,MaxT=202000
    integer,allocatable :: neighbour_matrix(:,:)
    real :: x(N),y(N),delta_x(N)=0.0
contains

    subroutine x0_y0()
        implicit none
        integer :: i
        real :: x0,y0
        do i=1,N,1
            read(10,*) x(i),y(i)
        end do
    end subroutine x0_y0

    subroutine fnf(p,sigma)
        implicit none
        integer :: i,j,t,p
        real :: func_x(N),func_y(N),epsilon=0.2D0,gamma=0.7D0,beta=0.0001D0,sigma,effect
        do t=1,MaxT,1
            do i=1,N,1
                effect=0.0D0
                do j=1,N,1
                    effect=effect+neighbour_matrix(i,j)*(x(j)-x(i))
                end do
                func_x(i)=x(i)+h*((x(i)-(x(i)**3)/3.0-y(i)+epsilon*(sigma/(2.0*p))*effect)/epsilon)
                func_y(i)=y(i)+h*(gamma*x(i)-y(i)-beta)
                if(t==202000) then
                    write(40,*) i,x(i)
                end if
                if(t>200000.and.t<=202000) then
                    if(i==1) then
                        call delta_i(i,x(i),x(i+1),x(N))
                    else
                        call delta_i(i,x(i),x(i+1),x(i-1))
                    end if
                end if
                if(t==202000) then
                    write(20,*) i,delta_x(i)/2000.0
                    write(*,*) i,delta_x(i)/2000.0
                end if
            end do
            x=func_x
            y=func_y
        end do
        return
    end subroutine fnf

    subroutine delta_i(i,x1,x2,x3)
        implicit none
        integer :: i
        real :: x1,x2,x3
        delta_x(i)=delta_x(i)+(2.0*x1-x2-x3)**2
        return
    end subroutine delta_i

    !建立网状结构
    subroutine neighbour(K)
        implicit none
        integer :: K !单节点K个邻居
        integer :: i,j,number
        allocate(neighbour_matrix(N,N))
        !初始化网络矩阵
        neighbour_matrix=0
        do i=1,N-1,1
            do j=i+1,N,1
                if(abs(i-j)<=K/2.or.abs(i-j)>=N-K/2) then
                    neighbour_matrix(i,j)=1
                    neighbour_matrix(j,i)=1
                end if
            end do
        end do
    end subroutine neighbour

end module FHN


program main
    use FHN
    implicit none
    integer :: p
    real :: sigma
    open(10,file="data_x0_y0.txt")
    open(20,file="data_i_delta_i.txt")
    open(40,file="data_i_x.txt")
    call random_seed()
    call x0_y0()
    p=25
    sigma=0.38D0
    call neighbour(2*p)
    call fnf(p,sigma)
    deallocate(neighbour_matrix)
    close(10)
    close(20)
end program main

figure 5_B

module FHN
    implicit none
    real,parameter :: h=0.001D0
    integer,parameter :: N=300,MaxT=600000
    integer,allocatable :: neighbour_matrix(:,:)
    real :: x(N),y(N)
contains

    subroutine x0_y0()
        implicit none
        integer :: i
        real :: x0,y0
        do i=1,N,1
            read(10,*) x(i),y(i)
        end do
    end subroutine x0_y0

    subroutine fnf(p,sigma)
        implicit none
        integer :: i,j,t,a,p,count_number(N)=0,N_pulse(N)=0
        real :: func_x(N),func_y(N),epsilon=0.2D0,gamma=0.7D0,beta=0.0001D0,sigma,effect
        do t=1,MaxT,1
            do i=1,N,1
                effect=0.0D0
                do j=1,N,1
                    effect=effect+neighbour_matrix(i,j)*(x(j)-x(i))
                end do
                func_x(i)=x(i)+h*((x(i)-(x(i)**3)/3.0-y(i)+epsilon*(sigma/(2.0*p))*effect)/epsilon)
                func_y(i)=y(i)+h*(gamma*x(i)-y(i)-beta)
                if(t>100000.and.mod(t,50)==0) then
                    if(i==216) then
                        write(40,*) t*h,x(i)
                    end if
                    if(func_x(i)>1.0) then
                        count_number(i)=count_number(i)+1
                    else if(count_number(i)>20.and.func_x(i)<-1.0) then
                        N_pulse(i)=N_pulse(i)+1
                        count_number(i)=0
                    end if
                end if
            end do
            x=func_x
            y=func_y
            if(t==600000) then
                do a=1,N,1
                    write(20,*) a,N_pulse(a)/500.0
                    write(*,*) a,N_pulse(a)
                end do
            end if
        end do
        return
    end subroutine fnf


    !建立网状结构
    subroutine neighbour(K)
        implicit none
        integer :: K !单节点K个邻居
        integer :: i,j,number
        allocate(neighbour_matrix(N,N))
        !初始化网络矩阵
        neighbour_matrix=0
        do i=1,N-1,1
            do j=i+1,N,1
                if(abs(i-j)<=K/2.or.abs(i-j)>=N-K/2) then
                    neighbour_matrix(i,j)=1
                    neighbour_matrix(j,i)=1
                end if
            end do
        end do
    end subroutine neighbour

end module FHN


program main
    use FHN
    implicit none
    integer :: p
    real :: sigma
    open(10,file="data_x0_y0.txt")
    open(20,file="data_i_f_i.txt")
    open(40,file="data_216_t.txt")
    call random_seed()
    call x0_y0()
    p=25
    sigma=0.38D0
    call neighbour(2*p)
    call fnf(p,sigma)
    deallocate(neighbour_matrix)
    close(10)
    close(20)
end program main

• i=39,M(500)=44

• i=187,M(500)=16

• i=190,M(500)=38

• i=216,M(500)=3

figure 5_C

module FHN
    implicit none
    real,parameter :: h=0.01D0
    integer,parameter :: N=300,MaxT=20000,L=90
    integer,allocatable :: neighbour_matrix(:,:)
    real :: x(N),y(N)
contains

    subroutine x0_y0()
        implicit none
        integer :: i
        real :: x0,y0
        do i=1,N,1
            read(10,*) x(i),y(i)
        end do
    end subroutine x0_y0

    subroutine fnf(p,sigma)
        implicit none
        integer :: i,j,t,p,a
        real :: func_x(N),func_y(N),epsilon=0.2D0,gamma=0.7D0,beta=0.0001D0,sigma,effect,xi_average(N)=0.0&
                &,XX(N,MaxT),xi(N)
        do t=1,MaxT+1,1
            do i=1,N,1
                effect=0.0D0
                do j=1,N,1
                    effect=effect+neighbour_matrix(i,j)*(x(j)-x(i))
                end do
                func_x(i)=x(i)+h*((x(i)-(x(i)**3)/3.0-y(i)+epsilon*(sigma/(2.0*p))*effect)/epsilon)
                func_y(i)=y(i)+h*(gamma*x(i)-y(i)-beta)
!                if(t>100000.and.t<=105000) then
!                    xi=x
!                    xi_average(i)=xi_average(i)+x(i)
!                end if
!                if(t==105001) then
!                    xi_mean(i)=xi_mean(i)+xi(i)-xi_average(i)/real(t-100000)
!                    xn_mean=xn_mean+xi(L)-xi_average(L)/real(t-100000)
!                    xx=(xi_mean(i)*xn_mean)
!                    yy1=((xn_mean)**2)
!                    yy2=((xi_mean(i))**2)
!                    R_i(i)=xx/sqrt(yy1*yy2)
!                    write(20,*) i,R_i(i)
!                end if
                XX(i,t)=x(i)
            end do
            x=func_x
            y=func_y
        end do
        call Rn_i(XX)
        return
    end subroutine fnf


    !建立网状结构
    subroutine neighbour(K)
        implicit none
        integer :: K !单节点K个邻居
        integer :: i,j,number
        allocate(neighbour_matrix(N,N))
        !初始化网络矩阵
        neighbour_matrix=0
        do i=1,N-1,1
            do j=i+1,N,1
                if(abs(i-j)<=K/2.or.abs(i-j)>=N-K/2) then
                    neighbour_matrix(i,j)=1
                    neighbour_matrix(j,i)=1
                end if
            end do
        end do
        return
    end subroutine neighbour

    subroutine Rn_i(XX)
        implicit none
        integer :: i,T_mid=10000
        real :: XX(N,MaxT),x_aver(N,MaxT),R_i(N)
        x_aver=0.0
        R_i=0.0
        do i=1,N
            x_aver(i,T_mid:MaxT)=XX(i,T_mid:MaxT)-real(sum(XX(i,T_mid:MaxT)))/real(MaxT-T_mid+1)
            write(*,*) sum(XX(i,T_mid:MaxT))
        end do
        do i=1,N
            R_i(i)=(sum(x_aver(i,T_mid:MaxT)*x_aver(90,T_mid:MaxT)))/real(MaxT-T_mid+1)/&
                    &sqrt(sum(x_aver(i,T_mid:MaxT)**2)/real(MaxT-T_mid+1)*sum(x_aver(90,T_mid:MaxT)**2)/real(MaxT-T_mid+1))
            write(20,*) i,R_i(i)
        end do


        return
    end subroutine Rn_i

end module FHN



program main
    use FHN
    implicit none
    integer :: p
    real :: sigma
    open(10,file="data_x0_y0.txt")
    open(20,file="data_i_R.txt")
    call random_seed()
    call x0_y0()
    p=25
    sigma=0.38D0
    call neighbour(2*p)
    call fnf(p,sigma)
    deallocate(neighbour_matrix)
    close(10)
    close(20)
end program main

figure 6

figure 6_A

module FHN
    implicit none
    real,parameter :: h=0.001D0
    integer,parameter :: N=300,MaxT=100000
    integer,allocatable :: neighbour_matrix(:,:)
    real :: x(N),y(N)
contains

    subroutine x0_y0()
        implicit none
        integer :: i
        real :: x0,y0
      do i=1,N,1
            read(10,*) x(i),y(i)
        end do
    end subroutine x0_y0

    subroutine fnf(p,sigma)
        implicit none
        integer :: i,j,t,p
        real :: func_x(N),func_y(N),epsilon=0.2D0,gamma=0.7D0,beta=0.001D0,sigma,effect
        do t=1,MaxT,1
            do i=1,N,1
                effect=0.0D0
                do j=1,N,1
                    effect=effect+neighbour_matrix(i,j)*(x(j)-x(i))
                end do
                func_x(i)=x(i)+h*((x(i)-(x(i)**3)/3.0-y(i)+epsilon*(sigma/(2.0*p))*effect)/epsilon)
                func_y(i)=y(i)+h*(gamma*x(i)-y(i)-beta)
                if(t==100000) then
                    write(20,*) i,x(i)
                    write(*,*) i,x(i)
                end if       
            end do
            x=func_x
            y=func_y
        end do
        return
    end subroutine fnf

    !建立网状结构
    subroutine neighbour(K)
        implicit none
        integer :: K !单节点K个邻居
        integer :: i,j,number
        allocate(neighbour_matrix(N,N))
        !初始化网络矩阵
        neighbour_matrix=0
        do i=1,N-1,1
            do j=i+1,N,1
                if(abs(i-j)<=K/2.or.abs(i-j)>=N-K/2) then
                    neighbour_matrix(i,j)=1
                    neighbour_matrix(j,i)=1
                end if
            end do
        end do
    end subroutine neighbour

end module FHN


program main
    use FHN
    implicit none
    integer :: p
    real :: sigma
    open(10,file="data_x0_y0.txt")
    open(20,file="data_i_x.txt")
    call random_seed()
    call x0_y0()
    p=25
    sigma=0.14D0
    call neighbour(2*p)
    call fnf(p,sigma)
    deallocate(neighbour_matrix)
    close(10)
    close(20)
end program main

figure 6_B

module FHN
    implicit none
    real,parameter :: h=0.001D0
    integer,parameter :: N=300,MaxT=450000
    integer,allocatable :: neighbour_matrix(:,:)
    real :: x(N),y(N)
contains

    subroutine x0_y0()
        implicit none
        integer :: i
        real :: x0,y0
      do i=1,N,1
            read(10,*) x(i),y(i)
        end do
    end subroutine x0_y0

    subroutine fnf(p,sigma)
        implicit none
        integer :: i,j,t,p,a
        real :: func_x(N),func_y(N),epsilon=0.2D0,gamma=0.7D0,beta=0.001D0,sigma,effect
        do t=1,MaxT,1
            do i=1,N,1
                effect=0.0D0
                do j=1,N,1
                    effect=effect+neighbour_matrix(i,j)*(x(j)-x(i))
                end do
                func_x(i)=x(i)+h*((x(i)-(x(i)**3)/3.0-y(i)+epsilon*(sigma/(2.0*p))*effect)/epsilon)
                func_y(i)=y(i)+h*(gamma*x(i)-y(i)-beta)
                do a=0,150000,5000
                    if(t==200000+a) then
                       write(20,*) i,x(i)
                       write(*,*) i,x(i)
                    end if
                end do       
            end do
            x=func_x
            y=func_y
        end do
        return
    end subroutine fnf

    !建立网状结构
    subroutine neighbour(K)
        implicit none
        integer :: K !单节点K个邻居
        integer :: i,j,number
        allocate(neighbour_matrix(N,N))
        !初始化网络矩阵
        neighbour_matrix=0
        do i=1,N-1,1
            do j=i+1,N,1
                if(abs(i-j)<=K/2.or.abs(i-j)>=N-K/2) then
                    neighbour_matrix(i,j)=1
                    neighbour_matrix(j,i)=1
                end if
            end do
        end do
    end subroutine neighbour

end module FHN


program main
    use FHN
    implicit none
    integer :: p
    real :: sigma
    open(10,file="data_x0_y0.txt")
    open(20,file="data_i_x.txt")
    call random_seed()
    call x0_y0()
    p=25
    sigma=0.14D0
    call neighbour(2*p)
    call fnf(p,sigma)
    deallocate(neighbour_matrix)
    close(10)
    close(20)
end program main

figure 6_C

module FHN
    implicit none
    real,parameter :: h=0.001D0
    integer,parameter :: N=300,MaxT=450000
    integer,allocatable :: neighbour_matrix(:,:)
    real :: x(N),y(N)
contains

    subroutine x0_y0()
        implicit none
        integer :: i
        real :: x0,y0
      do i=1,N,1
            read(10,*) x(i),y(i)
        end do
    end subroutine x0_y0

    subroutine fnf(p,sigma)
        implicit none
        integer :: i,j,t,p
        real :: func_x(N),func_y(N),epsilon=0.2D0,gamma=0.7D0,beta=0.001D0,sigma,effect
        do t=1,MaxT,1
            do i=1,N,1
                effect=0.0D0
                do j=1,N,1
                    effect=effect+neighbour_matrix(i,j)*(x(j)-x(i))
                end do
                func_x(i)=x(i)+h*((x(i)-(x(i)**3)/3.0-y(i)+epsilon*(sigma/(2.0*p))*effect)/epsilon)
                func_y(i)=y(i)+h*(gamma*x(i)-y(i)-beta)
                if(i==221.and.t>150000) then
                    write(20,*) t*h-150,x(i)
                    write(*,*) t*h,x(i)
                end if       
            end do
            x=func_x
            y=func_y
        end do
        return
    end subroutine fnf

    !建立网状结构
    subroutine neighbour(K)
        implicit none
        integer :: K !单节点K个邻居
        integer :: i,j,number
        allocate(neighbour_matrix(N,N))
        !初始化网络矩阵
        neighbour_matrix=0
        do i=1,N-1,1
            do j=i+1,N,1
                if(abs(i-j)<=K/2.or.abs(i-j)>=N-K/2) then
                    neighbour_matrix(i,j)=1
                    neighbour_matrix(j,i)=1
                end if
            end do
        end do
    end subroutine neighbour

end module FHN


program main
    use FHN
    implicit none
    integer :: p
    real :: sigma
    open(10,file="data_x0_y0.txt")
    open(20,file="data_i_x.txt")
    call random_seed()
    call x0_y0()
    p=25
    sigma=0.14D0
    call neighbour(2*p)
    call fnf(p,sigma)
    deallocate(neighbour_matrix)
    close(10)
    close(20)
end program main

figure 6_D

module FHN
    implicit none
    real,parameter :: h=0.001D0
    integer,parameter :: N=300,MaxT=450000
    integer,allocatable :: neighbour_matrix(:,:)
    real :: x(N),y(N)
contains

    subroutine x0_y0()
        implicit none
        integer :: i
        real :: x0,y0
      do i=1,N,1
            read(10,*) x(i),y(i)
        end do
    end subroutine x0_y0

    subroutine fnf(p,sigma)
        implicit none
        integer :: i,j,t,p
        real :: func_x(N),func_y(N),epsilon=0.2D0,gamma=0.7D0,beta=0.001D0,sigma,effect
        do t=1,MaxT,1
            do i=1,N,1
                effect=0.0D0
                do j=1,N,1
                    effect=effect+neighbour_matrix(i,j)*(x(j)-x(i))
                end do
                func_x(i)=x(i)+h*((x(i)-(x(i)**3)/3.0-y(i)+epsilon*(sigma/(2.0*p))*effect)/epsilon)
                func_y(i)=y(i)+h*(gamma*x(i)-y(i)-beta)
                if(i==221.and.t>150000.and.mod(t,50)==0) then
                        write(20,*) x(i),y(i)
                        write(*,*) t*h,x(i),y(i)
                end if          
            end do
            x=func_x
            y=func_y
        end do
        return
    end subroutine fnf

    !建立网状结构
    subroutine neighbour(K)
        implicit none
        integer :: K !单节点K个邻居
        integer :: i,j,number
        allocate(neighbour_matrix(N,N))
        !初始化网络矩阵
        neighbour_matrix=0
        do i=1,N-1,1
            do j=i+1,N,1
                if(abs(i-j)<=K/2.or.abs(i-j)>=N-K/2) then
                    neighbour_matrix(i,j)=1
                    neighbour_matrix(j,i)=1
                end if
            end do
        end do
    end subroutine neighbour

end module FHN


program main
    use FHN
    implicit none
    integer :: p
    real :: sigma
    open(10,file="data_x0_y0.txt")
    open(20,file="data_i_x.txt")
    call random_seed()
    call x0_y0()
    p=25
    sigma=0.14D0
    call neighbour(2*p)
    call fnf(p,sigma)
    deallocate(neighbour_matrix)
    close(10)
    close(20)
end program main

figure 7

figure 7_A

module FHN
    implicit none
    real,parameter :: h=0.001D0
    integer,parameter :: N=300,MaxT=102000
    integer,allocatable :: neighbour_matrix(:,:)
    real :: x(N),y(N),delta_x(N)=0.0
contains

    subroutine x0_y0()
        implicit none
        integer :: i
        real :: x0,y0
        do i=1,N,1
            read(10,*) x(i),y(i)
        end do
    end subroutine x0_y0

    subroutine fnf(p,sigma)
        implicit none
        integer :: i,j,t,p
        real :: func_x(N),func_y(N),epsilon=0.2D0,gamma=0.7D0,beta=0.0001D0,sigma,effect
        do t=1,MaxT,1
            do i=1,N,1
                effect=0.0D0
                do j=1,N,1
                    effect=effect+neighbour_matrix(i,j)*(x(j)-x(i))
                end do
                func_x(i)=x(i)+h*((x(i)-(x(i)**3)/3.0-y(i)+epsilon*(sigma/(2.0*p))*effect)/epsilon)
                func_y(i)=y(i)+h*(gamma*x(i)-y(i)-beta)
                if(t==102000) then
                    write(40,*) i,x(i)
                end if
                if(t>100000.and.t<=102000) then
                    if(i==1) then
                        call delta_i(i,x(i),x(i+1),x(N))
                    else
                        call delta_i(i,x(i),x(i+1),x(i-1))
                    end if
                end if
                if(t==102000) then
                    write(20,*) i,delta_x(i)/2000.0
                    write(*,*) i,delta_x(i)/2000.0
                end if
            end do
            x=func_x
            y=func_y
        end do
        return
    end subroutine fnf

    subroutine delta_i(i,x1,x2,x3)
        implicit none
        integer :: i
        real :: x1,x2,x3
        delta_x(i)=delta_x(i)+(2.0*x1-x2-x3)**2
        return
    end subroutine delta_i

    !建立网状结构
    subroutine neighbour(K)
        implicit none
        integer :: K !单节点K个邻居
        integer :: i,j,number
        allocate(neighbour_matrix(N,N))
        !初始化网络矩阵
        neighbour_matrix=0
        do i=1,N-1,1
            do j=i+1,N,1
                if(abs(i-j)<=K/2.or.abs(i-j)>=N-K/2) then
                    neighbour_matrix(i,j)=1
                    neighbour_matrix(j,i)=1
                end if
            end do
        end do
    end subroutine neighbour

end module FHN


program main
    use FHN
    implicit none
    integer :: p
    real :: sigma
    open(10,file="data_x0_y0.txt")
    open(20,file="data_i_delta_i.txt")
    open(40,file="data_i_x.txt")
    call random_seed()
    call x0_y0()
    p=25
    sigma=0.14D0
    call neighbour(2*p)
    call fnf(p,sigma)
    deallocate(neighbour_matrix)
    close(10)
    close(20)
end program main

figure 7_B

module FHN
    implicit none
    real,parameter :: h=0.001D0
    integer,parameter :: N=300,MaxT=600000
    integer,allocatable :: neighbour_matrix(:,:)
    real :: x(N),y(N)
contains

    subroutine x0_y0()
        implicit none
        integer :: i
        real :: x0,y0
        do i=1,N,1
            read(10,*) x(i),y(i)
        end do
    end subroutine x0_y0

    subroutine fnf(p,sigma)
        implicit none
        integer :: i,j,t,a,p,count_number(N)=0,N_pulse(N)=0
        real :: func_x(N),func_y(N),epsilon=0.2D0,gamma=0.7D0,beta=0.0001D0,sigma,effect
        do t=1,MaxT,1
            do i=1,N,1
                effect=0.0D0
                do j=1,N,1
                    effect=effect+neighbour_matrix(i,j)*(x(j)-x(i))
                end do
                func_x(i)=x(i)+h*((x(i)-(x(i)**3)/3.0-y(i)+epsilon*(sigma/(2.0*p))*effect)/epsilon)
                func_y(i)=y(i)+h*(gamma*x(i)-y(i)-beta)
                if(t>100000.and.mod(t,50)==0) then
                    if(i==216) then
                        write(40,*) t*h,x(i)
                    end if
                    if(func_x(i)>1.0) then
                        count_number(i)=count_number(i)+1
                    else if(count_number(i)>20.and.func_x(i)<-1.0) then
                        N_pulse(i)=N_pulse(i)+1
                        count_number(i)=0
                    end if
                end if
            end do
            x=func_x
            y=func_y
            if(t==600000) then
                do a=1,N,1
                    write(20,*) a,N_pulse(a)/500.0
                    write(*,*) a,N_pulse(a)
                end do
            end if
        end do
        return
    end subroutine fnf


    !建立网状结构
    subroutine neighbour(K)
        implicit none
        integer :: K !单节点K个邻居
        integer :: i,j,number
        allocate(neighbour_matrix(N,N))
        !初始化网络矩阵
        neighbour_matrix=0
        do i=1,N-1,1
            do j=i+1,N,1
                if(abs(i-j)<=K/2.or.abs(i-j)>=N-K/2) then
                    neighbour_matrix(i,j)=1
                    neighbour_matrix(j,i)=1
                end if
            end do
        end do
    end subroutine neighbour

end module FHN


program main
    use FHN
    implicit none
    integer :: p
    real :: sigma
    open(10,file="data_x0_y0.txt")
    open(20,file="data_i_f_i.txt")
    open(40,file="data_216_t.txt")
    call random_seed()
    call x0_y0()
    p=25
    sigma=0.14D0
    call neighbour(2*p)
    call fnf(p,sigma)
    deallocate(neighbour_matrix)
    close(10)
    close(20)
end program main

• i=216,M(500)=58

figure 7_C

module FHN
    implicit none
    real,parameter :: h=0.01D0
    integer,parameter :: N=300,MaxT=20000,L=90
    integer,allocatable :: neighbour_matrix(:,:)
    real :: x(N),y(N)
contains

    subroutine x0_y0()
        implicit none
        integer :: i
        real :: x0,y0
        do i=1,N,1
            read(10,*) x(i),y(i)
        end do
    end subroutine x0_y0

    subroutine fnf(p,sigma)
        implicit none
        integer :: i,j,t,p,a
        real :: func_x(N),func_y(N),epsilon=0.2D0,gamma=0.7D0,beta=0.0001D0,sigma,effect,xi_average(N)=0.0&
                &,XX(N,MaxT),xi(N)
        do t=1,MaxT+1,1
            do i=1,N,1
                effect=0.0D0
                do j=1,N,1
                    effect=effect+neighbour_matrix(i,j)*(x(j)-x(i))
                end do
                func_x(i)=x(i)+h*((x(i)-(x(i)**3)/3.0-y(i)+epsilon*(sigma/(2.0*p))*effect)/epsilon)
                func_y(i)=y(i)+h*(gamma*x(i)-y(i)-beta)
!                if(t>100000.and.t<=105000) then
!                    xi=x
!                    xi_average(i)=xi_average(i)+x(i)
!                end if
!                if(t==105001) then
!                    xi_mean(i)=xi_mean(i)+xi(i)-xi_average(i)/real(t-100000)
!                    xn_mean=xn_mean+xi(L)-xi_average(L)/real(t-100000)
!                    xx=(xi_mean(i)*xn_mean)
!                    yy1=((xn_mean)**2)
!                    yy2=((xi_mean(i))**2)
!                    R_i(i)=xx/sqrt(yy1*yy2)
!                    write(20,*) i,R_i(i)
!                end if
                XX(i,t)=x(i)
            end do
            x=func_x
            y=func_y
        end do
        call Rn_i(XX)
        return
    end subroutine fnf


    !建立网状结构
    subroutine neighbour(K)
        implicit none
        integer :: K !单节点K个邻居
        integer :: i,j,number
        allocate(neighbour_matrix(N,N))
        !初始化网络矩阵
        neighbour_matrix=0
        do i=1,N-1,1
            do j=i+1,N,1
                if(abs(i-j)<=K/2.or.abs(i-j)>=N-K/2) then
                    neighbour_matrix(i,j)=1
                    neighbour_matrix(j,i)=1
                end if
            end do
        end do
        return
    end subroutine neighbour

    subroutine Rn_i(XX)
        implicit none
        integer :: i,T_mid=10000
        real :: XX(N,MaxT),x_aver(N,MaxT),R_i(N)
        x_aver=0.0
        R_i=0.0
        do i=1,N
            x_aver(i,T_mid:MaxT)=XX(i,T_mid:MaxT)-real(sum(XX(i,T_mid:MaxT)))/real(MaxT-T_mid+1)
            write(*,*) sum(XX(i,T_mid:MaxT))
        end do
        do i=1,N
            R_i(i)=(sum(x_aver(i,T_mid:MaxT)*x_aver(90,T_mid:MaxT)))/real(MaxT-T_mid+1)/&
                    &sqrt(sum(x_aver(i,T_mid:MaxT)**2)/real(MaxT-T_mid+1)*sum(x_aver(90,T_mid:MaxT)**2)/real(MaxT-T_mid+1))
            write(20,*) i,R_i(i)
        end do


        return
    end subroutine Rn_i

end module FHN



program main
    use FHN
    implicit none
    integer :: p
    real :: sigma
    open(10,file="data_x0_y0.txt")
    open(20,file="data_i_R.txt")
    call random_seed()
    call x0_y0()
    p=25
    sigma=0.38D0
    call neighbour(2*p)
    call fnf(p,sigma)
    deallocate(neighbour_matrix)
    close(10)
    close(20)
end program main

figure 8

figure 8_A

module FHN
    implicit none
    real,parameter :: h=0.001D0
    integer,parameter :: N=300,MaxT=450000
    integer,allocatable :: neighbour_matrix(:,:)
    real :: x(N),y(N)
contains

    subroutine x0_y0()
        implicit none
        integer :: i
        real :: x0,y0
        do i=1,N,1
            read(10,*) x(i),y(i)
        end do
    end subroutine x0_y0

    subroutine fnf(p,sigma)
        implicit none
        integer :: i,j,t,p,a
        real :: func_x(N),func_y(N),epsilon=0.2D0,gamma=0.7D0,beta=0.0001D0,sigma,effect
        do t=1,MaxT,1
            do i=1,N,1
                effect=0.0D0
                do j=1,N,1
                    effect=effect+neighbour_matrix(i,j)*(x(j)-x(i))
                end do
                func_x(i)=x(i)+h*((x(i)-(x(i)**3)/3.0-y(i)+epsilon*(sigma/(2.0*p))*effect)/epsilon)
                func_y(i)=y(i)+h*(gamma*x(i)-y(i)-beta)
                do a=0,150000,5000
                    if(t==200000+a) then
                       write(20,*) i,x(i)
                       write(*,*) i,x(i)
                    end if
                end do            
            end do
            x=func_x
            y=func_y
        end do
        return
    end subroutine fnf

    !建立网状结构
    subroutine neighbour(K)
        implicit none
        integer :: K !单节点K个邻居
        integer :: i,j,number
        allocate(neighbour_matrix(N,N))
        !初始化网络矩阵
        neighbour_matrix=0
        do i=1,N-1,1
            do j=i+1,N,1
                if(abs(i-j)<=K/2.or.abs(i-j)>=N-K/2) then
                    neighbour_matrix(i,j)=1
                    neighbour_matrix(j,i)=1
                end if
            end do
        end do
    end subroutine neighbour

end module FHN


program main
    use FHN
    implicit none
    integer :: p
    real :: sigma
    open(10,file="data_x0_y0.txt")
    open(20,file="data_i_x.txt")
    call random_seed()
    call x0_y0()
    p=15
    sigma=0.19D0
    call neighbour(2*p)
    call fnf(p,sigma)
    deallocate(neighbour_matrix)
    close(10)
    close(20)
end program main

figure 8_B

module FHN
    implicit none
    real,parameter :: h=0.001D0
    integer,parameter :: N=300,MaxT=450000
    integer,allocatable :: neighbour_matrix(:,:)
    real :: x(N),y(N)
contains

    subroutine x0_y0()
        implicit none
        integer :: i
        real :: x0,y0
      do i=1,N,1
            read(10,*) x(i),y(i)
        end do
    end subroutine x0_y0

    subroutine fnf(p,sigma)
        implicit none
        integer :: i,j,t,p,a
        real :: func_x(N),func_y(N),epsilon=0.2D0,gamma=0.7D0,beta=0.0001D0,sigma,effect
        do t=1,MaxT,1
            do i=1,N,1
                effect=0.0D0
                do j=1,N,1
                    effect=effect+neighbour_matrix(i,j)*(x(j)-x(i))
                end do
                func_x(i)=x(i)+h*((x(i)-(x(i)**3)/3.0-y(i)+epsilon*(sigma/(2.0*p))*effect)/epsilon)
                func_y(i)=y(i)+h*(gamma*x(i)-y(i)-beta)
                do a=0,150000,5000
                    if(t==200000+a) then
                       write(20,*) i,x(i)
                       write(*,*) i,x(i)
                    end if
                end do            
            end do
            x=func_x
            y=func_y
        end do
        return
    end subroutine fnf

    !建立网状结构
    subroutine neighbour(K)
        implicit none
        integer :: K !单节点K个邻居
        integer :: i,j,number
        allocate(neighbour_matrix(N,N))
        !初始化网络矩阵
        neighbour_matrix=0
        do i=1,N-1,1
            do j=i+1,N,1
                if(abs(i-j)<=K/2.or.abs(i-j)>=N-K/2) then
                    neighbour_matrix(i,j)=1
                    neighbour_matrix(j,i)=1
                end if
            end do
        end do
    end subroutine neighbour

end module FHN


program main
    use FHN
    implicit none
    integer :: p
    real :: sigma
    open(10,file="data_x0_y0.txt")
    open(20,file="data_i_x.txt")
    call random_seed()
    call x0_y0()
    p=15
    sigma=0.19D0
    call neighbour(2*p)
    call fnf(p,sigma)
    deallocate(neighbour_matrix)
    close(10)
    close(20)
end program main