Processing math: 58%
Review

Dyadic Brain - A Biological Model for Deliberative Inference

  • The human brain is arguably the most complex information processing system. It operates by acquiring data from the environment, recognizing patterns of events’ occurrence, anticipating their re-occurrence and in turn generating appropriate behavioral responses. Through the lenses of the free-energy principle any self-organizing system that is at equilibrium with its environment must minimize its free energy either by manipulating the environmental sensory input or by manipulating its internal states thus altering the recognition density of the outside stimuli. However, several sets of challenges interfere with the human brain's ability to learn and adapt in such a theoretically optimal fashion. These may include, and are not limited to, functional inconsistencies related to attention and memory processes, the functions of “fast” and “slow” thinking and responding, and the ability of emotional states to generate unintended behavioral outcomes that are less adaptive or inappropriate. This paper will review literature on the subject of how ideal learning viewed from the free-energy principle perspective may be affected by the above mentioned limitations and will suggest a model of information processing that may have developed as a way of overcoming these challenges. This neurobiological model stipulates that a neuronal network is formed in response to environmental input and is paralleled by at least one and possibly multiple networks that activate intrinsically and represent “virtual responses” to a situation that demands a behavioral response. This model accounts for how the brain generates a multiplicity of potential behavioral responses and may “choose” the one that seems most appropriate and also explains the uncanny ability of humans to socialize and collaborate. Implications for understanding humans’ ability to learn from others, deliberate on opposing constructs and access and utilize information outside of individual minds are also discussed.

    Citation: Iliyan Ivanov, Kristin Whiteside. Dyadic Brain - A Biological Model for Deliberative Inference[J]. AIMS Neuroscience, 2017, 4(4): 169-188. doi: 10.3934/Neuroscience.2017.4.169

    Related Papers:

    [1] Jin Li . Linear barycentric rational collocation method to solve plane elasticity problems. Mathematical Biosciences and Engineering, 2023, 20(5): 8337-8357. doi: 10.3934/mbe.2023365
    [2] Wenting Li, Ailing Jiao, Wei Liu, Zhaoying Guo . High-order rational-type solutions of the analogous (3+1)-dimensional Hirota-bilinear-like equation. Mathematical Biosciences and Engineering, 2023, 20(11): 19360-19371. doi: 10.3934/mbe.2023856
    [3] Ishtiaq Ali . Bernstein collocation method for neutral type functional differential equation. Mathematical Biosciences and Engineering, 2021, 18(3): 2764-2774. doi: 10.3934/mbe.2021140
    [4] Andras Balogh, Tamer Oraby . Stochastic games of parental vaccination decision making and bounded rationality. Mathematical Biosciences and Engineering, 2025, 22(2): 355-388. doi: 10.3934/mbe.2025014
    [5] Baojian Hong . Bifurcation analysis and exact solutions for a class of generalized time-space fractional nonlinear Schrödinger equations. Mathematical Biosciences and Engineering, 2023, 20(8): 14377-14394. doi: 10.3934/mbe.2023643
    [6] Bruno Buonomo, Alberto d’Onofrio, Deborah Lacitignola . Rational exemption to vaccination for non-fatal SIS diseases: Globally stable and oscillatory endemicity. Mathematical Biosciences and Engineering, 2010, 7(3): 561-578. doi: 10.3934/mbe.2010.7.561
    [7] Alessia Andò, Simone De Reggi, Davide Liessi, Francesca Scarabel . A pseudospectral method for investigating the stability of linear population models with two physiological structures. Mathematical Biosciences and Engineering, 2023, 20(3): 4493-4515. doi: 10.3934/mbe.2023208
    [8] Max-Olivier Hongler, Roger Filliger, Olivier Gallay . Local versus nonlocal barycentric interactions in 1D agent dynamics. Mathematical Biosciences and Engineering, 2014, 11(2): 303-315. doi: 10.3934/mbe.2014.11.303
    [9] Mario Lefebvre . An optimal control problem without control costs. Mathematical Biosciences and Engineering, 2023, 20(3): 5159-5168. doi: 10.3934/mbe.2023239
    [10] Dimitri Breda, Davide Liessi . A practical approach to computing Lyapunov exponents of renewal and delay equations. Mathematical Biosciences and Engineering, 2024, 21(1): 1249-1269. doi: 10.3934/mbe.2024053
  • The human brain is arguably the most complex information processing system. It operates by acquiring data from the environment, recognizing patterns of events’ occurrence, anticipating their re-occurrence and in turn generating appropriate behavioral responses. Through the lenses of the free-energy principle any self-organizing system that is at equilibrium with its environment must minimize its free energy either by manipulating the environmental sensory input or by manipulating its internal states thus altering the recognition density of the outside stimuli. However, several sets of challenges interfere with the human brain's ability to learn and adapt in such a theoretically optimal fashion. These may include, and are not limited to, functional inconsistencies related to attention and memory processes, the functions of “fast” and “slow” thinking and responding, and the ability of emotional states to generate unintended behavioral outcomes that are less adaptive or inappropriate. This paper will review literature on the subject of how ideal learning viewed from the free-energy principle perspective may be affected by the above mentioned limitations and will suggest a model of information processing that may have developed as a way of overcoming these challenges. This neurobiological model stipulates that a neuronal network is formed in response to environmental input and is paralleled by at least one and possibly multiple networks that activate intrinsically and represent “virtual responses” to a situation that demands a behavioral response. This model accounts for how the brain generates a multiplicity of potential behavioral responses and may “choose” the one that seems most appropriate and also explains the uncanny ability of humans to socialize and collaborate. Implications for understanding humans’ ability to learn from others, deliberate on opposing constructs and access and utilize information outside of individual minds are also discussed.



    The propagation of electromagnetic waves in some materials is usually modeled by the famous Maxwell's equations [4] with various proper medium responses. These significant responses reflect the material's properties, such as the magnetic permeability and electric permittivity with respect to the location and the frequency of the propagating field. When it turns to high intensity radiation situation, not only the medium quantities may depend on the magnitude of the propagating field, but also the response will become nonlinear. In nonlinear optics, one may often focus on the propagation of monochromatic waves, such as continuous high intensity laser beams. In this case, some reasonable assumptions (see [3]) simplify the Maxwell's models to a nonlinear Helmholtz (NLH) equation [5,13]

    ΔE(X)+ω20c2n2(X,|E|)E(X)=f(X),n2(X,|E|)=n20(X)+2n0(X)n2(X)|E|2, (1)

    where E(X) denotes the electric field, X=[x1,,xD] is the spatial coordinate (xD is called the longitudinal direction which also denotes as z as that in Fig. 1 and the rest coordinates are regarded as transverse directions), f(X) is a given function, ω0 is the frequency, c is the speed of light in vacuum, Δ=2x1++2xD is the D-dimensional Laplacian operator, n0(X) and n2(X) are the linear index of refraction and the Kerr coefficient, respectively. And n0,n2 are always assumed to be real which means that all mediums are transparent. The above NLH equation can be used to govern the propagation of linearly polarized, time-harmonic electromagnetic waves in Kerr-type dielectrics, which can produce some important nonlinear optical effects, such as the optical bistability [6,30] and spatial solitons [20].

    Figure 1.  Kerr medium.

    Let the Kerr medium be surrounded by the linear homogeneous medium in which n0=next0,n2=0. We introduce the linear wave number k0=ω0next0/c and denote the normalized quantities

    v(X)=[n0(X)next0]2,ε(X)=2n2(X)n0(X)(next0)2.

    Then, equation (1) can be rewritten as

    ΔE(X)+k20[v(X)+ε(X)|E(X)|2]E(X)=f(X). (2)

    Since the Kerr medium coefficient n2 is discontinuous at the interfaces z=0 and z=Zmax, and n0 may also be discontinuous at these interfaces, it naturally leads to the discontinuities of the coefficients v and ε in equation (2). Thus, for the homogeneous Kerr medium, the coefficients v and ε are piecewise constants as follows

    v={1,z<0,vint ,0zZmax,1,z>Zmax,ε={0,z<0,εint ,0zZmax,0,z>Zmax.

    Remark 1. When the Kerr medium is inhomogeneous which means the whole Kerr medium will be cut into pieces by some other mediums, the linear medium for example. In this case, the discontinuities of the coefficients v and ε in equation (2) will be more complex, we refer [2,3] for the detail.

    When the electric field E and the material coefficients n0 and n2 vary only in one direction z, the model (2) reduces to the 1D nonlinear Helmholtz equation

    d2Edz2+k20(v+ε|E|2)E=f,z(0,Zmax). (3)

    To solve (3) in the interval [0,Zmax], some boundary conditions are needed. It is well-known that, at the interfaces z=0 and z=Zmax, the electronic field E and its first derivative are continuous [3],

    E(0+)=E(0),dEdz(0+)=dEdz(0),E(Z+max)=E(Zmax),dEdz(Z+max)=dEdz(Zmax).

    According to this fact, the so-called two-way boundary conditions which are also used in [3,14,15] are as follows

    (ddz+ik0)E|z=0=2ik0,(ddzik0)E|z=Zmax=0, (4)

    where i=1 is the imaginary unit.

    Remark 2. According to [3], the above boundary condition (4) is developed from the inhomogeneous Sommerfeld type relation and its complete form is

    (ddz+ik0)E|z=0=2ik0E0inc,(ddzik0)E|z=Zmax=2ik0EZmaxinc,

    where E0inc and EZmaxinc are the incident waves which shoot into the computational domain from z= and z=+, respectively. In this paper, we assume that the laser beam only shoots in the Kerr medium from z=, so EZmaxinc=0. And the problem (3) can be rescaled by ˜E=E/E0inc,˜ε=ε|E0inc|2. Without loss of generality, we assume E0inc=1.

    In 2D case, assuming the linear homogeneous medium in R2Ω0 is truncated by a finite domain Ω and Ω0Ω is fully filled with the Kerr medium (see Fig. 2), then (2) is followed by,

    ΔE+k20(v+ε|E|2)E=f,(x,y)Ω:=(ax,bx)×(ay,by),En+ik0E=g,(x,y)Γ:=Γ1Γ2Γ3Γ4, (5)
    Figure 2.  Computational domain in 2D.

    where n denotes the unit outer normal vector to Γ and g=Eincn+ik0Einc, ε0 in ΩΩ0, ε0 in Ω0.

    Many researches are done for the NLH equation. Through transforming it into a phase-amplitude equation, Chen and Mills proposed an approach to obtain the closed form solution of the NLH equation with a single nonlinear layer [7] and multilayered structures [8]. And by using the multidimensional generalization of the nonlinear Schrödinger equation, the exact solution of the NLH equation in special case was also considered in [18]. Moreover, in [11,12], the existence and asymptotic behavior of the real-valued standing wave solution of the NLH equations were analyzed. On the other hand, numerical methods for the NLH equations are also investigated. Fibich and the collaborators studied the NLH equation in [1,2,3,14,15]: In [14], the authors constructed a two-way artificial boundary condition for the NLH equation to ensure that not only the backscattered waves generate no reflection but also the correct value of the incoming wave can be imposed; the NLH equation was also solved by using nonorthogonal expansions in [15]; by coupling with a new technique of the separation of variables, a fourth order finite difference scheme was developed in [1], and the algorithm was also extended to the three-dimensional axially symmetric problem; in [2,3], the authors solved the NLH equations by developing an efficient Newton's iteration method to deal with the strong nonlinearity. In addition, a finite element method which can approximate the discontinuous coefficient problem is constructed in [22,23]. Recently, Wu and Zou proved the existence and uniqueness of the NLH equation, and also analyzed the stability and the error estimate with explicit wave numbers for the finite element approximation in [29]. And the author proposed a robust modified Newton's method in [31].

    There are many difficulties when the NLH equation is approximated by using numerical schemes. Firstly, since the NLH equation is a strong nonlinear problem, we need to search a robust iteration method for solving it. Secondly, similar to the Helmholtz equation, the solution of this problem is highly oscillating with a high number (see [19]). Moreover, it usually contains discontinuous coefficients due to the different propagating mediums. There are also other issues to be solved, such as the strong indefinite linear system generated from this equation and so on. But, in this paper, we mainly focus on the case that it admits highly oscillated solutions. It is well-known that, the FDM is one of the most frequently used numerical methods due to its simple structure. Usually, the FDM is constructed through approximating the derivative terms in the original equation with some difference quotients which are obtained by the Taylor's expansion directly [21,32,33], this may lead the FDM to suffer from some disadvantages, low computation accuracy for example. To simulate the Helmholtz equation with high wave numbers in one dimension, a new finite difference scheme was proposed in [25]. Being different from the classical one directly based on the Taylor expansion, the new finite difference method is constructed by a rearranged formula which can contain more regularity information of the solution, and thus more accurate approximations were achieved. Recently, this kind of schemes were applied to the higher dimensional problems in [16,17,26,27,28]. But all of these problems investigated above are linear problems. In this article, we will extend this idea to solve the NLH equation. Through some iteration methods, the NLH equation will be linearized as a linear one at each iterative step firstly. Several iteration methods are considered, including the classical ones and the error correction iteration method [24] in which the original iteration solution was modified by a residual. Then, based on the above resulted linear problem, the new finite difference scheme is constructed, which is naturally suitable for the problem with discontinuous coefficients to match the different propagating mediums (Kerr and linear mediums).

    The rest of this paper is organized as follows. In the next section, we apply several iteration methods including the error correction one to linearize the NLH equations. Then, in Section 3, after constructing the new finite difference scheme for the 1D linearized equation, we extend the scheme to the 2D problem by using the ADI technique. To test the efficiency of the numerical schemes, some numerical experiments are performed in Section 4. We finally make conclusions in Section 5.

    To solve the nonlinear equation, an iteration method is needed. In this section, we introduce several kinds of iteration methods for solving the NLH equation. For convenience, we rewrite the problem as,

    LE+k20ε|E|2E=f, (6)

    where L:=Δ+k20v denotes the linear operator.

    Frozen-nonlinearity iteration may be the simplest iteration method in which the nonlinear term is frozen as a known quantity. For example, by replacing |E|2 with the previous iteration solution, (6) is linearized by

    LEl+1+k20ε|El|2El+1=f, (7)

    where El donates the lth iteration solution. However, it is well known that the frozen-nonlinearity iteration is only the first order convergent. An alternative iteration method is the Newton's method. Letting

    λ(E,¯E):=LE+k20εE2¯Ef=0, (8)

    then there hold,

    λE(E,¯E)=L+2k20εE¯E,λ¯E(E,¯E)=k20εE2, (9)

    where ¯E is the conjugation of E. Assuming E=s+El(s is small enough) and using the Taylor's expansion, we can approximate (8) as

    λ(El,¯El)+λE(El,¯El)s+λ¯E(El,¯El)ˉs=0. (10)

    Letting s=EElEl+1El and substituting (8) and (9) into (10), we can get the Newton's iteration method,

    LEl+1+2k20ε|El|2El+1+k20ε(El)2¯El+1=f+2k20ε|El|2El. (11)

    Furthermore, a modified Newton's method was proposed in [31] by replacing ¯El+1 with ¯El in (11), that is

    LEl+1+2k20ε|El|2El+1=f+k20ε|El|2El. (12)

    In this paper, we also employ the error correction method in [24] for solving the nonlinear Helmholtz equation. Next, we will show the process. For simplicity, the above three iteration methods can be rewritten as a general formula

    L˜El+1+N˜El+1+M¯˜El+1=˜f, (13)

    where

    N:={k20ε|El|2,in (7),2k20ε|El|2,in (11),2k20ε|El|2,in (12),M:={0,in (7),k20ε(El)2,in (11),0,in (12),˜f:={f,in (7),f+2k20ε|El|2El,in (11),f+k20ε|El|2El,in (12).

    Assuming that μl+1 is the error between ˜El+1 and the exact solution E, i.e., μl+1=E˜El+1, then, subtracting (13) from the original equation (6), we get the error equation as follows

    Lμl+1+k20ε|E|2EN˜El+1M¯˜El+1=f˜f. (14)

    Since

    |E|2E=|˜El+1+μl+1|2(˜El+1+μl+1)=(2|˜El+1|2+|μl+1|2)μl+1+(˜El+1)2¯μl+1+|˜El+1|2˜El+1+2|μl+1|2˜El+1+(μl+1)2¯˜El+1,

    freezing the terms (μl+1)2 and |μl+1|2 at μl in the above equation, (14) yields

    Lμl+1+Pμl+1+Q¯μl+1=fμ, (15)

    where

    P=2k20ε|˜El+1|2+k20ε|μl|2,Q=k20ε(˜El+1)2,fμ=f˜f(2k20ε|μl|2+k20ε|˜El+1|2N)˜El+1[k20ε(μl)2M]¯˜El+1.

    Furthermore, for the error on the boundary, there holds

    μl+1n+ik0μl+1=0. (16)

    Obviously, by solving (15), the original solution ˜El+1 could be modified by a more accurate one El+1=˜El+1+μl+1 which is expected to be useful in decreasing the iteration number. And the iterative procedure with the error correction can be concluded in the following Algorithm 1.

    Algorithm 1
    Step 1: Give a initial guess E0 and set μ0=0;
    Step 2: Solve equation (13) to obtain ˜E1;
    Step 3: Solve equation (15) to obtain μ1;
    Step 4: Set E1=˜E1+μ1, if E1E02E02<δ(threshold of the iteration), stop; else, set E0=E1,μ0=μ1 and go back to Step 2.

     | Show Table
    DownLoad: CSV

    After applying the iteration methods introduced in the above section, the NLH equation is linearized to linear problems at each iterative step. To implement these methods, a spatial discretization scheme is needed. From the stability analysis in [29], we know that the solution of the NLH equation satisfying ||E||Hjc0kj10(j=0,1,2) with c0 being a positive constant and Hj being the classical norm in the Sobolev space. This implies the oscillation of the field E when k0 is large. In this section, we will develop a family of efficient finite difference schemes for solving (13) and (15) with corresponding boundary conditions, respectively. Since these two equations have the same form, by taking (13) for an example, we will exhibit the construction of the new finite difference scheme.

    Recalling the iteration formula for the 1D NLH equation

    d2˜El+1dz2+(k20v+N)˜El+1+M¯˜El+1=˜f,z(0,Zmax), (17)
    (ddz+ik0)˜El+1|z=0=2ik0,(ddzik0)˜El+1|z=Zmax=0. (18)

    When the frozen-nonlinearity iteration (7) and the modified Newton's method (12) are used, i.e., M=0, we can rewrite (17) as

    d2˜El+1dz2=τ˜El+1+˜f, (19)

    where τ=(k20v+N). (19) has the same form as the linear Helmholtz equation, so we can extend the idea in [25,27] to solve it.

    Since the parameters τ and ˜f may be discontinuous at the interfaces of the Kerr medium and the linear medium, we divide the computational domain by a mesh in which these interfaces are specially selected as mesh points. Assume that {zm}(1mN) with z1=0,zN=Zmax are the mesh points. And for any interior point zm, setting hm=zmzm1,h+m=zm+1zm and

    τm={τm,zm[zm1,zm],τ+m,zm[zm,zm+1],˜fm={˜fm,zm[zm1,zm],˜f+m,zm[zm,zm+1].

    Taylor's expansion tells that

    ˜El+1m+1=˜El+1m++(h+m)2k(2k)!(˜El+1m)(2k)+(h+m)2k+1(2k+1)!(˜El+1m)(2k+1)+, (20)
    ˜El+1m1=˜El+1m++(hm)2k(2k)!(˜El+1m)(2k)(hm)2k+1(2k+1)!(˜El+1m)(2k+1)+. (21)

    And according to (19), there hold

    (˜El+1m)(2k)=τnm˜El+1m+ks=1τksm˜f(2s2)m, (22)
    (˜El+1m)(2k+1)=τnm(˜El+1m)(1)+ks=1τksm˜f(2s1)m. (23)

    Substituting (22)-(23) into (20) and (21), it yields

    ˜El+1m+1=G+m˜El+1m+H+m(˜El+1m)(1)++k=1[L+k;m(˜f+m)(2k2)+X+k;m(˜f+m)(2k1)], (24)
    ˜El+1m1=Gm˜El+1mHm(˜El+1m)(1)++k=1[Lk;m(˜fm)(2k2)Xk;m(˜fm)(2k1)], (25)

    where

    G±m=G(τ±m,h±m),H±m=H(τ±m,h±m),L±k;m=Lk(τ±m,h±m),X±k;m=Xk(τ±m,h±m),

    with

    G(ρ,h):=ehρ+ehρ2,H(ρ,h):=1ρehρehρ2,Lk(ρ,h):=1ρk[ehρ+ehρ2k1s=0(hρ)2s(2s)!],Xk(ρ,h):=1ρk+1/2[ehρehρ2k1s=0(hρ)2s+1(2s+1)!].

    Then, eliminating (˜El+1m)(1) in (24) and (25), we get the scheme for the interior point as follows

    H+m˜El+1m1(H+mGm+HmG+m)˜El+1m+Hm˜El+1m+1=+Hm+k=1[L+k;m(˜f+m)(2k2)+X+k;m(˜f+m)(2k1)]+H+m+k=1[Lk;m(˜fm)(2k2)Xk;m(˜fm)(2k1)]. (26)

    For the boundary points z=0,z=Zmax, directly letting m=1 and m=N in (24) and (25), respectively, we have

    ˜El+12=G+1˜El+11+H+1(˜El+11)(1)++k=1[L+k;1(˜f+1)(2k2)+X+k;1(˜f+1)(2k1)],˜El+1N1=GN˜El+1NHN(˜El+1N)(1)++k=1[Lk;N(˜fN)(2k2)Xk;N(˜fN)(2k1)].

    Then, substituting the boundary condition (18) into the above formulas, we have the numerical schemes for z=0 and z=Zmax as follows

    (G+1ik0H+1)˜El+11+˜El+12=2ik0H+1++k=1[L+k;1(˜f+1)(2k2)+X+k;1(˜f+1)(2k1)], (27)
    ˜El+1N1(GNik0HN)˜El+1N=+k=1[Lk;N(˜fN)(2k2)Xk;N(˜fN)(2k1)]. (28)

    Obviously, taking different k in (26)-(28) could obtain different finite difference schemes. For example, letting k=1, the new finite difference schemes are

    H+m˜El+1m1(H+mGm+HmG+m)˜El+1m+Hm˜El+1m+1=+Hm[L+1;m˜f+m+2X+1;mhm+h+m(˜f+m+1˜f+m1)]+H+m[L1;m˜fm2X1;mhm+h+m(˜fm+1˜fm1)],m=2,3,,N1,(G+1ik0H+1)˜El+11+˜El+12=2ik0H+1+L+1;1˜f+1+X+1;1h+1(˜f+2˜f+1), (29)
    ˜El+1N1(GNik0HN)˜El+1N=L1;N˜fNX1;NhN(˜fN˜fN1).

    Remark 3. In fact, a more accurate numerical scheme could be developed. For example, when the frozen-nonlinearity iteration method is used, (19) is the equation with variable coefficient, it yields,

    d2˜El+1dz2=τ(z)˜El+1+˜f,

    where τ(z)=k20[v+εEl(z)¯El(z)].

    In this case, like (22)-(23), we can get more precise formulas

    (˜El+1m)(2n)=τnm˜El+1m+τn1m˜fm,(˜El+1m)(2n+1)=τnm(˜El+1m)(1)+(C01+C23++C2n22n1)τn1mτ(1)m˜El+1m+τn1m˜f(1)m,

    where C2n22n1=(2n1)!/(2n2)!.

    Then, according to the Taylor's series, we get

    ˜El+1m=A+m˜Em+B+m˜E(1)m+C+m˜fm+D+m˜f(1)m,˜El1m=Am˜Em+Bm˜E(1)m+Cm˜fm+Dm˜f(1)m,

    where

    A±m=A(τ±m,h±m),B±m=B(τ±m,h±m),C±m=C(τ±m,h±m),D±m=D(τ±m,h±m),A(ν,h)=12[eνh+eνh]+ν(1)8{1(ν)3[eνheνh2νh]}+ν(1)8{h(ν)2[eνh+eνh2]+h2ν[eνheνh]},B(ν,h)=12ν[eνheνh],C(ν,h)=12ν[eνh+eνh],D(ν,h)=12νν[eνheνh].

    Thus, eliminating the terms ˜E(1)m, we can obtain

    1Bm˜Em1(A+mB+mAmBm)˜Em+1B+m˜Em+1=(C+mB+mCmBm)˜fm+(D+mB+mDmBm)˜f(1)m. (30)

    Finally, by approximating τ(1)m,˜f(1)m with difference quotients in (30), we can get a more accurate scheme. Moreover, the corresponding schemes for boundary points can also be constructed easily like (27)-(28).

    When the Newton's iteration method (11) is considered, (17) needs to be separated into the real and imaginary parts due to the existence of ¯˜El+1, which is followed by

    d2Rdz2=ˆRR+ˆII+˜fR, (31)
    (dRdzk0I)|z=0=0,(dRdz+k0I)|z=Zmax=0, (32)
    d2Idz2=II+RR+˜fI, (33)
    (dIdz+k0R)|z=0=2k0,(dIdzk0R)|z=Zmax=0, (34)

    where

    R=real(˜El+1),I=imag(˜El+1),˜fR=real(˜f),˜fI=imag(˜f),ˆR=real(k20v+N+M),ˆI=imag(Mk20vN),I=real(k20v+NM),R=imag(M+k20v+N).

    Taking the real part equation (31) for an example, following the same process for (19), we have, at any interior point zm, that

    R(2k)m=ˆRkmRm+ˆImks=1ˆRksmI(2s2)m+ks=1ˆRksm˜f(2s2)R;m, (35)
    R(2k+1)m=ˆRkmR(1)m+ˆImks=1ˆRksmI(2s1)m+ks=1ˆRksm˜f(2s1)R;m, (36)

    and

    Rm+1=G+mRm+H+mR(1)m++k=1[L+k;m(ˆI+mIm+˜f+R;m)(2k2)+X+k;m(ˆI+mIm+˜f+R;m)(2k1)], (37)
    Rm1=GmRmHmR(1)m++k=1[Lk;m(ˆImIm+˜fR;m)(2k2)Xk;m(ˆImIm+˜fR;m)(2k1)], (38)

    where

    G±m=G(ˆR±m,h±m),H±m=H(ˆR±m,h±m),L±k;m=Lk(ˆR±m,h±m),X±k;m=Xk(ˆR±m,h±m).

    Eliminating R(1)m in (37) and (38), we get

    H+mRm1(H+mG+HmG+)Rm+HmRm+1=+k=1(HmF+k+H+mFk), (39)

    where

    F+k=ˆI+m(L+k;mI(2k2)m+X+k;mI(2k1)m)+L+k;m(˜f+R;m)(2k2)+X+k;m(˜f+R;m)(2k1)Fk=ˆIm(Lk;mI(2k2)mXk;mI(2k1)m)+Lk;m(˜fR;m)(2k2)X+k;m(˜fR;m)(2k1).

    Similarly, setting m=1 and m=N in (37) and (38), respectively, and according to the boundary condition (32), we get

    G+1R1+R2k0H+1I1=+k=1[ˆI+1(L+k;1I(2k2)1+X+k;1I(2k1)1)+L+k;1(˜f+R;1)(2k2)+X+k;1(˜f+R;1)(2k1)], (40)
    RN1GNRNk0HNIN=+k=1[ˆIN(Lk;NI(2k2)NXk;NI(2k1)N)+Lk;N(˜fR;N)(2k2)X+k;N(˜fR;N)(2k1)]. (41)

    Obviously, by retaining different terms in the right hand side of (39)-(41), we can also get a series of finite difference schemes for the real part equation (31)-(32). For example, taking k=1 in (39), it yields

    H+mRm1(H+mG+HmG+)Rm+HmRm+1=(HmˆI+mL+1;m+H+mˆImL1;m)Im+(HmˆI+mX+1;mH+mˆImX1;m)I(1)m+HmL+1;m(˜f+R;m)+H+mL1;m(˜fR;m)HmX+1;m(˜f+R;m)(1)H+mX1;m(˜fR;m)(1).

    By approximating ϕ(1)m(ϕ=I,˜fR) with ϕm+1ϕm1h+m+hm in the above formula, we can get a finite difference scheme for (31) at the interior point zm(m=2,,N1),

    A[R,I]=BFR, (42)

    where

    A=(A1,A2,A3,A4,A5,A6),R=(Rm1,Rm,Rm+1),I=(Im1,Im,Im+1),B=(B1,B2,B3,B4,B5,B6),FR=(˜f+R;m1,˜f+R;m,˜f+R;m+1,˜fR;m1,˜fR;m,˜fR;m+1),

    with

    A1=H+m,A2=H+mGmHmG+m,A3=Hm,A4=A6=HmˆI+mX+1;mH+mˆImX1;mh+m+hm,A5=HmˆI+mL+1;mH+mˆImL1;m,B1=B3=HmX+1;mh+m+hm,B2=HmL+1;m,B4=B6=H+mX1;mh+m+hm,B5=H+mL1;m.

    Similarly, letting k=1 in (40)-(41), we have

    G+1R1+R2k0H+1I1=ˆI+1(L+1;1I1+X+1;1I(1)1)+L+1;1˜f+R;1+X+1;1(˜f+R;1)(1),RN1GNRNk0HNIN=ˆIN(L1;NINX1;NI(1)N)+L1;N˜fR;NX+1;N(˜fR;N)(1).

    By using the boundary condition (34) and approximating (˜f+R,1)(1) and (˜fR,N)(1) with ˜f+R,2˜f+R,1h+1 and ˜fR,N˜fR,N1hN, respectively, we get a finite difference scheme at the boundary points z=0 and z=Zmax,

    (G+1k0ˆI+1X1;1)R1+R2(k0H+1+L+1;1ˆI+1)I1=(L+1;1X+1;1h+1)˜f+R,1+X+1;1h+1˜f+R,2+2k0ˆI+1X1;1, (43)
    RN1(GNk0X1;NˆIN)RN(k0HN+L1;NˆIN)IN=(L1;NX1;NhN)˜fR;N+X1;NhN˜fR,N1. (44)

    For the imaginary part (33)-(34), the same procedure can be applied to develop the new finite difference scheme at any interior point zm(m=2,,N1),

    A[I,R]=BFI, (45)

    where

    FI=(˜f+I;m1,˜f+I;m,˜f+I;m+1,˜fI;m1,˜fI;m,˜fI;m+1),A1=T+m,A2=T+mSmTmS+m,A3=Tm,A4=A6=TmR+mY+1;mT+mRmY1;mh+m+hm,A5=TmR+mM+1;mT+mRmM1;m,B1=B3=TmY+1;mh+m+hm,B2=TmM+1;m,B4=B6=T+mY1;mh+m+hm,B5=T+mM1;m,

    with

    S±m=G(I±m,h±m),T±m=H(I±m,h±m),M±1;m=L1(I±m,h±m),Y±1;m=X1(I±m,h±m).

    And the schemes for boundary points z=0 and z=Zmax are

    (S+1+k0R+1Y+1;1)I1+I2(M+1;1R+1k0T+1)R1=(M+1;1Y+1;1h+1)˜f+I,1+Y+1;1h+1˜f+I,2+2k0T+1, (46)
    IN1(SN+k0Y1;NRN)IN(M1;NRNk0TN)RN=(M1;NY1;NhN)˜fI;N+Y1;NhN˜fI,N1. (47)

    Obviously, (42)-(47) constitute a finite difference scheme for the system of equations (31)-(34). According to the above process, it can be found that, through translating the high order terms E(2n) and E(2n+1) into the lower ones, much more terms in the Taylor's expansion could be included in the new finite difference scheme. Therefore, the new scheme is expected to achieve much better computational accuracy. Moreover, the case of discontinuous coefficients is considered fully and naturally in these schemes.

    In this section, we extend the new finite difference scheme to the 2D problem by applying the ADI method [9,10]. Similar to (31)-(34), the 2D equation (5) need to be divided into real and imaginary parts like

    ΔR=ˆRR+ˆII+˜fR,(x,y)Ω, (48)
    Ryk0I=g1R,(x,y)Γ1,Ryk0I=g2R,(x,y)Γ2,Rxk0I=g3R,(x,y)Γ3,Rxk0I=g4R,(x,y)Γ4, (49)
    ΔI=II+RR+˜fI,(x,y)Ω, (50)
    Iy+k0R=g1I,(x,y)Γ1,Iy+k0R=g2I,(x,y)Γ2,Ix+k0R=g3I,(x,y)Γ3,Ix+k0R=g4I,(x,y)Γ4, (51)

    where glR=real(gl),glI=imag(gl)(l=1,2,3,4).

    According to (5) and Fig. 2, the parameters ˆR,ˆI,I,R are also discontinuous at the interface Ω0. So, we divide the computational domain Ω into Nx×Ny parts such that Ω0 (including its four vertexes) is overlapped with some mesh points. And let (xm,yn)(1mNx,1nNy) be the mesh points and hm=xmxm1,h+m=xm+1xm,kn=ynyn1,k+n=yn+1yn be the mesh sizes.

    It is well-known that the ADI method is used to simulate a high-dimensional problem by solving a series of one-dimensional problems. Based on this, by directly separating the real part equation (48) into two 1D equations in x,y directions at the interior point (xm,yn), we have

    2Rm,nx2=γxˆRm,nRm,n+γxˆIm,nIm,n+˜fxR;m,n, (52)
    2Rm,ny2=γyˆRm,nRm,n+γyˆIm,nIm,n+˜fyR;m,n, (53)

    where γx+γy=1,˜fxR+˜fyR=˜fR.

    Similar to (42), the new finite difference schemes for (52) and (53) can be directly got as follows

    H+xRm1,n(H+xGx+HxG+x)Rm,n+HxRm+1,n(HxγxˆI+m,nL+x,1+H+xγxˆIm,nLx,1)Im,nHxγxˆI+m,nX+x,1+H+xγxˆIm,nXx,1h+m+hm(Im+1,nIm1,n)=(HxL+x,1+H+xLx,1)˜fxR;m,n,H+yRm,n1(H+yGy+HyG+y)Rm,n+HyRm,n+1(HyγyˆI+m,nL+y,1+H+yγyˆIm,nLy,1)Im,n (54)
    HyγyˆI+m,nX+y,1+H+yγyˆIm,nXy,1k+n+kn(Im,n+1Im,n1)=(HyL+y;1+H+yLy;1)˜fyR;m,n, (55)

    where

    G±x=G(γxˆR±m,n,h±m),H±x=H(γxˆR±m,n,h±m),L±x,1=L1(γxˆR±m,n,h±m),X±x,1=X1(γxˆR±m,n,h±m),G±y=G(γyˆR±m,n,k±n),H±y=H(γyˆR±m,n,k±n),L±y;1=L1(γyˆR±m,n,k±n),X±y;1=X1(γyˆR±m,n,k±n).

    Combining (54) and (55), we get the new finite difference scheme for (48) at the interior point (xm,yn)(2mNx1,2nNy1),

    A[R,I]=˜fR;m,n, (56)

    where

    A=(A1,A2,A3,A4,A5,A6,A7,A8,A9,A10),R=(Rm,n1,Rm1,n,Rm,n,Rm+1,n,Rm,n+1),I=(Im,n1,Im1,n,Im,n,Im+1,n,Im,n+1),

    with

    A1=H+yHyL+y;1+H+yLy;1,A2=H+xHxL+x,1+H+xLx,1,A3=H+xGx+HxG+xHxL+x,1+H+xLx,1H+yGy+HyG+yHyL+y;1+H+yLy;1,A4=HxHxL+x,1+H+xLx,1,A5=HyHyL+y;1+H+yLy;1,A6=HyγyˆI+m,nX+y,1+H+yγyˆIm,nXy,1(k+n+kn)(HyL+y;1+H+yLy;1),A7=HxγxˆI+m,nX+x,1+H+xγxˆIm,nXx,1(h+m+hm)(HxL+x,1+H+xLx,1),A8=HxγxˆI+m,nL+x,1+H+xγxˆIm,nLx,1HxL+x,1+H+xLx,1HyγyˆI+m,nL+y,1+H+yγyˆIm,nLy,1HyL+y;1+H+yLy;1,A9=HxγxˆI+m,nX+x,1+H+xγxˆIm,nXx,1(h+m+hm)(HxL+x,1+H+xLx,1),A10=HyγyˆI+m,nX+y,1+H+yγyˆIm,nXy,1(k+n+kn)(HyL+y;1+H+yLy;1).

    Similar to the interior points, the new finite difference scheme for each boundary point is also constructed by developing two schemes in x and y directions, respectively. For example, when (xm,yn)Γ1{Γ1Γ3,Γ1Γ4}, the scheme (54) with n=1 can be seen as the new finite difference scheme in x direction, and in y direction, according to (53), we get

    Rm,2=G+yRm,1+H+yR(1)m,1y+L+y,1˜fyR;m,1+γyˆI+m,1(L+y;1Im,1+X+y,1I(1)m,1y).

    Then, substituting the corresponding boundary condition into the above formula, we have

    (G+y+k0γyˆI+m,1X+y,1)Rm,1+Rm,2+(k0H+yγyˆI+m,1L+y,1)Im,1=L+y;1˜fyR;m,1H+yg1R;m,1γyˆI+m,1X+y;1g1I;m,1. (57)

    So, combining (54)(n=1) and (57), the new finite difference scheme for the boundary points on Γ1{Γ1Γ3,Γ1Γ4} is

    AbRb=Fb, (58)

    where

    Rb=(Rm1,1,Rm,1,Rm+1,1,Im1,1,Im,1,Im+1,1,Rm,2),Ab=(A1,A2,A3,A4,A5,A6,A7),A1=H+xHxL+x,1+H+xLx,1,A3=HxHxL+x,1+H+xLx,1,A2=H+xGx+HxG+xHxL+x,1+H+xLx,1G+y+k0γyˆI+m,1X+y,1L+y;1,A4=HxγxˆI+m,1X+x,1+H+xγxˆIm,1Xx,1(h+m+hm)(HxL+x,1+H+xLx,1),A5=k0H+yγyˆI+m,1L+y,1L+y;1HxγxˆI+m,1L+x,1+H+xγxˆIm,1Lx,1HxL+x,1+H+xLx,1,A6=HxγxˆI+m,1X+x,1+H+xγxˆIm,1Xx,1(h+m+hm)(HxL+x,1+H+xLx,1),A7=1L+y;1,Fb=˜fR;m,1H+yL+y;1g1R;m,1γyˆI+m,1X+y;1L+y;1g1I;m,1.

    Similar to (57), on other three boundaries (excluding the vertexes), we also have: for the boundary points on Γ2{Γ2Γ3,Γ2Γ4} (y direction):

    Rm,Ny1(Gy+k0γyˆIm,NyXy;1)Rm,Ny+(k0HyγyˆIm,NyLy;1)Im,Ny=Ly;1˜fyR;m,NyHyg2R;m,NyγyˆIm,NyXy;1g2I;m,Ny, (59)

    for the boundary points on Γ3{Γ3Γ1,Γ3Γ2}(x direction):

    (G+x+k0γxˆI+1,nX+x;1)R1,n+R2,n+(k0H+xγxˆI+1,nL+x,1)I1,n=L+x;1˜fxR;1,nH+xg3R;1,nγxˆI+1,nX+x;1g3I;1,n, (60)

    for the boundary points on Γ4{Γ4Γ1,Γ4Γ2}(x direction):

    RNx1,n(Gx+k0γxˆINx,nXx;1)RNx,n+(k0HxγxˆINx,nLx,1)INx,n=Lx,1˜fxR;Nx,nHxg4R;Nx,nγxˆINx,nXx;1g4I;Nx,n. (61)

    Thus, applying the same process, the new finite difference schemes for these boundary points can also be written as (58) with different Ab,Rb and Fb. The details are: for the points on Γ2{Γ2Γ3,Γ2Γ4}, there hold

    Rb=(Rm1,Ny,Rm,Ny,Rm+1,Ny,Im1,Ny,Im,Ny,Im+1,Ny,Rm,Ny1),A1=H+xHxL+x,1+H+xLx,1,A3=HxHxL+x,1+H+xLx,1,A2=H+xGx+HxG+xHxL+x,1+H+xLx,1Gy+k0γyˆIm,NyXy;1Ly;1,A4=HxγxˆI+m,NyX+x,1+H+xγxˆIm,NyXx,1(h+m+hm)(HxL+x,1+H+xLx,1),A5=k0HyγyˆIm,NyLy;1Ly;1HxγxˆI+m,NyL+x,1+H+xγxˆIm,NyLx,1HxL+x,1+H+xLx,1,A6=HxγxˆI+m,NyX+x,1+H+xγxˆIm,NyXx,1(h+m+hm)(HxL+x,1+H+xLx,1),A7=1Ly;1,Fb=˜fR;m,NyHyg2R;m,NyLy;1γyˆIm,NyXy;1g2I;m,NyLy;1,

    and for the points on Γ3{Γ3Γ1,Γ3Γ2}, there hold

    Rb=(R1,n1,R1,n,R1,n+1,I1,n1,I1,n,I1,n+1,R2,n),A1=H+yHyL+y;1+H+yLy;1,A3=HyHyL+y;1+H+yLy;1,A2=G+x+k0γxˆI+1,nX+x;1L+x;1H+yGy+HyG+yHyL+y;1+H+yLy;1,A4=HyγyˆI+1,nX+y,1+H+yγyˆI1,nXy,1(k+n+kn)(HyL+y;1+H+yLy;1),A5=k0H+xγxˆI+1,nL+x,1L+x;1HyγyˆI+1,nL+y,1+H+yγyˆI1,nLy,1HyL+y;1+H+yLy;1,A6=HyγyˆI+1,nX+y,1+H+yγyˆI1,nXy,1(k+n+kn)(HyL+y;1+H+yLy;1),A7=1L+x;1,Fb=˜fR;1,nH+xg3R;1,nL+x;1γxˆI+1,nX+x;1g3I;1,nL+x;1,

    and for the points on \Gamma_4\backslash\{\Gamma_4\cap\Gamma_1,\Gamma_4\cap\Gamma_2\} , there hold

    \begin{align*} &\mathbf{R_b} = (R_{N_x,n-1},R_{N_x,n},R_{N_x,n+1},I_{N_x,n-1},I_{N_x,n},I_{N_x,n+1},R_{N_x-1,n}),\\ &A_1 = \frac{H_{y}^{+}}{H_{y}^{-} L_{y ; 1}^{+}+H_{y}^{+} L_{y ; 1}^{-}} ,\quad A_3 = \frac{H_{y}^{-}}{H_{y}^{-} L_{y ; 1}^{+}+H_{y}^{+} L_{y ; 1}^{-}} ,\\ &A_2 = -\frac{G_{x}^{-}+k_{0} \gamma_{x} \widehat{\mathcal{I}}_{N_{x}, n} X_{x;1}^{-}}{L_{x, 1}^{-}} -\frac{H_{y}^{+} G_{y}^{-}+H_{y}^{-} G_{y}^{+}}{H_{y}^{-} L_{y ; 1}^{+}+H_{y}^{+} L_{y ; 1}^{-}} , \\ &A_4 = \frac{H_{y}^{-} \gamma_{y} \widehat{\mathcal{I}}_{N_x, n}^{+} X_{y, 1}^{+}+H_{y}^{+} \gamma_{y} \widehat{\mathcal{I}}_{N_x, n}^{-} X_{y, 1}^{-}}{(k_{n}^{+}+k_{n}^{-})(H_{y}^{-} L_{y ; 1}^{+}+H_{y}^{+} L_{y ; 1}^{-})}, \end{align*}
    \begin{align*} &A_5 = \frac{k_{0}H_{x}^{-}-\gamma_{x} \widehat{\mathcal{I}}_{N_{x}, n} L_{x, 1}^{-}}{L_{x, 1}^{-}} -\frac{H_{y}^{-} \gamma_{y} \widehat{\mathcal{I}}_{N_x, n}^{+} L_{y, 1}^{+}+H_{y}^{+} \gamma_{y} \widehat{\mathcal{I}}_{N_x, n}^{-} L_{y, 1}^{-}}{H_{y}^{-} L_{y ; 1}^{+}+H_{y}^{+} L_{y ; 1}^{-}},\\ &A_6 = -\frac{H_{y}^{-} \gamma_{y} \widehat{\mathcal{I}}_{N_x, n}^{+} X_{y, 1}^{+}+H_{y}^{+} \gamma_{y} \widehat{\mathcal{I}}_{N_x, n}^{-} X_{y, 1}^{-}}{(k_{n}^{+}+k_{n}^{-})(H_{y}^{-} L_{y ; 1}^{+}+H_{y}^{+} L_{y ; 1}^{-})}, \\ &A_7 = \frac{1}{L_{x, 1}^{-}},\quad F_b = \widetilde{f}_{R ; N_{x}, n}-\frac{H_{x}^{-} g_{4 R ; N_{x}, n}}{L_{x, 1}^{-}} -\frac{\gamma_{x} \widehat{\mathcal{I}}_{N_{x}, n}^{-} X_{x ; 1}^{-} g_{4 I ; N_{x}, n}}{L_{x, 1}^{-}}. \end{align*}

    According to (57), (59)-(61), the new finite schemes for four vertexes can also be obtained. For example setting m = 1 and n = 1 in (57) and (60), respectively, and combining them together, we get the finite difference scheme at the vertex \Gamma_1\cap\Gamma_3 ,

    \begin{align} \mathbf{A_v\cdot R_v} = F_v, \end{align} (62)

    where

    \begin{align*} &\mathbf{A_v} = (A_1,A_2,A_3,A_4),\quad \mathbf{R_v} = (R_{1,1},R_{2,1},R_{1,2},I_{1,1}),\\ &A_1 = -\frac{G_{y}^{+}+k_{0} \gamma_{y} \widehat{\mathcal{I}}_{1, 1}^{+} X_{y, 1}^{+}}{L_{y ; 1}^{+} } -\frac{G_{x}^{+}+k_{0} \gamma_{x} \widehat{\mathcal{I}}_{1, 1}^{+} X_{x;1}^{+}}{L_{x ; 1}^{+} },\\ &A_2 = \frac{1}{L_{x ; 1}^{+} },\quad A_3 = \frac{1}{L_{y ; 1}^{+} },\quad A_4 = \frac{k_{0} H_{y}^{+}-\gamma_{y} \widehat{\mathcal{I}}_{1, 1}^{+} L_{y, 1}^{+}}{L_{y ; 1}^{+} } +\frac{k_{0} H_{x}^{+}-\gamma_{x} \widehat{\mathcal{I}}_{1, 1}^{+} L_{x, 1}^{+}}{L_{x ; 1}^{+} },\\ &F_v = \widetilde{f}_{R ; 1, 1}-\frac{H_{y}^{+} }{L_{y ; 1}^{+} }g_{1 R ; 1, 1} -\frac{\gamma_{y} \widehat{\mathcal{I}}_{1, 1}^{+} X_{y ; 1}^{+}}{L_{y ; 1}^{+} } g_{1 I ; 1, 1} -\frac{H_{x}^{+} g_{3 R ; 1, 1}}{L_{x ; 1}^{+} } -\frac{\gamma_{x} \widehat{\mathcal{I}}_{1, 1}^{+} X_{x ; 1}^{+} g_{3 I ; 1, 1}}{L_{x ; 1}^{+} }. \end{align*}

    Similarly, for the rest three vertexes, their new finite difference schemes can be also concluded in (62) with different \mathbf{A_v,R_v} and F_v : for the vertex \Gamma_1\cap\Gamma_4 , there hold

    \begin{align*} &\mathbf{R_v} = (R_{N_x-1,1},R_{N_x,1},R_{N_x,2},I_{N_x,1}),\\ &A_1 = \frac{1}{L_{x, 1}^{-} },\quad A_3 = \frac{1}{L_{y ; 1}^{+} },\\ &A_2 = -\frac{G_{y}^{+}+k_{0} \gamma_{y} \widehat{\mathcal{I}}_{N_x, 1}^{+} X_{y, 1}^{+}}{L_{y ; 1}^{+} } -\frac{G_{x}^{-}+k_{0} \gamma_{x} \widehat{\mathcal{I}}_{N_{x}, 1} X_{x;1}^{-}}{L_{x, 1}^{-} } ,\\ &A_4 = \frac{k_{0} H_{y}^{+}-\gamma_{y} \widehat{\mathcal{I}}_{N_x, 1}^{+} L_{y, 1}^{+}}{L_{y ; 1}^{+} } +\frac{k_{0}H_{x}^{-}-\gamma_{x} \widehat{\mathcal{I}}_{N_{x}, 1} L_{x, 1}^{-}}{L_{x, 1}^{-} } ,\\ &F_v = \widetilde{f}_{R ; N_x, 1}-\frac{H_{y}^{+} }{L_{y ; 1}^{+} }g_{1 R ; N_x, 1} -\frac{\gamma_{y} \widehat{\mathcal{I}}_{N_x, 1}^{+} X_{y ; 1}^{+}}{L_{y ; 1}^{+} } g_{1 I ; N_x, 1}\\ &-\frac{H_{x}^{-} g_{4 R ; N_{x}, 1}}{L_{x, 1}^{-} } -\frac{\gamma_{x} \widehat{\mathcal{I}}_{N_{x}, 1}^{-} X_{x ; 1}^{-} g_{4 I ; N_{x}, 1}}{L_{x, 1}^{-} }, \end{align*}

    and for the vertex \Gamma_2\cap\Gamma_3 , there hold

    \begin{align*} &\mathbf{R_v} = (R_{1,N_y-1},R_{1,N_y},R_{2,N_y},I_{1,N_y}),\\ &A_1 = \frac{1}{L_{y ; 1}^{-} },\quad A_3 = \frac{1}{L_{x ; 1}^{+} }, \end{align*}
    \begin{align*} &A_2 = -\frac{G_{y}^{-}+k_{0} \gamma_{y} \widehat{\mathcal{I}}_{1, N_{y}}^{-} X_{y;1}^{-}}{L_{y ; 1}^{-} } -\frac{G_{x}^{+}+k_{0} \gamma_{x} \widehat{\mathcal{I}}_{1, N_y}^{+} X_{x;1}^{+}}{L_{x ; 1}^{+} },\\ &A_4 = \frac{k_{0} H_{y}^{-}-\gamma_{y} \widehat{\mathcal{I}}_{1, N_{y}}^{-} L_{y ; 1}^{-}}{L_{y ; 1}^{-} } \frac{k_{0} H_{x}^{+}-\gamma_{x} \widehat{\mathcal{I}}_{1, N_y}^{+} L_{x, 1}^{+}}{L_{x ; 1}^{+} } ,\\ &F_v = \widetilde{f}_{R ; 1, N_{y}}-\frac{H_{y}^{-} g_{2 R ; 1, N_{y}}}{L_{y ; 1}^{-} } -\frac{\gamma_{y} \widehat{\mathcal{I}}_{1, N_{y}}^{-} X_{y ; 1}^{-} g_{2 I ; 1, N_{y}}}{L_{y ; 1}^{-} }\\ &\qquad-\frac{H_{x}^{+} g_{3 R ; 1, N_y}}{L_{x ; 1}^{+} } -\frac{\gamma_{x} \widehat{\mathcal{I}}_{1, N_y}^{+} X_{x ; 1}^{+} g_{3 I ; 1, N_y}}{L_{x ; 1}^{+} }, \end{align*}

    and for the vertex \Gamma_2\cap\Gamma_4 , there hold

    \begin{align*} &\mathbf{R_v} = (R_{N_x,N_y-1},R_{N_x-1,N_y},R_{N_x,N_y},I_{N_x,N_y}),\\ &A_1 = \frac{1}{L_{y ; 1}^{-} },\quad A_2 = \frac{1}{L_{x, 1}^{-} },\\ &A_3 = -\frac{G_{y}^{-}+k_{0} \gamma_{y} \widehat{\mathcal{I}}_{N_x, N_{y}}^{-} X_{y;1}^{-}}{L_{y ; 1}^{-} } -\frac{G_{x}^{-}+k_{0} \gamma_{x} \widehat{\mathcal{I}}_{N_{x}, N_y} X_{x;1}^{-}}{L_{x, 1}^{-} } ,\\ &A_4 = \frac{k_{0} H_{y}^{-}-\gamma_{y} \widehat{\mathcal{I}}_{N_x, N_{y}}^{-} L_{y ; 1}^{-}}{L_{y ; 1}^{-} } +\frac{k_{0}H_{x}^{-}-\gamma_{x} \widehat{\mathcal{I}}_{N_{x}, N_y} L_{x, 1}^{-}}{L_{x, 1}^{-} },\\ &F_v = \widetilde{f}_{R ; N_x, N_{y}}-\frac{H_{y}^{-} g_{2 R ; N_x, N_{y}}}{L_{y ; 1}^{-} } -\frac{\gamma_{y} \widehat{\mathcal{I}}_{N_x, N_{y}}^{-} X_{y ; 1}^{-} g_{2 I ; N_x, N_{y}}}{L_{y ; 1}^{-} }\\ &\qquad-\frac{H_{x}^{-} g_{4 R ; N_{x}, N_y}}{L_{x, 1}^{-} } -\frac{\gamma_{x} \widehat{\mathcal{I}}_{N_{x}, N_y}^{-} X_{x ; 1}^{-} g_{4 I ; N_{x}, N_y}}{L_{x, 1}^{-} }. \end{align*}

    For the imaginary part (50)-(51), the new finite difference scheme can be produced in the same way. And the details are omitted here.

    Remark 4. In fact, for the 2D equation, when the frozen-nonlinearity iteration and the modified Newton's method are used, it is not necessary to separate (5) into real and imaginary parts. In this case, the 2D equation can be divided into two 1D equations directly. Furthermore, since two 1D equations are separated from a 2D equation, we assume that \widetilde{f}_{xR},\widetilde{f}_{yR} are not discontinuous, and \varepsilon\neq0 in \widetilde{f} at the interface \partial\Omega_0 of the Kerr medium and the linear medium.

    In this section, we will show some numerical tests to verify the efficiency of the scheme proposed in the above section. And we set \delta = 10^{-8} in Algorithm 1.

    Firstly, let k_0,\; v and \varepsilon be given, and set the exact solution of the 1D equation to be (see [7])

    \begin{align} E& = \lambda(z) e^{i \varphi(z)}, \end{align} (63)

    where \lambda(z) and \varphi(z) are both real and satisfy

    \begin{align} \varphi(z)& = \varphi(0)+W \int_{0}^{z} \frac{1}{\beta(t)} d t, \end{align} (64)
    \begin{align} \lambda(z)& = \sqrt{\beta(z)}, \end{align} (65)

    with

    \begin{align} \beta(z)& = \beta_{2}+(\beta_{1}-\beta_{2}) \operatorname{cn}^{2}\left[\sqrt{\frac{\varepsilon}{2v}}\left(\beta_{1}-\beta_{3}\right)^{1 / 2} k_0\sqrt{v}(Z_{\max}-z)\bigg | \frac{\beta_{1}-\beta_{2}}{\beta_{1}-\beta_{3}}\right], \end{align} (66)
    \begin{align} \beta_{1}& = \frac{W}{k_{0}}, \end{align} (67)
    \begin{align} \beta_{2}& = \frac{v}{\varepsilon}\left\{\left[\left(1+\frac{\varepsilon W}{2 vk_{0}}\right)^{2}+\frac{2 \varepsilon W}{v^{2} k_{0}}\right]^{1 / 2}-\left(1+\frac{\varepsilon W}{2 vk_{0}}\right)\right\}, \end{align} (68)
    \begin{align} \beta_{3}& = -\frac{v}{\varepsilon}\left\{\left[\left(1+\frac{\varepsilon W}{2 vk_{0}}\right)^{2}+\frac{2 \varepsilon W}{v^{2} k_{0}}\right]^{1 / 2}+\left(1+\frac{\varepsilon W}{2 vk_{0}}\right)\right\}, \end{align} (69)

    and \text{cn}[\cdot|\cdot] being the Jacobi elliptic function and W\in[0,k_0] being a parameter need to be determined.

    Moreover, at z = 0 , there hold

    \begin{align} \frac{1}{2} \varepsilon \beta^{2}(0)+\left(v-1\right) \beta(0)+4-\left(v+3\right) \frac{W}{k_{0}}-\frac{1}{2}\frac{ \varepsilon W^{2}}{k_{0}^{2}} = 0, \end{align} (70)
    \begin{align} \frac{d\lambda}{d z}\bigg|_{z = 0} = 2 k_{0} \sin \varphi(0), \end{align} (71)

    with

    \begin{align} \left(\frac{d \lambda}{d z}\right)^{2}+\frac{W^{2}}{\lambda^{2}}+k_0^{2}v \lambda^{2}+\frac{1}{2} k_0^{2}\varepsilon \lambda^{4} = A. \end{align} (72)

    Thus, putting (66) with z = 0 into (70) and solving the nonlinear equation, we can determine W . Then, we obtain \beta(z_m)(1\leqslant m\leqslant N) by (66), A by A = k_{0}\left(v+1\right) W+\frac{ \varepsilon}{2} W^{2} , \lambda'(0) by (72), \varphi(0) by (71), and \varphi(z_m) by (64) one by one. Finally, putting \lambda(z_m) = \sqrt{\beta(z_m)} and \varphi(z_m) into (63), the exact solution is determined.

    To test the accuracy of the proposed scheme, with v = 1, k_{0} = 8, Z_{\max } = 10 and the initial guess E^{0}(z) = e^{i k_{0} z} , we firstly solve the 1D equation (3) with \varepsilon = 0.01 and 0.1 by using the new scheme based on the frozen-nonlinearity iteration method with k_0h = 1 . We can see that the numerical solutions can highly match the reference ones (see Figs. 3-4). Then, selecting frozen-nonlinearity iteration, the errors in l^{\infty} -norm with \varepsilon = 0.01 are compared in Table 1 among different numerical schemes: the standard finite difference (SFD) scheme, the compact finite difference (CFD) scheme, the finite volume (FV) method proposed in [2] and two new schemes (26)( k = 1 ) and (29). Clearly, the newly proposed finite difference schemes can achieve the best accuracy. These all imply that the newly proposed finite difference scheme can approximate the high oscillation solution of the NLH equation effectively.

    Figure 3.  Numerical solutions with \varepsilon = 0.01 in 1D problem (Red: new scheme with k_0h = 1 ; Blue: reference solution).
    Figure 4.  Numerical solutions with \varepsilon = 0.1 for the 1D problem (Red: new scheme with k_0h = 1 ; Blue: reference solution).
    Table 1.  Errors in l^\infty -norm for the 1D problem with \varepsilon = 0.01 .
    N 100 200 400 800 1600
    k_0=10
    SFD 2.14 1.05 2.69e-1 6.71e-2 1.67e-3
    FV[2] 1.59 5.03e-1 1.29e-1 3.23e-2 8.09e-3
    CFD 5.38e-1 3.70e-2 3.27e-3 6.67e-4 2.72e-4
    Scheme (29) 1.26e-3 2.99e-4 7.43e-5 1.89e-5 5.16e-6
    Scheme (30) 2.16e-4 5.68e-5 1.43e-5 3.68e-6 1.16e-6
    k_0=20
    SFD 2.31 2.13 1.80 5.33e-1 1.34e-1
    FV[2] 2.00 2.00 9.76e-1 2.60e-1 6.52e-2
    CFD 2.17 1.02 7.16e-2 5.46e-3 8.00e-4
    Scheme (29) 8.56e-3 1.57e-3 3.75e-4 9.57e-5 2.70e-5
    Scheme (30) 1.29e-3 3.16e-4 7.56e-5 1.87e-5 4.95e-6
    k_0=40
    SFD 1.24 2.35 2.13 2.03 1.03
    FV[2] - 2.00 1.99 1.70 5.16e-1
    CFD 1.22 2.36 1.76 1.40e-1 9.86e-3
    Scheme (29) 1.12e-2 9.70e-3 1.80e-3 4.49e-4 1.31e-4
    Scheme (30) 5.64e-3 2.68e-3 6.86e-4 1.05e-4 2.61e-5
    k_0=80
    SFD 1.07 1.05 2.32 2.29 2.02
    FV[2] - - 2.00 1.98 1.97
    CFD 1.04 1.21 2.31 1.99 0.29
    Scheme (29) 7.38e-3 8.68e-3 4.92e-3 1.99e-3 1.52e-3
    Scheme (30) 1.47e-3 1.05e-3 2.48e-4 2.56e-4 2.04e-4

     | Show Table
    DownLoad: CSV

    Then, under the same computational environments, but with Z_{\max } = 1 , the iteration numbers are compared among four iteration methods for varying \varepsilon and k_0 in Table 2. It can be found that, among the classical ones, the Newton's iteration method converges fastest, but the error correction has obvious advantage in decreasing the iteration number especially for the cases of large \varepsilon and k_0 . The values of k_0 or \varepsilon seem have little influence on the iteration number when the error correction method is convergent.

    Table 2.  Iteration numbers of different iteration methods for the 1D problem.
    k_0 10 20 40 80 160 320 640 1280
    \varepsilon=0.01
    Frozen-nonlinearity 5 5 6 7 9 12 17 38
    Error Correction 3 4 4 4 4 4 5 6
    Modified Newton 5 6 7 8 10 14 22 -
    Newton's method 4 4 5 5 6 8 11 -
    \varepsilon=0.02
    Frozen-nonlinearity 5 6 7 9 12 19 45 -
    Error Correction 4 4 4 4 5 5 7 9
    Modified Newton 6 7 8 10 14 23 - -
    Newton's method 4 5 5 6 8 11 - -
    \varepsilon=0.04
    Frozen-nonlinearity 6 8 9 13 22 55 - -
    Error Correction 4 4 5 5 6 8 13 -
    Modified Newton 7 9 10 16 25 - - -
    Newton's method 5 5 6 8 10 - - -
    \varepsilon=0.06
    Frozen-nonlinearity 7 9 12 18 35 - - -
    Error Correction 4 5 5 6 7 10 - -
    Modified Newton 8 9 14 20 39 - - -
    Newton's method 5 6 7 10 - - - -
    \varepsilon=0.08
    Frozen-nonlinearity 8 10 14 20 89 - - -
    Error Correction 5 5 6 7 9 - - -
    Modified Newton 9 11 17 25 - - - -
    Newton's method 5 6 8 11 - - - -
    \varepsilon=0.1
    Frozen-nonlinearity 9 10 16 35 - - - -
    Error Correction 5 6 6 8 12 - - -
    Modified Newton 10 13 18 36 - - - -
    Newton's method 6 7 9 - - - - -

     | Show Table
    DownLoad: CSV

    Furthermore, we simulate the optical bistability by using the proposed finite difference scheme. Firstly, letting the transmittance T: = E(Z_{\text{max}}) , we solve (3)-(4) with varying \varepsilon and plot |T|^2 in Fig. 5, the result is very similar to that in [2]. It worth noting that, since the Newton' s method is sensitive to the initial guess, we start the Newton's method with a initial guess which is obtained by solving (3)-(4) with the frozen-nonlinearity iteration method when \varepsilon = 0.01,E^0 = 0 . The solutions corresponding to the rest \varepsilon are obtained by increasing the value of \varepsilon step by step through some proper \Delta\varepsilon and the solution with \varepsilon_i is selected as the initial guess for \varepsilon_{i+1} = \varepsilon_i+\Delta\varepsilon . Then, we select much smaller \Delta\varepsilon to compute the switchback-type multi-solution near \varepsilon = 0.724 . In Fig. 6, we choose the solution with \varepsilon = 0.722 (point A) and \varepsilon = 0.726 (point E) as the initial guess to calculate the solutions corresponding to the lower branch and the upper branch, respectively. It worth to note that, in the lower branch, when we increase the value of \varepsilon at point F(\varepsilon\thickapprox0.72524) , the value of |T|^2 will directly jump to the one corresponding to point H and then varies following the route H\rightarrow E . Similarly, in the upper branch, decreasing the value of \varepsilon at point S(\varepsilon\thickapprox0.72376) leads to the value of |T|^2 follows the route S\rightarrow G\rightarrow A . Obviously, |E|^2 has three different values when \varepsilon\in(0.72376,0.72524) . To simualte the middle branch, we firstly set \varepsilon = 0.724226 and solve (3)-(4) to obtain the solution corresponding to point C with the initial guess which is selected as the linear combination of the solution at point B ( \varepsilon = 0.724226 in the lower branch) and the solution at point D ( \varepsilon = 0.724226 in the upper branch). Then, the middle branch is captured by selecting the solution corresponding to point C as the initial guess. Successfully approximating the optical bistability also indicates the robustness of the new finite difference method.

    Figure 5.  |T|^2 with respect to \varepsilon for the 1D problem.
    Figure 6.  Switchback-type non-uniqueness of |T|^2 near \varepsilon = 0.724 for the 1D problem.

    Now, we turn to a 2D problem. Setting \Omega = [-1 / 2,1 / 2]\times[-1 / 2,1 / 2], \Omega_{0} = [-1 / 4,1 / 4]\times[-1 / 4,1 / 4] and v = 1 in (6). In this case, we set the parameter \gamma_x = \gamma_y = 1/2 in (51)-(52) in our new finite difference scheme. Firstly, we examine the accuracy of the proposed finite difference scheme by taking \varepsilon = k_0^{-2} and the exact solution [29]

    \begin{align*} E = \frac{5 \sqrt{2} e^{\mathrm{i} y \sqrt{k_0^{2}+25}}}{\sqrt{\varepsilon} k_0 \cosh (5 x)}. \end{align*}

    In Fig. 7, we exhibit the exact solution and the numerical solution obtained by the new finite difference scheme with k_0 = 100,h = 1/400 , it can be observed that the solution is highly oscillating, but the numerical solution can match the exact solution well.

    Figure 7.  Solutions for the 2D problem with k_0 = 100,h = 1/400 (Left: numerical solution; Right: exact solution).

    To simulate the transmission and collision of the nonparaxial solitons which are also considered in [3,29], we solve the NLH equation (5) with two different incident waves

    \begin{align*} E_{\text {inc }}^1& = \frac{20 \sqrt{2} e^{i y \sqrt{k_{0}^{2}+400}}}{\cosh (20 x)},\\ E_{\text {inc }}^2& = \frac{20 \sqrt{2} e^{i y \sqrt{k_{0}^{2}+400}}}{\cosh (20 x)}+\frac{20 \sqrt{2} e^{i \sqrt{k_{0}^{2}+400}(y / 2-\sqrt{3} x / 2)}}{\cosh [20(x / 2+\sqrt{3} y / 2)]}. \end{align*}

    And the source term is set as

    \begin{align*} f = \left\{\begin{array}{c} -\Delta E_{\text {inc }}^l-k_{0}^{2} E_{\text {inc }}^l,(x, y) \in \Omega \backslash \Omega_{0}, \\ 0,\; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; (x, y) \in \Omega_{0}, \end{array}\right. l = 1,2. \end{align*}

    The intensities of the incident field and the total field for different cases are shown in Fig. 8 and Fig. 9. When only one incident wave E_{\text {inc }}^1 comes into \Omega from south (see Fig. 8(A)), it can be found that the incident wave can pass through the Kerr medium almost without any change when \varepsilon = k_0^{-2} (see Fig. 8(B)), while the field is totally different when \varepsilon = 0 (see Fig. 8(C)). When there are two solitons with the field E_{\text {inc }}^2 come into the medium from south and east respectively, they meet and collide in \Omega_0 . In the case of \varepsilon = k_0^{-2} , these two solitons are also nearly unchange when they pass through the Kerr medium (see Fig. 9(B)). It is worth to note that, comparing with Fig. 8(B), the amplitudes near (0,-1/2) and (1/2,0) have bigger values, which means that the backward scattering becomes stronger when two solitons colliding. Similarly, the total field changes a lot when the same medium is fulfilled in \Omega_0 and \Omega\setminus\Omega_0 (see Fig. 9(C)). Therefore, from the transmission and collision of the nonparaxial solitons examples, we can conclude that, the Helmholtz equation and the NLH equation are much different despite the nonlinear coefficient is very small (here \varepsilon = 10^{-4} ). And the scheme studied in this paper can release these clearly.

    Figure 8.  Transmission of a single soliton.
    Figure 9.  Collision of two solitons.

    In this paper, we construct a kind of new finite difference schemes for solving the nonlinear Helmholtz equation based on some iteration methods. Numerical results indicate that, the proposed scheme not only can approximate the high oscillation solution with better computational accuracy, but also can be used to simulate some physical phenomenons in the Kerr medium, such as the optical bistability and the collision of nonparaxial solitons. Moreover, without any extra consideration, this finite difference scheme also provides a route to deal with the problems with discontinuous coefficients or source terms. Thus, it can be extended to much more complex cases, such as the multi-layered Kerr mediums propagating problem and the nonlinear Maxwell's equations.

    [1] Friston K (2010) The free-energy principle: a unified brain theory? Nat Rev Neurosci 11: 127-138. doi: 10.1038/nrn2787
    [2] Jarzynski C (2012) Nonequilibrium equality for free energy differences. Phys Rev Lett 78: 2690-2693.
    [3] Landauer R (1996) Spatial variation of currents and fields due to localized scatterers in metallic conduction. IBM J Res Dev 1: 223-231.
    [4] Sengupta B, Stemmler MB, Friston KJ (2013) Information and efficiency in the nervous system--a synthesis. PLoS Comput Biol 9: e1003157. doi: 10.1371/journal.pcbi.1003157
    [5] Simons DJ, Chabris CF (2011) What people believe about how memory works: a representative survey of the U.S. population. PLoS One 6: e22757. doi: 10.1371/journal.pone.0022757
    [6] Friston K, FitzGerald T, Rigoli F, et al. (2016) Active inference and learning. Neurosci Biobehav Rev 68: 862-879. doi: 10.1016/j.neubiorev.2016.06.022
    [7] Friston K, Rigoli F, Ognibele D, et al. (2015) Active inference and epistemic value. Cogn Neurosci 6: 187-214. doi: 10.1080/17588928.2015.1020053
    [8] Ognibene D, Chinellato E, Sarabia M, et al. (2013) Contextual action recognition and target localization with an active allocation of attention on a humanoid robot. Bioinspir Biomim 8: 035002. doi: 10.1088/1748-3182/8/3/035002
    [9] Kahneman D (2011) Thinking, fast and slow. New York, NY.
    [10] Clark A, Chalmers DJ (1998) The extended mind. Analysis 58: 7-19. doi: 10.1093/analys/58.1.7
    [11] Kahneman D, Frederick S (2002) Representativeness revisited: Attribute substitution in intuitive judgment. In Gilovich T, Griffin D, Kahneman, D (Eds), Heuristics and biases: The psychology of intuitive judgment, New York, NY: 49-81.
    [12] Kahneman D (2003) Maps of bounded rationality: Psychology for behavioral economics. American Economic Review 93: 1449-1475. doi: 10.1257/000282803322655392
    [13] Simon HA (1955) A behavioral model of rational choice. Quarterly Journal of Economics 69: 99-118. doi: 10.2307/1884852
    [14] Lerner RG, Trigg GL (1991) Encyclopedia of Physics (2nd Edition). New York, NY: VCH Publishers.
    [15] Parker CB (1994) McGraw-Hill Encyclopedia of Physics (2nd Edition). New York, NY: McGraw-Hill.
    [16] Penrose R (2007) The Road to Reality: A Complete Guide to the Laws of the Universe. New York, NY: Vintage Books.
    [17] Panksepp J (2005) Affective consciousness: core emotional feelings in animals and humans. Conscious Cogn 14: 30-80. doi: 10.1016/j.concog.2004.10.004
    [18] Boltzmann L (1877) U ̈ er die beziehung zwischen dem zweiten haupt- satz der mechanischen wa ̈rmetheorie und der wahrscheinlichkeitsrech- nung respektive den sa ̈tzen u ̈ber das wa ̈rmegleichgewicht. [On the relationship between the second law of the mechanical theory of heat and the probability calculus]. Wiener Berichte 76: 373-435.
    [19] Wiener N (1961) Cybernetics-or control and communication in the animal and the machine. New York, NY.
    [20] Hirsh JB, Mar RA, Peterson JB (2012) Psychological entropy: A framework for understanding uncertainty-related anxiety. Psychol Rev 119: 304-320. doi: 10.1037/a0026767
    [21] Keysers C (2009) Mirror neurons. Curr Biol 19: 971-973. doi: 10.1016/j.cub.2009.08.026
    [22] Molenberghs P, Cunnington R, Mattingley JB (2009) Is the mirror neuron system involved in imitation? A short review and meta-analysis. Neurosci Biobehav Rev 33: 975-980.
    [23] Rizzolatti G, Craighero L (2004) The mirror-neuron system. Ann Rev Neurosci 27: 169-192. doi: 10.1146/annurev.neuro.27.070203.144230
    [24] Filimon F, Rieth CA, Sereno MI, et al. (2015) Observed, executed, and imagined action representations can be decoded from ventral and dorsal areas. Cereb Cortex 25: 3144-3158. doi: 10.1093/cercor/bhu110
    [25] Trapp K, Spengler S, Wüstenberg T, et al. (2014) Imagining triadic interactions simultaneously activates mirror and mentalizing systems. Neuroimage 98: 314-323. doi: 10.1016/j.neuroimage.2014.05.003
    [26] Felleman DJ, Van Essen DC (1991) Distributed hierarchical processing in the primate cerebral cortex. Cereb Cortex 1: 1-47.
    [27] von Stein A, Sarnthein J (2000) Different frequencies for different scales of cortical integration: from local gamma to long range alpha/theta synchronization. Int J Psychophysiol 38: 301-313. doi: 10.1016/S0167-8760(00)00172-0
    [28] Réka Albert AB (2002) Statistical mechanics of complex networks. Rev Mod Phys 74: 47-97. doi: 10.1103/RevModPhys.74.47
    [29] Amit DJ (1995) The Hebbian paradigm reintegrated: local reverberations as internal representations. Behav Brain Sci 18: 631.
    [30] Goldman-Rakic PS (1995) Cellular basis of working memory. Neuron 14: 477-485. doi: 10.1016/0896-6273(95)90304-6
    [31] Wang XJ (2001) Synaptic reverberation underlying mnemonic persistent activity. Trends Neurosci 24: 455-463. doi: 10.1016/S0166-2236(00)01868-3
    [32] Trantham-Davidson H, Neely LC, Lavin A, et al. (2004) Mechanisms underlying differential D1 versus D2 dopamine receptor regulation of inhibition in prefrontal cortex. J Neurosci 24: 10652-10659. doi: 10.1523/JNEUROSCI.3179-04.2004
    [33] Stephens GJ, Silbert LJ, Hasson U (2010) Speaker–listener neural coupling underlies successful communication. Proc Natl Acad Sci U S A 107: 14425-14430. doi: 10.1073/pnas.1008662107
    [34] Hasson U, Ghazanfa AA, Galantucci B, et al. (2012) Brain-to-brain coupling: a mechanism for creating and sharing a social world. Trends Cogn Sci 16: 114-121. doi: 10.1016/j.tics.2011.12.007
    [35] Benoit RG, Gilbert SJ, Frith CD, et al. (2012) Rostral prefrontal cortex and the focus of attention in prospective memory. Cereb Cortex 22: 1876-1886. doi: 10.1093/cercor/bhr264
    [36] Nader K, Schafe GE, Le Doux JE (2000) Fear memories require protein synthesis in the amygdala for reconsolidation after retrieval. Nature 406: 722-726. doi: 10.1038/35021052
    [37] Przybyslawski J, Sara SJ (1997) Reconsolidation of memory after its reactivation. Behav Brain Res 84: 241-246. doi: 10.1016/S0166-4328(96)00153-2
    [38] Sara SJ (2000) Retrieval and reconsolidation: toward a neurobiology of remembering. Learn Mem 7: 73-84. doi: 10.1101/lm.7.2.73
    [39] Colarusso CA, Nemiroff RA (1981) Adult Development: A New Dimension in Psychodynamic Theory and Practice. New York, NY: Plenum.
    [40] Pennisi E (2010) Conquering by copying. Science 328: 165-167. doi: 10.1126/science.328.5975.165
  • This article has been cited by:

    1. Jin Li, Barycentric rational collocation method for semi-infinite domain problems, 2023, 8, 2473-6988, 8756, 10.3934/math.2023439
    2. Jin Li, Yongling Cheng, Barycentric rational interpolation method for solving KPP equation, 2023, 31, 2688-1594, 3014, 10.3934/era.2023152
    3. Jin Li, Kaiyan Zhao, Xiaoning Su, Selma Gulyaz, Barycentric Interpolation Collocation Method for Solving Fractional Linear Fredholm-Volterra Integro-Differential Equation, 2023, 2023, 2314-8888, 1, 10.1155/2023/7918713
    4. Jin Li, Yongling Cheng, Barycentric rational interpolation method for solving fractional cable equation, 2023, 31, 2688-1594, 3649, 10.3934/era.2023185
    5. Jin Li, Yongling Cheng, Barycentric rational interpolation method for solving 3 dimensional convection–diffusion equation, 2024, 304, 00219045, 106106, 10.1016/j.jat.2024.106106
    6. Jin Li, Yongling Cheng, Barycentric rational interpolation method for solving time-dependent fractional convection-diffusion equation, 2023, 31, 2688-1594, 4034, 10.3934/era.2023205
    7. Jin Li, Yongling Cheng, Spectral collocation method for convection-diffusion equation, 2024, 57, 2391-4661, 10.1515/dema-2023-0110
    8. Jin Li, Linear barycentric rational interpolation method for solving Kuramoto-Sivashinsky equation, 2023, 8, 2473-6988, 16494, 10.3934/math.2023843
  • Reader Comments
  • © 2017 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(4892) PDF downloads(1109) Cited by(0)

Figures and Tables

Figures(2)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog