Hysteresis-driven pattern formation in reaction-diffusion-ODE systems

Figure 1.  Typical configurations of the zero sets of the kinetic functions

Figure 2.  (A) Phase plane of $\gamma^{-1}U_{xx}+q \bar{u}(U)=0$. The blue trajectory $(U(x),U_x(x))$ connects the points $(u_0,0)$ and $(u_e,0)$ and is a solution of the boundary value problem satisfying $U_x(0)=U_x(1)=0$. (B) A monotone increasing stationary solution with jump at $\bar{u}$ and layer position at $\bar{x}$

Figure 3.  The time-maps for the kinetic functions $f(u,v)=1.4v-u$, $g(u,v)=u-(v^3-6.3v^2+10v)$ and the jump $\bar{u}=3.1$. Here, it holds $Q \bar{u}(u_2)<0$ and u_min=0.8957. We see that T_u¹ is decreasing, whereas $T_{\bar{u}}^{2}$ is increasing, which leads to $\frac{\mathrm{d}}{\mathrm{d} u_{0}} T_{\bar{u}}(u_0)<0$. Furthermore, we observe that $\lim_{u_0\to {\bar{u}}}T_{\bar{u}}{}(u_0)=0$ and $\lim_{u_0\to u_{\rm{min}}}T_{\bar{u}}{}(u_0)=\infty$

Figure 4.  Simulations of the generic model with hysteresis for different types of perturbations of a stationary solution. The plots show initial conditions (dotted lines) and the approached stationary solution (continuous lines) after a sufficiently large time $ t_{end} $

Figure 5.  The layer position $ \bar{x}( \bar{u}) $ and the Interval $ I^0 $. (A) and (B) Plots for the kinetic functions eq. (60) and for diffusion coefficients $ 1/\gamma $ with $ \gamma=50 $ in A and $ \gamma=200 $ in B. (C) Plots for the kinetic functions eq. (61) and the diffusion coefficient $ 1/200 $

Figure 6.  The phase planes of $\frac{1}{\gamma} U_{x x}+q_{H}(U)=0$ and $\frac{1}{\gamma} U_{x x}+q{T}(U)=0$ are overlapping. In blue we see a periodic solution with jump at u. We cannot determine the mode of a periodic solution in the phase plane. It corresponds to how often the trajectory has been traveled through. In red we see a irregular solution with three different jumps

Figure 7.  An irregular solution $ \big(U(x), V(x)\big) $ with three jumps $ \bar{u}^{{1}}, \bar{u}^{{2}}, \bar{u}^{{3}} $, which is monotone increasing restricted to $ [0, x^1] $. We see that for continuity of $ U(x) $ we need to have $ u_{e}^1=u_0^2 $ and $ u_{e}^2=u_0^3 $ fulfilled. Furthermore, we see how the partition of the interval is determined: $ x^1=T(\bar{u}^{{1}}, u_0^1), x^2=x^1+T(\bar{u}^{{2}}, u_0^2) $ and $ 1=x^2+T(\bar{u}^{{3}}, u_0^3) $. The layer positions are given by $ \bar{x}^{{1}}=T_1(\bar{u}^{{1}}, u_0^1) $ and $ \bar{x}^{{3}}=x^2+T_1(\bar{u}^{{3}}, u_0^3) $, because $ U(x) $ is increasing on the corresponding subintervals and $ \bar{x}^{{2}}=x^1+T_2(\bar{u}^{{2}}, u_0^2)=x^2-T_1(\bar{u}^{{2}}, u_0^3) $

Figure 8.  Simulations of model (1)-(3) for admissible kinetic functions, diffusion coefficient $ 1/\gamma=1/1000 $ and initial conditions of type (59) having four discontinuities. The $ u $-component is plotted in blue, whereas the $ v $-component is red. The initial condition $ \big(u(0, x), v(0, x)\big)=\big(u_0(x), v_0(x)\big) $ is indicated by dotted lines and the stationary solution $ \big(u(t_{end}, x), v(t_{end}, x)\big) $ is indicated by continuous bold lines. Here, $ t_{end} $ is a sufficiently large timepoint, such that the solution $ \big(u(t, x), v(t, x)\big) $ does not change in time anymore. A-C the kinetic functions are given by (60) D the kinetic functions are given by (62)