双稳机制的嵌合态

Bistability-mechanism induced chimera state in one-dimensional paced excitable ring with nonlocal couplings

阅读

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

复现

figure 1

• u、v初始值：https://chinazhang-my.sharepoint.com/:f:/g/personal/rep_rebirth_zh_zhangwenhao_icu/EmQf09T44ktCrMuKjAYl2rIBxKiALnDCv6W48ns0834-mw?e=wOCtFx

module BOM
    implicit none
    integer,parameter :: R=33,MaxT=50000,L=200
    integer,allocatable :: neighbour_matrix(:,:)
    real,parameter :: h=0.02,PI=3.1415926
contains

    subroutine fnf()
        implicit none
        integer :: i,j,t,num
        real :: epsilon=0.04,a=0.84,b=0.07,A1=0.01,f=4.0,sigma=0.132,effect,u(L),v(L),fun_u(L),fun_v(L),x,y,Z(L),U_max(L),U_min(L)
        U_max=-10.0
        U_min=10.0
        !do i=1,L,1
        !    call random_number(x)
        !    call random_number(y)
        !    u(i)=2.0*x
        !    v(i)=2.0*y
        !    write(10,*) u(i),v(i)
        !end do
        do i=1,L,1
            read(10,*) u(i),v(i)
        end do
        do t=1,MaxT,1
            do i=1,L,1
                effect=0.0
                do j=1,L
                    effect=effect+neighbour_matrix(i,j)*(u(j)-u(i))
                end do
                fun_u(i)=u(i)+h*(1.0/epsilon*u(i)*(1.0-u(i))*(u(i)-(v(i)+b)/a)+sigma/(2.0*R)*effect+A1*sin(2.0*PI*f*t*h))
                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
                call parameter_Z(u,v,i,Z)
                if(mod(t,50)==0) then
                    write(20,*) i,t*h,u(i)
                    write(40,*) i,t*h,Z(i)
                end if
            end do
            call U_max_min(u,U_max,U_min)
            u=fun_u
            v=fun_v
        end do
        do j=1,L,1
            !if(U_max(j)>1.2) then
            !    write(60,*) j,U_max(j)
            !    write(80,*) j,U_min(j)
            !end if
            write(60,*) j,U_max(j)
            write(80,*) j,U_min(j)
        end do
        return
    end subroutine fnf

    !局域有序参数Z
    subroutine parameter_Z(u,v,i,Z)
        implicit none
        real :: u(L),v(L),Z(L),Ans_Z_X,Ans_Z_Y
        integer :: i,j
        Ans_Z_X=0.0
        Ans_Z_Y=0.0
        do j=1,L,1
            Ans_Z_X=Ans_Z_X+neighbour_matrix(i,j)*(cos(atan(v(j)/u(j))))
            Ans_Z_Y=Ans_Z_Y+neighbour_matrix(i,j)*(sin(atan(v(j)/u(j))))
        end do
        Z(i)=sqrt((Ans_Z_X/(2.0*R))**2+(Ans_Z_Y/(2.0*R))**2)
        return
    end subroutine parameter_Z

    !建立网状结构
    subroutine neighbour(N,K)
        implicit none
        integer :: N,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 BOM


program main
    use BOM
    implicit none
    open(10,file="u0_article.txt",access='sequential',position = 'rewind',action = 'read')
    open(20,file="data_u.txt")
    open(40,file="data_z.txt")
    open(60,file="data_u_max.txt")
    open(80,file="data_u_min.txt")
    call neighbour(L,2*R)
    call random_seed()
    call fnf()
    deallocate(neighbour_matrix)
    close(20)
end program main

figure 1_A

figure 1_B

figure 1_C

module BOM
    implicit none
    integer,parameter :: R=33,MaxT=50000,L=200
    integer,allocatable :: neighbour_matrix(:,:)
    real,parameter :: h=0.02,PI=3.1415926
contains

    subroutine fnf()
        implicit none
        integer :: i,j,t,num
        real :: epsilon=0.04,a=0.84,b=0.07,A1=0.01,f=4.0,sigma=0.132,effect,u(L),v(L),fun_u(L),fun_v(L),x,y,Z(L),U_max(L),U_min(L)
        U_max=-10.0
        U_min=10.0
        do i=1,L,1
            read(10,*) u(i),v(i)
        end do
        do t=1,MaxT,1
            do i=1,L,1
                effect=0.0
                do j=1,L
                    effect=effect+neighbour_matrix(i,j)*(u(j)-u(i))
                end do
                fun_u(i)=u(i)+h*(1.0/epsilon*u(i)*(1.0-u(i))*(u(i)-(v(i)+b)/a)+sigma/(2.0*R)*effect+A1*sin(2.0*PI*f*t*h))
                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
            end do
            if(t>=45000) then
                call U_max_min(u,U_max,U_min)
            end if
            u=fun_u
            v=fun_v
        end do
        do j=1,L,1
            if(U_max(j)<0.2) then
                write(60,*) j,U_max(j)
                write(80,*) j,U_min(j)
            end if
        end do
        return
    end subroutine fnf


    !Umax与Umin
    subroutine U_max_min(u,U_max,U_min)
        implicit none
        real :: U_max(L),U_min(L),u(L)
        integer :: i
        do i=1,L,1
            if(U_max(i)<u(i)) then
                U_max(i)=u(i)
            end if
            if(U_min(i)>u(i)) then
                U_min(i)=u(i)
            end if
        end do
        return
    end subroutine U_max_min


    !建立网状结构
    subroutine neighbour(N,K)
        implicit none
        integer :: N,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 BOM


program main
    use BOM
    implicit none
    open(10,file="u0_article.txt",access='sequential',position = 'rewind',action = 'read')
    open(60,file="data_u_max.txt")
    open(80,file="data_u_min.txt")
    call neighbour(L,2*R)
    call random_seed()
    call fnf()
    deallocate(neighbour_matrix)
    close(20)
end program main

figure 1_D

module BOM
    implicit none
    integer,parameter :: R=33,MaxT=50000,L=200
    integer,allocatable :: neighbour_matrix(:,:)
    real,parameter :: h=0.02,PI=3.1415926
contains

    subroutine fnf()
        implicit none
        integer :: i,j,t,M(L),stepT,count_number(L)
        real :: epsilon=0.04,a=0.84,b=0.07,A1=0.01,f=4.0,sigma=0.132,effect,u(L),v(L),fun_u(L),fun_v(L)
        M(L)=0
        count_number(L)=0
        do i=1,L,1
            read(10,*) u(i),v(i)
        end do
        do t=1,MaxT,1
            do i=1,L,1
                effect=0.0
                do j=1,L
                    effect=effect+neighbour_matrix(i,j)*(u(j)-u(i))
                end do
                fun_u(i)=u(i)+h*(1.0/epsilon*u(i)*(1.0-u(i))*(u(i)-(v(i)+b)/a)+sigma/(2.0*R)*effect+A1*sin(2.0*PI*f*t*h))
                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(t>=1) then
                    if(u(i)>0.2) then
                        count_number(i)=count_number(i)+1
                    else if(count_number(i)>30.and.u(i)<0.2) then
                        M(i)=M(i)+1
                        count_number(i)=0
                    end if    
                    if(i==4) then
                        continue
                    end if    
                end if
            end do
            u=fun_u
            v=fun_v
        end do
        stepT=MaxT
        call parameter_Omega(M,stepT)
        return
    end subroutine fnf


    !平均相速度omega
    subroutine parameter_Omega(M,stepT)
        implicit none
        real :: omega(L)
        integer :: i,stepT,M(L)
        do i=1,L,1
            omega(i)=(2.0*PI*M(i))/(stepT*h)
            write(20,*) i,omega(i)
        end do
        return
    end subroutine parameter_Omega

    !建立网状结构
    subroutine neighbour(N,K)
        implicit none
        integer :: N,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 BOM


program main
    use BOM
    implicit none
    open(10,file="u0_article.txt",access='sequential',position = 'rewind',action = 'read')
    open(20,file="data.txt")
    open(60,file="data_M.txt")
    open(80,file="data_u_min.txt")
    call neighbour(L,2*R)
    call random_seed()
    call fnf()
    deallocate(neighbour_matrix)
    close(20)
end program main

figure 2

module BOM
    implicit none
    integer,parameter :: R=33,MaxT=50000,L=200
    integer,allocatable :: neighbour_matrix(:,:)
    real,parameter :: h=0.02,PI=3.1415926
contains

    subroutine fnf()
        implicit none
        integer :: i,j,t,num
        real :: epsilon=0.04,a=0.84,b=0.07,A1=0.01,f=4.0,sigma=0.132,effect,u(L),v(L),fun_u(L),fun_v(L),x,y,Z(L),U_max(L),U_min(L)
        U_max=-10.0
        U_min=10.0
        do i=1,L,1
            read(10,*) u(i),v(i)
        end do
        do t=1,MaxT,1
            do i=1,L,1
                effect=0.0
                do j=1,L
                    effect=effect+neighbour_matrix(i,j)*(u(j)-u(i))
                end do
                fun_u(i)=u(i)+h*(1.0/epsilon*u(i)*(1.0-u(i))*(u(i)-(v(i)+b)/a)+sigma/(2.0*R)*effect+A1*sin(2.0*PI*f*t*h))
                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==140.and.t>=45000) then
                    write(20,*) t*h,u(i)
                    write(40,*) t*h,v(i)
                end if
            end do
            u=fun_u
            v=fun_v
        end do
        return
    end subroutine fnf

    !建立网状结构
    subroutine neighbour(N,K)
        implicit none
        integer :: N,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 BOM


program main
    use BOM
    implicit none
    open(10,file="u0_article.txt",access='sequential',position = 'rewind',action = 'read')
    open(20,file="data_u.txt")
    open(40,file="data_v.txt")
    call neighbour(L,2*R)
    call random_seed()
    call fnf()
    deallocate(neighbour_matrix)
    close(20)
    close(40)
end program main

figure 2_A

figure 2_B

figure 2_C

figure 2_D

figure 2_E

figure 2_F

figure 3

figure 3_A

figure 3_B

module BOM
    implicit none
    real,parameter :: PI=3.1415926
contains

    subroutine fnf(u,v)
        implicit none
        integer :: i,j,t,num
        real :: epsilon=0.04,a=0.84,b=0.07,fun_u,fun_v,u,v
        real :: theta,tantheta,sintheta,radius
        fun_u=1.0/epsilon*u*(1.0-u)*(u-(v+b)/a)
        if(u<1.0/3.0)then
            fun_v=0.0-v
        else if(u>=1.0/3.0.and.u<=1)then
            fun_v=1-6.75*u*(u-1)**2-v
        else if(u>1.0)then
            fun_v=1.0-v
        end if
        tantheta=fun_v/fun_u !tantheta>0，向量方向指向一、三象限
        radius=sqrt(fun_u**2+fun_v**2)
        sintheta=fun_v/radius !sintheta>0,向量方向指向一、二象限
        if(sintheta>=0.0) then
            if(tantheta>=0.0) then
                theta=atan(tantheta) !指向第一象限
            else
                theta=atan(tantheta)+PI !指向第二象限
            end if
        else
            if(tantheta>=0) then
                theta=atan(tantheta)+PI !指向第三象限
            else
                theta=atan(tantheta) !指向第四象限
            end if
        end if
        write(20,*) u,v,theta,1
        return
    end subroutine fnf


end module BOM


program main
    use BOM
    implicit none
    integer :: a,b
    open(20,file="data_vector.txt")
    do a=-5,29,1
        do b=-5,25,1
            call fnf(a*0.05,b*0.05)
        end do
    end do
    close(20)
end program main

figure 4

module BOM
    implicit none
    integer,parameter :: R=33,MaxT=1500,L=200
    integer,allocatable :: neighbour_matrix(:,:)
    real,parameter :: h=0.02,PI=3.1415926
contains

    subroutine fnf()
        implicit none
        integer :: i,j,t,num
        real :: epsilon=0.04,a=0.84,b=0.07,A1=0.01,f=4.0,sigma=0.132,effect,u(L),v(L),fun_u(L),fun_v(L)
        u=1.01
        v=0.5
        do t=1,MaxT,1
            do i=1,L,1
                effect=0.0
                do j=1,L
                    effect=effect+neighbour_matrix(i,j)*(u(j)-u(i))
                end do
                fun_u(i)=u(i)+h*(1.0/epsilon*u(i)*(1.0-u(i))*(u(i)-(v(i)+b)/a)+sigma/(2.0*R)*effect+A1*sin(2.0*PI*f*t*h))
                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==10) then
                    write(20,*) t*h,u(i)
                    write(40,*) t*h,v(i)
                end if
            end do
            u=fun_u
            v=fun_v
        end do
        return
    end subroutine fnf

    !建立网状结构
    subroutine neighbour(N,K)
        implicit none
        integer :: N,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 BOM


program main
    use BOM
    implicit none
    open(20,file="data_u.txt")
    open(40,file="data_v.txt")
    call neighbour(L,2*R)
    call random_seed()
    call fnf()
    deallocate(neighbour_matrix)
    close(20)
    close(40)
end program main

figure 4_A

figure 4_B

figure 4_C

module BOM
    implicit none
    integer,parameter :: R=33,MaxT=1500,L=200
    integer,allocatable :: neighbour_matrix(:,:)
    real,parameter :: h=0.02,PI=3.1415926
contains

    subroutine fnf(u0,v0)
        implicit none
        integer :: i,j,t,state
        real :: epsilon=0.04,a=0.84,b=0.07,A1=0.01,f=4.0,sigma=0.132,effect,u(L),v(L),fun_u(L),fun_v(L),U_max(L),U_min(L),u0,v0
        U_max=-10.0
        U_min=10.0
        u=u0
        v=v0
        do t=1,MaxT,1
            do i=1,L,1
                effect=0.0
                do j=1,L
                    effect=effect+neighbour_matrix(i,j)*(u(j)-u(i))
                end do
                fun_u(i)=u(i)+h*(1.0/epsilon*u(i)*(1.0-u(i))*(u(i)-(v(i)+b)/a)+sigma/(2.0*R)*effect+A1*sin(2.0*PI*f*t*h))
                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
            end do
            if(t>=1000) then
                call U_max_min(u,U_max,U_min)
            end if
            u=fun_u
            v=fun_v
        end do
        do j=10,10,1
            if(U_max(j)>1.2) then
                state = 1
            else
                state = 0
            end if
            write(20,*) u0,v0,state
            write(*,*) u0,v0,state
        end do
        return
    end subroutine fnf

    !建立网状结构
    subroutine neighbour(N,K)
        implicit none
        integer :: N,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

    !Umax与Umin
    subroutine U_max_min(u,U_max,U_min)
        implicit none
        real :: U_max(L),U_min(L),u(L)
        integer :: i
        do i=1,L,1
            if(U_max(i)<u(i)) then
                U_max(i)=u(i)
            end if
            if(U_min(i)>u(i)) then
                U_min(i)=u(i)
            end if
        end do
        return
    end subroutine U_max_min

end module BOM


program main
    use BOM
    implicit none
    integer :: a,b
    open(20,file="data_u0_v0.txt")
    call neighbour(L,2*R)
    do a=1,200,1
        do b=1,200,1
            call fnf(a*0.01,b*0.01)
        end do
    end do
    deallocate(neighbour_matrix)
    close(20)
end program main

figure 5

module BOM
    implicit none
    integer,parameter :: R=33,MaxT=50000,L=200
    integer,allocatable :: neighbour_matrix(:,:)
    real,parameter :: h=0.02,PI=3.1415926
contains

    subroutine fnf()
        implicit none
        integer :: i,j,t,num
        real :: epsilon=0.04,a=0.84,b=0.07,A1=3.0,f=4.0,sigma=0.132,effect,u(L),v(L),fun_u(L),fun_v(L),x,y,Z(L),U_max(L),U_min(L)
        U_max=-10.0
        U_min=10.0
        do i=1,L,1
            read(10,*) u(i),v(i)
        end do
        do t=1,MaxT,1
            do i=1,L,1
                effect=0.0
                do j=1,L
                    effect=effect+neighbour_matrix(i,j)*(u(j)-u(i))
                end do
                fun_u(i)=u(i)+h*(1.0/epsilon*u(i)*(1.0-u(i))*(u(i)-(v(i)+b)/a)+sigma/(2.0*R)*effect+A1*sin(2.0*PI*f*t*h))
                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
                call parameter_Z(u,v,i,Z)
                if(mod(t,50)==0) then
                    write(20,*) i,t*h,u(i)
                    write(40,*) i,t*h,Z(i)
                end if
            end do
            if(t>=45000) then
                call U_max_min(u,U_max,U_min)
            end if
            u=fun_u
            v=fun_v
        end do
        do j=1,L,1
            if(U_max(j)<0.2) then
                write(61,*) j,U_max(j),U_min(j)
            else if(U_max(j)>0.2.and.U_max(j)<1.2) then
                write(62,*) j,U_max(j),U_min(j)
            else
                write(63,*) j,U_max(j),U_min(j)
            end if
        end do
        return
    end subroutine fnf

    !局域有序参数Z
    subroutine parameter_Z(u,v,i,Z)
        implicit none
        real :: u(L),v(L),Z(L),Ans_Z_X,Ans_Z_Y
        integer :: i,j
        Ans_Z_X=0.0
        Ans_Z_Y=0.0
        do j=1,L,1
            Ans_Z_X=Ans_Z_X+neighbour_matrix(i,j)*(cos(atan(v(j)/u(j))))
            Ans_Z_Y=Ans_Z_Y+neighbour_matrix(i,j)*(sin(atan(v(j)/u(j))))
        end do
        Z(i)=sqrt((Ans_Z_X/(2.0*R))**2+(Ans_Z_Y/(2.0*R))**2)
        return
    end subroutine parameter_Z

    !Umax与Umin
    subroutine U_max_min(u,U_max,U_min)
        implicit none
        real :: U_max(L),U_min(L),u(L)
        integer :: i
        do i=1,L,1
            if(U_max(i)<u(i)) then
                U_max(i)=u(i)
            end if
            if(U_min(i)>u(i)) then
                U_min(i)=u(i)
            end if
        end do
        return
    end subroutine U_max_min


    !建立网状结构
    subroutine neighbour(N,K)
        implicit none
        integer :: N,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 BOM


program main
    use BOM
    implicit none
    open(10,file="u0_article.txt",access='sequential',position = 'rewind',action = 'read')
    open(20,file="data_u.txt")
    open(40,file="data_z.txt")
    open(61,file="data_u_max_min_blue.txt")
    open(62,file="data_u_max_min_green.txt")
    open(63,file="data_u_max_min_red.txt")
    call neighbour(L,2*R)
    call random_seed()
    call fnf()
    deallocate(neighbour_matrix)
    close(20)
end program main

figure 5_A

figure 5_B

figure 5_C

figure 5_D

figure 6

figure 6_A

module chimera
    implicit none
    integer,parameter :: L=200,MaxT=1500,R=33
    integer,allocatable :: neighbour_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=3.0,f=4.0,sigma=0.132,effect,u(L),v(L),fun_u(L),fun_v(L)
        integer :: i,j,t
        u=0.9
        v=1.0
        do t=1,MaxT,1
            do i=1,L,1
                effect=0
                do j=1,L
                    effect=effect+neighbour_matrix(i,j)*(u(j)-u(i))
                end do
                fun_u(i)=u(i)+h*(1.0/epsilon*u(i)*(1.0-u(i))*(u(i)-(v(i)+b)/a)+sigma/2.0*R*effect+A1*sin(2.0*PI*f*t*h))
                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)
                    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 neighbour(N,K)
        implicit none
        integer :: N,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 chimera


program main
    use chimera
    implicit none
    integer :: i
    open(10,file="neighbour_matrix.txt")
    open(20,file="data_u.txt")
    open(40,file="data_v.txt")
    call neighbour(L,2*R)
    call fnf()
    deallocate(neighbour_matrix)
end program main

figure 6_B

figure 6_C

figure 7

module BOM
    implicit none
    integer,parameter :: R=33,MaxT=1500,L=200
    integer,allocatable :: neighbour_matrix(:,:)
    real,parameter :: h=0.02,PI=3.1415926
contains

    subroutine fnf(A1,f)
        implicit none
        integer :: i,j,t,num=0
        real :: epsilon=0.04,a=0.84,b=0.07,A1,f,sigma=0.132,effect,u(L),v(L),fun_u(L),fun_v(L),x,y,Z(L),U_max(L),U_min(L)
        U_max=-10.0
        U_min=10.0
        do i=1,L,1
            read(10,*) u(i),v(i)
        end do
        do t=1,MaxT,1
            do i=1,L,1
                effect=0.0
                do j=1,L
                    effect=effect+neighbour_matrix(i,j)*(u(j)-u(i))
                end do
                fun_u(i)=u(i)+h*(1.0/epsilon*u(i)*(1.0-u(i))*(u(i)-(v(i)+b)/a)+sigma/(2.0*R)*effect+A1*sin(2.0*PI*f*t*h))
                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
            end do
            if(t>=1000) then
                call U_max_min(u,U_max,U_min)
            end if
            u=fun_u
            v=fun_v
        end do
        do j=1,1,1
            if(U_max(j)<0.2) then
                num = 0
            else if(U_max(j)>0.2.and.U_max(j)<1.2) then
                num = 1
            else
                num = 2
            end if
            write(40,*) A1,f,num
            write(*,*) A1,f,num
        end do
        return
    end subroutine fnf


    !Umax与Umin
    subroutine U_max_min(u,U_max,U_min)
        implicit none
        real :: U_max(L),U_min(L),u(L)
        integer :: i
        do i=1,L,1
            if(U_max(i)<u(i)) then
                U_max(i)=u(i)
            end if
            if(U_min(i)>u(i)) then
                U_min(i)=u(i)
            end if
        end do
        return
    end subroutine U_max_min


    !建立网状结构
    subroutine neighbour(N,K)
        implicit none
        integer :: N,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 BOM


program main
    use BOM
    implicit none
    integer :: a,b
    open(40,file="data_A_f.txt")
    call neighbour(L,2*R)
    call random_seed()
    do a=1,50,1
        do b=1,50,1
            open(10,file="u0_article.txt",access='sequential',position = 'rewind',action = 'read')
            call fnf(a*0.1,b*0.1)
            close(10)
        end do
    end do
    deallocate(neighbour_matrix)
    close(40)
end program main

figure 8

figure 8_随机网络

module BOM
    implicit none
    integer,parameter :: MaxT=50000,L=200
    integer,allocatable :: matrix(:,:)
    real,parameter :: h=0.02,PI=3.1415926
contains

    subroutine fnf()
        implicit none
        integer :: i,j,t,num
        real :: epsilon=0.04,a=0.84,b=0.07,A1=0.01,f=4.0,sigma=0.02,effect,u(L),v(L),fun_u(L),fun_v(L),x,y,Z(L),U_max(L),U_min(L)
        U_max=-10.0
        U_min=10.0
        do i=1,L,1
            call random_number(x)
            call random_number(y)
            u(i)=2.0*x
            v(i)=2.0*y
            write(10,*) u(i),v(i)
        end do
        do t=1,MaxT,1
            do i=1,L,1
                effect=0.0
                do j=1,L
                    effect=effect+matrix(i,j)*(u(j)-u(i))
                end do
                fun_u(i)=u(i)+h*(1.0/epsilon*u(i)*(1.0-u(i))*(u(i)-(v(i)+b)/a)+0.002*effect+A1*sin(2.0*PI*f*t*h))
                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(mod(t,50)==0) then
                    write(20,*) i,t*h,u(i)
                end if
            end do
            if(t>=45000) then
                call U_max_min(u,U_max,U_min)
            end if
            u=fun_u
            v=fun_v
        end do
        do j=1,L,1
            if(U_max(j)<0.2) then
                write(41,*) j,U_max(j),U_min(j)
            else if(U_max(j)>0.2.and.U_max(j)<1.2) then
                write(42,*) j,U_max(j),U_min(j)
            else
                write(43,*) j,U_max(j),U_min(j)
            end if
        end do
        return
    end subroutine fnf

    !Umax与Umin
    subroutine U_max_min(u,U_max,U_min)
        implicit none
        real :: U_max(L),U_min(L),u(L)
        integer :: i
        do i=1,L,1
            if(U_max(i)<u(i)) then
                U_max(i)=u(i)
            end if
            if(U_min(i)>u(i)) then
                U_min(i)=u(i)
            end if
        end do
        return
    end subroutine U_max_min


    !建立网状结构
    subroutine network(N,P)
        implicit none
        integer :: N,i,j,a,b
        real :: x,P
        allocate(matrix(N,N))
        !初始化网络矩阵
        matrix=0
        !赋值节点关系,随机判断
        do i=N,2,-1
            do j=i-1,1,-1
                !产生随机数
                call random_number(x)
                if(x<p) then
                    matrix(j,i)=1
                    matrix(i,j)=1
                end if
            end do
        end do
    end subroutine network

end module BOM


program main
    use BOM
    implicit none
    open(20,file="data_u.txt")
    open(41,file="data_umax_umin_blue.txt")
    open(42,file="data_umax_umin_green.txt")
    open(43,file="data_umax_umin_red.txt")
    call random_seed()
    call network(L,0.4)
    call fnf()
    deallocate(matrix)
    close(20)
    close(41)
    close(42)
    close(43)
end program main

figure 8_小世界网络

module BOM
    implicit none
    integer,parameter :: R=22,MaxT=50000,L=200
    integer,allocatable :: matrix(:,:)
    real,parameter :: h=0.02,PI=3.1415926
contains

    subroutine fnf()
        implicit none
        integer :: i,j,t,num
        real :: epsilon=0.04,a=0.84,b=0.07,A1=0.01,f=4.0,sigma=0.132,effect,u(L),v(L),fun_u(L),fun_v(L),U_max(L),U_min(L)
        U_max=-10.0
        U_min=10.0
        do i=1,L,1
            call random_number(u(i))
            call random_number(v(i))
            u(i)=2.0*u(i)
            v(i)=2.0*v(i)
        end do
        !do i=1,L,1
        !    read(10,*) u(i),v(i)
        !end do
        do t=1,MaxT,1
            do i=1,L,1
                effect=0.0
                do j=1,L
                    effect=effect+matrix(i,j)*(u(j)-u(i))
                end do
                fun_u(i)=u(i)+h*(1.0/epsilon*u(i)*(1.0-u(i))*(u(i)-(v(i)+b)/a)+0.02*effect+A1*sin(2.0*PI*f*t*h))
                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(mod(t,50)==1) then
                    write(20,*) i,t*h,u(i)
                end if
            end do
            if(t>=49000) then
                call U_max_min(u,U_max,U_min)
            end if
            u=fun_u
            v=fun_v
        end do
        do j=1,L,1
            if(U_max(j)<0.2) then
                write(41,*) j,U_max(j),U_min(j)
            else if(U_max(j)>0.2.and.U_max(j)<1.2) then
                write(42,*) j,U_max(j),U_min(j)
            else
                write(43,*) j,U_max(j),U_min(j)
            end if
        end do
        return
    end subroutine fnf

    !Umax与Umin
    subroutine U_max_min(u,U_max,U_min)
        implicit none
        real :: U_max(L),U_min(L),u(L)
        integer :: i
        do i=1,L,1
            if(U_max(i)<u(i)) then
                U_max(i)=u(i)
            end if
            if(U_min(i)>u(i)) then
                U_min(i)=u(i)
            end if
        end do
        return
    end subroutine U_max_min


    !建立网状结构
    subroutine network(N,K,P_con)
        implicit none
        integer :: N,K !单节点K个邻居
        integer :: i,j,a,b,number
        real :: x,P,P_con
        integer :: y
        allocate(matrix(N,N))
        !初始化网络矩阵
        matrix=0
        !赋值节点关系
        !分析可知：邻边矩阵对应|i-j|>=N-K/2或者|i-j|<=K/2
        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
                    matrix(i,j)=1
                    matrix(j,i)=1
                end if
            end do
        end do
        !均匀随机化
        do i=1,N,1
            call random_number(P)
            if(P<=P_con) then
100             call random_number(x)
                y=nint((x-0.51/N)*N+1)!取N的概率比其它少0.02
                number=0
                do a=1,N,1
                    if(matrix(i,a)==1) then
                        number=number+1
                    end if
                end do
                if(number<N-1) then !判断该节点是否可以加边，可以就加，不行就遍历下一个节点
                    if(matrix(i,y)==0.and.i/=y) then
                        matrix(i,y)=1
                        matrix(y,i)=1
                    else
                        goto 100
                    end if
                else
                    cycle
                end if
            end if
        end do
    end subroutine network

end module BOM


program main
    use BOM
    implicit none
    !open(10,file="u0_article.txt",access='sequential',position = 'rewind',action = 'read')
    open(20,file="data_u.txt")
    open(41,file="data_umax_umin_blue.txt")
    open(42,file="data_umax_umin_green.txt")
    open(43,file="data_umax_umin_red.txt")
    call random_seed()
    call network(L,2*R,1.0)
    call fnf()
    deallocate(matrix)
    !close(10)
    close(20)
    close(41)
    close(42)
    close(43)
end program main

figure 8_BA无标度网络

module BOM
    implicit none
    integer,parameter :: MaxT=50000,L=200
    integer,allocatable :: matrix(:,:)
    real,parameter :: h=0.02,PI=3.1415926
contains

    subroutine fnf()
        implicit none
        integer :: i,j,t,num
        real :: epsilon=0.04,a=0.84,b=0.07,A1=0.01,f=4.0,sigma=0.02,effect,u(L),v(L),fun_u(L),fun_v(L),x,y,Z(L),U_max(L),U_min(L)
        U_max=-10.0
        U_min=10.0
        do i=1,L,1
            read(10,*) u(i),v(i)
        end do
        do t=1,MaxT,1
            do i=1,L,1
                effect=0.0
                do j=1,L
                    effect=effect+matrix(i,j)*(u(j)-u(i))
                end do
                fun_u(i)=u(i)+h*(1.0/epsilon*u(i)*(1.0-u(i))*(u(i)-(v(i)+b)/a)+0.02*effect+A1*sin(2.0*PI*f*t*h))
                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(mod(t,50)==0) then
                    write(20,*) i,t*h,u(i)
                end if
            end do
            if(t>=45000) then
                call U_max_min(u,U_max,U_min)
            end if
            u=fun_u
            v=fun_v
        end do
        do j=1,L,1
            if(U_max(j)<0.2) then
                write(41,*) j,U_max(j),U_min(j)
            else if(U_max(j)>0.2.and.U_max(j)<1.2) then
                write(42,*) j,U_max(j),U_min(j)
            else
                write(43,*) j,U_max(j),U_min(j)
            end if
        end do
        return
    end subroutine fnf

    !Umax与Umin
    subroutine U_max_min(u,U_max,U_min)
        implicit none
        real :: U_max(L),U_min(L),u(L)
        integer :: i
        do i=1,L,1
            if(U_max(i)<u(i)) then
                U_max(i)=u(i)
            end if
            if(U_min(i)>u(i)) then
                U_min(i)=u(i)
            end if
        end do
        return
    end subroutine U_max_min


    !建立网状结构
    subroutine network(N,P_ER,M,K)
        implicit none
        integer :: N,K,M !K新节点度
        integer :: i,j,a,b
        real :: sum1,sum2,sum3
        real,allocatable :: sequence_1(:)
        real :: x,P_ER,P
        allocate(matrix(N+M,N+M))
        allocate(sequence_1(N+M-1))
        !初始化网络矩阵
        matrix=0
        !赋值节点关系
        !建立ER随机网络
        do i=N,2,-1
            do j=i-1,1,-1
                !产生随机数
                call random_number(x)
                if(x<P_ER) then
                    matrix(i,j)=1
                    matrix(j,i)=1
                end if
            end do
        end do
        !依次加节点，连边
        do i=1,M,1
            !添加节点
            matrix(:,N+i)=0
            matrix(N+i,:)=0
            !连边
            !sum1旧节点总度数
            sum1=0.0
            sum1=sum(matrix(:,:))
            !sum2旧节点每个节点度数，sequence旧节点的比例情况
            sequence_1=0.0
            do a=1,N+i-1,1
                sum2=0.0
                do b=1,N+i-1,1
                    if(matrix(a,b)==1) then
                        sum2=sum2+1
                    end if
                end do
                sequence_1(a)=sum2/sum1
            end do
            do j=1,K,1
100             call random_number(P) !与旧节点连接概率
                !每个节点度比例
                sum3=0.0 !判断随机数P所在位置
                do a=1,N+M-1,1
                    sum3=sum3+sequence_1(a)
                    if(P<=sum3) then
                        if(matrix(a,N+i)==0) then
                            matrix(a,N+i)=1
                            matrix(N+i,a)=1
                            exit
                        else
                            goto 100
                        end if
                    end if
                end do
            end do
        end do
        deallocate(sequence_1)
    end subroutine network

end module BOM


program main
    use BOM
    implicit none
    open(10,file="u0_article.txt",access='sequential',position = 'rewind',action = 'read')
    open(20,file="data_u.txt")
    open(41,file="data_umax_umin_blue.txt")
    open(42,file="data_umax_umin_green.txt")
    open(43,file="data_umax_umin_red.txt")
    call network(2,1.0,L-2,2)
    call fnf()
    deallocate(matrix)
    close(10)
    close(20)
    close(41)
    close(42)
    close(43)
end program main