1. Introduction
We consider the Kobayashi–Warren–Carter system, introduced in [Reference Kobayashi, Warren and Carter14, Reference Kobayashi, Warren and Carter15, Reference Warren, Kobayashi and Carter29], to model evolutions of structures in a multi-grain problem. It is a kind of phase-field system in a domain Ω in Rn, formally a gradient flow of the energy


Here, $\alpha_0(v)\geq0$ is a given function, typically
$\alpha_0(v)=rv^2$, with a constant r > 0, and F(v) is a single-well potential, typically
$F(v)=a^2(v-1)^2$, with a constant a > 0. The functional
$E_\mathrm{sMM}^\varepsilon$ is often called a single-well Modica–Mortola functional. The Kobayashi–Warren–Carter system is regarded as a gradient flow with respect to L 2-inner product

where $\alpha_w\geq0$ and τ > 0 are weights. The function
$\alpha_w=\alpha_w(v)$ is given, but it may depend on
$v\in L^2(\Omega)$, so the above inner product is a Riemann metric on the tangent bundle
$TL^2(\Omega)$. A typical form of
$\alpha_w(v)$ equals
$\alpha_w(v)=\tau_0v^2$, where τ 0 is a positive constant. We consider the gradient flow of
$E_\mathrm{KWC}^\varepsilon$ under this metric, and its explicit form is

An explicit form in [Reference Kobayashi, Warren and Carter14] corresponds to the case when $s=\varepsilon$,
$\nu=\varepsilon^2$,
$\tau_1=\tau\varepsilon$,
$F(v)=(v-1)^2$,
$\alpha_0(v)=v^2$,
$\alpha_w(v)=\tau_0 v^2/\varepsilon$ with
$\tau_0 \gt 0$. In other words,

The function v represents an order parameter, where v = 1 corresponds to a grain region and where v away from 1 corresponds to grain boundaries. The function u represents a structure-like averaged angle in each grain.
We are interested in a singular limit problem for Eqs. (1.3) and (1.4) as $\varepsilon\downarrow0$. It turns out that the correct scaling of time should be
$\tau=\tau_1/\varepsilon$, while τ 1 is independent of ɛ. Since the system Eqs. (1.3)–(1.4) is regarded as a gradient flow of
$E_\mathrm{KWC}^\varepsilon(u,v)$ of Eq. (1.1), we are tempted to expect that the limit flow is the gradient flow of its limit energy
$E_\mathrm{KWC}^0$ which was obtained in our papers [Reference Giga, Okamoto, Sakakibara and Uesaka8, Reference Giga, Okamoto and Uesaka9]. Surprisingly, this conjecture is wrong. The limit flow contains a fractional time derivative. In this paper, we consider the problem in a one-dimensional setting. Moreover, we consider a special but typical case when the problem is essentially reduced to a single equation for v in Eq. (1.3) because handing Eq. (1.4) is technically involved since it is a total variation flow type equation. This reduced problem becomes a linear problem and is easy to discuss.
We consider Eqs. (1.3)–(1.4), where Ω is an interval $\Omega=(-L,L)$, and impose the Dirichlet boundary condition for u and the Neumann boundary condition for v. More precisely,

while

We set

We expect that the function

with b > 0, solves Eq. (1.4). Since Eq. (1.4) is of total variation flow type, the definition of a solution is not obvious. Fortunately, under a suitable assumption of v, say $v(0,t)\leq v(x,t)$ for all
$x\in(-L,L)$, t > 0, the function ub solves Eq. (1.4) under Eq. (1.7), as shown in the following lemma by setting
$\beta=\alpha_0(v)$.
Lemma 1.1. (A stationary solution)
Assume that $\beta\in C[-L,L]$ satisfies

Then, ub solves

under Eq. (1.7).
We stress that the notion of a solution of the equation for u in lemma 1.1 and also Eq. (1.4) under the Dirichlet condition Eq. (1.7) is not obvious and will be discussed in §4.
Problem Eqs. (1.3)–(1.4) are reduced to

under the boundary condition

and the initial condition

where $1_{x \gt 0}$ is a characteristic function of
$(0,+\infty)$, i.e., the Heaviside function, so that its distributional derivative equals the Dirac δ function. We note that if we assume that
$v_0^\varepsilon$ is even and non-decreasing for x > 0, i.e.,
$v_0^\varepsilon(x)=v_0^\varepsilon(-x)$ and
$v_{0x}^\varepsilon(x)x\geq0$ for
$x\in(-L,L)$, then ub is indeed a stationary solution of Eq. (1.4) with Eq. (1.7). This follows from lemma 1.1 since we have
$v^\varepsilon(0,t)\leq v^\varepsilon(x,t)$ for all
$x\in(-L,L)$, t > 0 by the maximum principle. By a scaling transformation
$y=x/\varepsilon$, Eq. (1.10) becomes

in $(-L/\varepsilon,L/\varepsilon)\times(0,\infty)$ for
$V=V^\varepsilon(y,t)=v^\varepsilon(\varepsilon y,t)$, where v ɛ is a solution of Eq. (1.10). Thus, we expect that this limit
$V=\lim_{\varepsilon\to0}V^\varepsilon$ solves Eq. (1.13) on R and is bounded. Since the solution v ɛ of Eq. (1.10) is expected to converge to 1 except at x = 0, we are interested in the behaviour of
$\xi^\varepsilon(t)=v^\varepsilon(0,t)$. More precisely, we would like to find the equation which
$\xi=\lim_{\varepsilon\to0}\xi^\varepsilon$ solves.
We set

Since $F(v)=a^2(v-1)^2$ to get
$G(\zeta)=a(\zeta-1)^2/2$, we have

As proved in [Reference Giga, Okamoto and Uesaka9], this energy is obtained as a value at Ξ of Γ-limit of $E_\mathrm{sMM}^\varepsilon(v)$ in Eq. (1.2) as ɛ → 0 if v converges to a set-valued function Ξ (under the graph convergence) of the form

for $\zeta\in(0,1)$.
A key observation is to derive an equation for ξ. For $\alpha\in\mathbf{R}$, β > 0, we set

where Γ denotes the gamma function. We consider well-prepared initial data in the sense that it is continuous and solves Eq. (1.13) outside y = 0.
Lemma 1.2. (Limit equation)
Assume that $\tau_1=1$. Let V be the bounded solution of Eq. (1.13) in
$\mathbf{R}\times(0,\infty)$ with initial data
$V(y,0)=1-ce^{-a|y|}$ with some
$c\in\mathbf{R}$. Then,
$\xi(t)=V(0,t)$ solves

with

Moreover, $m_a^\prime(t) \lt 0$ and
$m_a(t) \gt 0$, and
$\lim_{t\to\infty}m_a(t)=0$.
The assumption $\tau_1=1$ is just for the convenience of presentation. The formula for general τ 1 is obtained by rescaling the time variable t by
$t^\prime=t\tau_1$.
Note that in case a = 0, the left-hand side of Eq. (1.14) becomes the Caputo derivative $2\partial_t^{1/2}$. For even initial data, the solution V of Eq. (1.13) is even. Since Vy is odd in y, we see that the delta part of Vyy equals

where $V(\pm0,t)=\lim_{\gamma\to\pm0}V(\gamma,t)$ and δ denotes the Dirac delta function. Since V is smooth outside y = 0 and Vt has no delta part, Eq. (1.13) deduces that

Thus, Eq. (1.13) in $\mathbf{R}\times(0,\infty)$ is reduced to the Robin boundary problem in
$(-\infty,0)\times(0,\infty)$ with

The Caputo derivative appears in the equation of the boundary value.
Corollary 1.3. Let $w=w(x,t)$ be the bounded solution of the heat equation

with the Robin boundary condition

Assume that $w(x,0)=-c$ for some
$c\in\mathbf{R}$; the boundary value
$\xi(t)=w(0,t)$ solves

In other words, $\partial_t^{1/2}\xi=-b\xi$, where
$\partial_t^{1/2}$ is the Caputo half derivative.
It is well known that the fractional Laplace operator $(-\Delta)^{1/2}$ arises as the Dirichlet–Neumann map of the Laplace equation. Here, the Caputo derivative
$\partial_t^{1/2}$ is obtained as the Dirichlet–Neumann map of the heat equation. Formally, it is easy to guess since the Robin boundary condition yields
$\partial_t^{1/2}\xi+b\xi=0$ by replacing
$\partial_x$ with
$\partial_t^{1/2}$, which is natural since
$\partial_t=\partial_x^2$ for w. In a seminal paper, Caffarelli and Silvestre [Reference Caffarelli and Silvestre3] show that
$(-\Delta)^\gamma$ (
$0 \lt \gamma \lt 1$) is obtained as the Dirichlet–Neumann map for the degenerate Laplace equation. We remark that
$\partial_t^\gamma$ is obtained as the Dirichlet–Neumann map of a degenerate heat equation
$w_t-(-x)^\alpha w_{xx}=0$ with
$\gamma=1/(2-\alpha)$ as in [Reference Caffarelli and Silvestre3]; see remark 2.4 at the end of §2. We do not pursue this problem in this paper.
Since Eq. (1.13) is linear and of constant coefficients, we can use the Laplace transform to obtain the desired equation. As one expects, we are able to prove the convergence of V ɛ.
Lemma 1.4. (Convergence)
Let V be the bounded solution of Eq. (1.13) in $\mathbf{R}\times(0,\infty)$ with initial data V 0, which is bounded and uniformly continuous on R. Let vɛ be the solution of Eq. (1.10) under Eqs. (1.11) and (1.12). Assume that
$V_0^\varepsilon=\left.V^\varepsilon\right|_{t=0}\in C\left[-L/\varepsilon,L/\varepsilon\right]$ converges to V0 in the sense that

Then, Vɛ converges to V locally uniformly in $\mathbf{R}\times[0,\infty)$. In particular,
$\xi^\varepsilon(t)=v^\varepsilon(0,t)=V^\varepsilon(0,t)$ converges to
$\xi(t)=V(0,t)$ locally uniformly in
$[0,\infty)$. If we assume that

then

for any T > 0. Moreover,

Although the statement is for linear equations, we have to be a little bit careful since the domains where solutions are defined depend on ɛ. We extend the initial data $V_0^\varepsilon$ to a whole space. If the solution W ɛ to the whole space problem with the extended data
$\overline{V_0^\varepsilon}$ solves the Neumann problem in
$(-L/\varepsilon,L/\varepsilon)$, it is very convenient. We take this strategy. However,
$\overline{V_0^\varepsilon}$ may not be close to V 0 outside
$(-L/\varepsilon,L/\varepsilon)$, so the maximum principle does not imply the convergence of the solution W ɛ starting from
$\overline{V_0^\varepsilon}$ to V even locally. We regularize initial data and prove that
$\{W^\varepsilon\}$ is bounded and equi-continuous in
$\mathbf{R}\times[0,T]$,
$0 \lt T \lt \infty$. By Arzelà–Ascoli theorem and a diagonal argument, it converges to some function by taking a subsequence locally uniformly in
$\mathbf{R}\times[0,T]$. The limit is identified with a solution V, so it is a full convergence. This eventually yields that V ɛ converges to V locally uniformly in
$\mathbf{R}\times[0,\infty)$, so it yields the convergence of ξ ɛ. Because of the term
$V\partial_y(1_{y \gt 0})$, the actual proof is more involved. We divide the initial data into even part and odd part. For the odd part, the problem is reduced to the problem without the
$V\partial_y(1_{y \gt 0})$ term. For the even part, the problem is reduced to the Robin boundary problem and the proof is more involved. Note that we do not assume behaviour of
$V_0(y)$ so
$|y|\to\infty$ to derive locally uniform convergence.
To obtain the uniform convergence of V ɛ to V in $(-L/\varepsilon,L/\varepsilon)\times(0,T)$, it seems that the behaviour of
$V_0(y)$ as
$|y|\to\infty$ should be controlled. We observe that our decay assumption for
$V_0(y)-1$ is preserved not only for V but also for
$\partial_yV$ for
$t \gt \delta \gt 0$. This allows us to compare V ɛ with V globally in space. All proofs are quite elementary, but we give full proofs for the reader’s convenience.
Applying lemmas 1.2 and 1.4, we can obtain a characterization of the limit equation.
Theorem 1.5 Assume that $\tau_1=1$. Let vɛ be the solution of Eq. (1.10) under Eqs. (1.11) and (1.12). Assume that
$v_0^\varepsilon$ is well-prepared in the sense

as $\varepsilon \to 0$, with some c independent of ɛ. Then,
$\xi^\varepsilon(t)$ converges to ξ locally uniformly in
$[0,\infty)$, and ξ solves Eq. (1.14). Moreover, the graph of vɛ converges to a set-valued function Ξ of the form

The convergence is in the sense of the Hausdorff distance of graphs over $[-L,L]\times[0,T]$ for any T > 0.
We also handle initial data not necessarily well-prepared. In this case, Eq. (1.14) is altered because there are a few lower-order terms; see Eq. (2.9). If $V(y,0)=1-ce^{-\mu|x|}$ with
$\mu\neq a$, we still get an explicit form corresponding to Eq. (1.14); see Eq. (2.16). In these cases, the solution ξ is explicitly represented by using the error function. We shall discuss these extensions for non-well-prepared data and the proof of lemma 1.2 in §2. We also calculate numerically how the solution of the ɛ-problem converges by comparing it with the explicit solution of Eq. (1.14) and more general Eq. (2.16).
We now come back to the singular limit problem of the Kobayashi–Warren–Carter system Eqs. (1.3)–(1.4), with $\tau=\tau_1/\varepsilon$,
$F(v)=a^2(v-1)^2$, and
$\alpha_0(v)=v^2$ as ɛ → 0. We give a purely formal argument. Eq. (1.4) is of the total variation flow type, and its well-posedness for given v is known when v is independent of time and
$\inf\alpha_w(v) \gt 0$; in this case, Eq. (1.4) is the gradient flow of the weighted total variation
$\int\alpha_0(v)|\nabla u|$ under
$\alpha_w(v)$-weighted L 2 inner product if we impose the natural boundary condition like the Neumann boundary condition. As in [Reference Giga, Giga and Kobayashi4], we consider Eq. (1.4) for piecewise constant functions. For an interval
$I=(p,q)$ and its given division
$p=p_0 \lt p_1 \lt \cdots \lt p_m=q$, we consider a piecewise constant function of form

We interpret a solution by mimicking the notion of a solution when v is independent of time and $\inf\alpha_w(v) \gt 0$ (under the periodic boundary condition for simplicity). It is of the form

where $\chi_j=\operatorname{sgn}(h_{j+1}-h_j)$; we identify
$p_m=p_0$ by periodicity. See figure 1. The function z is called a Cahn–Hoffman vector field.

Figure 1. Graph of u with $\chi_1=1, \chi_2=-1 \text{ and } \chi_3=-1$
For a given v, it is unclear whether the solution remains in a class of piecewise constant functions [Reference Giga, Giga and Kobayashi4], known as a facet-splitting problem. If v is constant, it is well known that the solution remains in a class of piecewise constant functions.
We postulate that our system Eq. (1.15) has a (spatially) piecewise constant solution. Integrating the first equation of Eq. (1.15) in $(p_{j-1},p_j)$ yields

If $v_\varepsilon\to\Xi$ in the graph topology as ɛ → 0, we arrive at

where the set-valued function Ξ is defined as

Since jump $b_j=|h_{j+1}-h_j|$ determines the ξj-equation as indicated in theorem 1.5, the equation for ξj (assuming
$\tau_1=1$) is expected to be

where $M_a f=\int_0^t m_a(t-s)f(s)\,ds$ for well-prepared initial data at least when
$h_{j+1}-h_j$ is independent of time t. We postulate that Eq. (1.17) is still valid when hj depends on time. Thus, the singular limit equation of the Kobayashi–Warren–Carter equations Eqs. (1.3)–(1.4) as ɛ → 0 (under the periodic boundary condition) is expected to be the system Eqs. (1.16)–(1.17) for hj and ξj. If we consider other boundary conditions, we always impose the homogeneous Neumann boundary condition for v, like Eq. (1.8). If we consider the Dirichlet boundary condition for u with time-independent data, the values of h 1 and hm are prescribed. We can impose the Neumann condition for u; in this case, imposing
$z(p_0)=z(p_m)=0$ in Eq. (1.15) is natural. It is rather standard [Reference Podlubny23] to construct a unique local-in-time solution for a system of a fractional differential equation and an ordinary differential equation.
We note that the solvability of the initial value problem for the original Kobayashi–Warren–Carter system Eqs. (1.5)–(1.6) is still an open problem, even in a one-dimensional setting. If we write it in the form of Eqs. (1.3)–(1.4), the difficulty stems from the fact that the weights αw and α 0 can vanish somewhere. In the literature, αw is assumed to be away from zero. If α 0 is allowed to vanish, the $\Delta u$ term is added in the right-hand side of Eq. (1.6). For example, instead of considering Eq. (1.6), we consider

with δ > 0, $\delta^\prime\geq0$ and
$\mu\geq0$ such that
$\delta^\prime+\mu \gt 0$. The existence of a solution to Eqs. (1.5) and (1.18) with its large time behaviour is established in [Reference Ito, Kenmochi and Yamazaki10–Reference Kenmochi and Yamazaki13, Reference Moll and Shirakawa19, Reference Moll, Shirakawa and Watanabe20, Reference Shirakawa and Watanabe26, Reference Shirakawa, Watanabe and Yamazaki27, Reference Watanabe and Shirakawa30] under several homogeneous boundary conditions. (The case δ = 0 is included in [Reference Ito, Kenmochi and Yamazaki11, Reference Watanabe and Shirakawa30].) Unfortunately, the uniqueness of their solution is only known in a one-dimensional setting under the relaxation term µ > 0 [Reference Ito, Kenmochi and Yamazaki10, Theorem 2.2]. The extension of these results to inhomogeneous boundary condition is not difficult. In [Reference Moll, Shirakawa and Watanabe21], under non-homogeneous Dirichlet boundary conditions, structured patterns of stationary (i.e., time-independent) solutions were studied. In a one-dimensional setting, they thoroughly characterized all stationary solutions. Our Γ-convergence result in [Reference Giga, Okamoto and Uesaka9] gives the convergence of minimizer of
$E_\mathrm{KWC}^\varepsilon$ to that of the limit energy as ɛ → 0. We do not know the convergence of stationary solutions.
This paper is organized as follows. In §2, we derive the limit equation for both well-prepared and non-well-prepared initial data. We also give the large time behaviour of a solution. Most of the calculations are very explicit. In §3, we prove lemma 1.4. Section 4 gives a rigorous definition of the Dirichlet problem for Eq. (1.4), assuming that v is given and $\alpha_w(v)\equiv1$. In §5, we give several numerical tests.
2. Derivation of equations with the fractional time derivative
The first goal of this section is to prove lemma 1.2. We then consider more general initial data. We in particular consider the case $V(y,0)=1-ce^{-\mu|y|}$ with µ > 0 not necessarily equal to a and derive the equation corresponding to Eq. (1.14); see Eq. (2.16). The solution ξ can be written explicitly (see Eq. (2.18)), and it is of the form

with a lower-order term

This asymptotic result is shown for general µ; see lemma 2.3 for a detailed statement. At the end of this section, we give remarks related to corollary 1.3. It turns out that corollary 1.3 follows from an expression of the Dirichlet–Neumann map of the heat equation. The Caputo derivatives can be naturally derived using the Dirichlet–Neumann map.
We begin with recalling several elemental properties of the Laplace transform

for a locally integrable function g in $[0,\infty)$. By definition,

and by definition of the Gamma function, we see

We now arrive at a well-known formula

for $\alpha,\beta \gt 0$.
Lemma 2.1. For a > 0, let $m_a(t)$ be

Then

Moreover, $m_a^\prime(t) \lt 0$ and
$m_a(t) \gt 0$ for t > 0, with
$\lim_{t\to\infty}m_a(t)=0$.
Proof. We recall that the Laplace transform of the convolution $(g_1*g_2)(t)=\int_0^t g_1(t-s)g_2(s)\,ds$ is given by

If we take $g_1\equiv1$, we see

since $\mathcal{L} [g_1](\lambda)=\lambda^{-1}$ by Eq. (2.2). Setting
$\alpha=a^2$, we apply Eqs. (2.2) and (2.4) to get

which yields the desired formula for $\mathcal{L}[m_a]$.
We differentiate $m_a/2$ to get

We observe that

Thus, $\lim_{t\to\infty}m_a(t)=0$ since
$\lim_{t\to\infty}f_{1/2}^\alpha(t)=0$, which implies that
$m_a(t) \gt 0$ for t > 0 since
$m_a^\prime(t) \lt 0$.
Proof of lemma 1.2
Since the property of ma has been proved in lemma 2.1, it suffices to derive Eq. (1.14). Studying the equation for $w=V-1$ instead of Eq. (1.13) is more convenient. Eq. (1.13) for w (with
$\tau_1=1$) becomes

The initial data $w(x,0)$ equal

with $w_0(x)=-ce^{-a|x|}$. Let
$\widehat{w}(x,\lambda)$ be the Laplace transform of w in the t variable, i.e.,

We note that

Taking the Laplace transform of Eq. (2.5), we arrive at

a second-order linear ordinary differential equation with a jump of the derivatives $\widehat{w}_x$. Since the coefficients are constants, we can solve Eq. (2.6) explicitly.
A bounded solution satisfying

is of form

with $A\in\mathbf{R}$. One should determine A so that Eq. (2.6) holds. We observe

so that

since $(\operatorname{sgn}x)^\prime=2\delta$. Eq. (2.6) imposes that

which implies

Since $\widehat{w}(0,\lambda)=A-c\lambda^{-1}$, we end up with

We set $\eta(t)=w(0,t)$ and observe that

We next extract the time derivative of η. We note that

so that $\widehat{\eta_t}=\lambda\widehat{\eta}+c$ since
$\eta(0)=-c$. Then,

Multiplying $\sqrt{\lambda+a^2}+b$ and subtracting
$(a+b)\widehat{\eta_t}$ from both sides, we get

The energy for ξ equals

Since $\eta=\xi-1$, we see
$\operatorname{grad}E_\mathrm{sMM}^{0,b}(\xi)=2(b+a)\eta+2b$. Thus,

By lemma 2.1 and Eq. (2.3), we now conclude that $\xi=\eta+1$ solves Eq. (1.14) by taking the inverse Laplace transform of Eq. (2.7).
In lemma 1.2, we only considered well-prepared initial data, meaning that $V(y,t)$ is a stationary solution of Eq. (1.13) in
$\mathbf{R}\backslash\{0\}$. We take this opportunity to write a general equation corresponding to Eq. (1.14) starting from general initial data. We set

and Eq. (1.14) is

For general initial data, our equation corresponding to Eq. (2.8) becomes more complicated than Eq. (2.8).
We set

which is the Green function of $-\partial_y^2+\sigma$ with
$\sigma=\lambda+a^2$, i.e.,

For $w_0\in L^\infty(\mathbf{R})$, we set

where $*_x$ denotes the convolution in space, i.e.,

Lemma 2.2. Let $\tau_1=1$. Let V be the bounded solution of Eq. (1.13) in
$\mathbf{R}\times(0,\infty)$ with initial data
$V(y,0)=1+w_0$, where w 0 is bounded and Lipschitz continuous in R. Then,
$\xi(t)=V(0,t)$ solves

Proof. If w 0 is Lipschitz, then we easily see that $|w_t|t^{1/2}$ is bounded (near t = 0) so that the integrand of

is integrable for all $x\in\mathbf{R}$ near s = 0; see lemma 3.4.
We argue in the same way to prove lemma 1.2. We begin with Eq. (2.6), where $w=V-1$. Its solution outside x = 0 is given as

with $A\in\mathbf{R}$. As before, we determine A and obtain that
$\eta=\xi-1$ satisfies

Since $\widehat{\eta_t}=\lambda\widehat{\eta}-\eta(0)$, we proceed with

Thus,

Taking the inverse Laplace transform, we obtain

the same as Eq. (2.9).
The term $\mathcal{L}^{-1} \left[2\sqrt{\lambda+a^2} g^a\right]$ has a more explicit form. Let
$E(x,t)$ be the Gauss kernel, i.e.,

We know that its Laplace transform (as a function of t) is

Let $E^a=e^{-a^2 t}E$, then

Since $\mathcal{L}\left[f_{1/2}^{a^2}\right]=(\lambda+a^2)^{-1/2}$ and
$\mathcal{L}[g_1*g_2] = \mathcal{L}[g_1] \mathcal{L}[g_2]$, we end up with

Thus,

Since $\sqrt{\lambda+a^2}g^a=\frac12\left(e^{-\sqrt{\lambda+a^2}|x|} *_x w_0 \right)(0,t)$, we conclude that

If $w_0=-ce^{-\mu|x|}$, we observe that

where $\operatorname{erfc}$ denotes the complementary error function, i.e.,

Indeed, we proceed with

By a direct calculation, we observe that

Thus,

where we invoked Eqs. (2.2) and (2.1). We set $q^\mu(t)=f_{1/2}^0(t)-\mu e^{\mu^2 t}\operatorname{erfc}\left(\mu\sqrt{t}\right)$ so that Eq. (2.13) becomes

Using Eq. (2.1) again, we now obtain Eq. (2.11).
The formula (2.13) gives another representation of ma. Indeed, since

we see

in particular, implies that $m_a(t)\leq 2f_{1/2}^{a^2}(t)$ and

Moreover, by Eq. (2.15), we see

Therefore, the positivity of ma (lemma 2.1) implies $q^a(t) \gt 0$ for t > 0. The formula Eq. (2.9) in lemma 2.2 becomes

if $w_0=-ce^{-\mu|x|}$. If
$\mu=a$, this is reduced to Eq. (1.14) (or Eq. (2.8)).
One can give an explicit form of a solution of Eq. (2.16) starting from $\xi(0)=1-c$. We substitute Eq. (2.12) into Eq. (2.10) to get

Since

by Eqs. (2.13) and (2.1), we have

However, the calculation of $q^b*q^\mu$ is quite involved, and it is easier to calculate
$\widehat{\eta}$ in Eq. (2.17) more. We proceed

It is easy to see that

As we already observed,

For the first term,

By definition,

We thus conclude that

Thus, $\xi=\eta+1$ is the solution of Eq. (2.16) with
$\xi(0)=1-c$.
From this solution formula, we can establish the solution’s large-time behaviour. We set

(i) The function
$\bar{\eta}$ is negative, monotonically decreasing, and converging to
$-b/(b+a)$ as
$t\to\infty$.
(ii) The estimate
\begin{align*} \left| \frac{\eta_e}{c(\mu+b)} \right| \leq e^{-a^2t} \int_t^\infty q^b(s)\,ds \leq e^{-a^2t} \frac1b \end{align*}
holds for t > 0. In particular,
\begin{align*} \lim_{t\to\infty}\eta_e e^{a^2t} = 0. \end{align*}
(i) Since
$q^b\geq0$, the monotonicity is clear. We observe that
\begin{align*} \int_0^\infty e^{-a^2s}q^b(s)\,ds = \mathcal{L} \left[q^b\right] (0+a^2) = \frac{1}{\sqrt{a^2}+b} \end{align*}
by Eq. (2.13). Thus,
$\lim_{t\to\infty} \bar{\eta}(t)=-b/(b+a)$.
(ii) Since
\begin{align*} \int_0^\infty e^{-\mu^2s}q^b(s)\,ds = \frac{1}{\mu+b}, \end{align*}
we observe that
\begin{align*} \eta_e &= -c(\mu+b) \int_t^\infty e^{-\mu^2s} q^b(s)\,ds\,e^{-(a^2-\mu^2)t} \\ &= -c(\mu+b) \int_t^\infty e^{\mu^2(t-s)} q^b(s)\,ds\,e^{-a^2t}. \end{align*}
Since
$e^{\mu^2(t-s)} \leq 1$ for
$s\geq t$, this implies
$\left| \eta_e/c(\mu+b)\right| \leq e^{-a^2t} \int_t^\infty q^b(s)\,ds \leq e^{-a^2t} \int_0^\infty q^b(s)\,ds = e^{-a^2t}/b$. The proof is now complete.
We conclude this section by giving remarks on corollary 1.3. It turns out that corollary 1.3 can be derived using the expression of the Dirichlet–Neumann map.
Remark 2.4. We consider the initial-boundary value problem for

Then, as in [Reference Caffarelli and Silvestre3], we obtain that

where $\gamma=1/(2-\alpha)$ with some constant
$c_\gamma \gt 0$, provided that α < 1. Indeed, as in [Reference Caffarelli and Silvestre3], let ψ be a solution of

Since the Laplace transform $\widehat{w}$ of w satisfies

we see, by scaling, that

Thus,

since $\eta(0)=0$. Thus,

If γ < 1, then $f_{1-\gamma}^0$ is integrable. As noted in [Reference Caffarelli and Silvestre3],
$\psi^\prime(0)$ exists (even for the degenerate case, i.e., α > 0) and
$\psi^\prime(0) \gt 0$. Thus,

at least for $\gamma=1/(2-\alpha) \lt 1$, i.e., α < 1.
Corollary 1.3 is easily derived by this result since $w_x=\partial_t^{1/2}w$ and
$w_x+bw=0$; note that
$c_{1/2}=1$.
Remark 2.5. The reader might be interested in how fractional partial differential equations like fractional diffusion equations are derived. We consider

where $f=f(y,t)$ is a given function. Then, by remark 2.4, the equation for
$\eta(y,t)=w(0,y,t)$ is formally obtained as

This type of equation is a kind of fractional diffusion that has been well-studied; see [Reference Kubica, Ryszewska and Yamamoto17, Reference Zacher32]. Here, we briefly recall only the well-posedness of its initial boundary value problem for Eq. (2.19) in a domain. In the framework of distributions, the well-posedness of its initial boundary value problems has been established in [Reference Sakamoto and Yamamoto25, Reference Zacher31] by using the Galerkin method. The unique existence of viscosity solutions for Eq. (2.19), including general nonlinear problems, has been established in [Reference Giga and Namba7, Reference Namba22] and also in [Reference Topp and Yangari28] for the whole space $\mathbf{R}^{n-1}$. The scope of equations these theories apply is different. However, it has been proved in [Reference Giga, Mitake and Sato6] that two notions of solutions (viscosity solution and distributional solution) agree for Eq. (2.19) when we consider the Dirichlet problem in a smooth bounded domain.
3. Convergence
The goal of this section is to prove lemma 1.4. For this purpose, we shall study a homogeneous version of Eq. (1.13) in a whole space $\mathbf{R}\times(0,\infty)$ and estimate the derivative of a solution as
$|y|\to\infty$. We begin with several estimates related to the heat semigroup in R. Let
$E(x,t)$ be the Gauss kernel, and write
$E_t(x)=E(x,t)$, i.e.,

and $E_t^a=e^{-a^2t}E_t$. Then,

solves

with initial data $f\in L^\infty(\mathbf{R})\cap C(\mathbf{R})$, i.e., f is bounded and continuous on R.
Proposition 3.1. Assume that $f\in C(\mathbf{R})\cap L^\infty(\mathbf{R})$ satisfies a decay condition

Then, for $T \gt \delta \gt 0$,

and

Proof. We may assume a = 0. We notice that

with some C independent of t and x, since $\partial_x E_t=-(x/2t)E_t$ and
$\sup_{y \gt 0}ye^{-y^2} \lt \infty$. Thus, it suffices to prove that

We divide

We notice that

with some constant C ʹ independent of $x,y\in\mathbf{R}$ and t < T. Since
$2|y|\leq|x|$ implies
$|x-y|\geq|x|-|y|\geq|x|/2$, we obtain that

For $II$, we proceed

We thus obtain Eq. (3.1). The proof is now complete.
We also need to estimate decay as $|x|\to\infty$ for

when $g(s)=h(s)\partial_x(1_{x \gt 0})$,
$h\in L^\infty(0,T)$. This quantity equals

Proposition 3.2. For $h\in L^\infty(0,T)$ and m > 0,

Proof. Again, we may assume that a = 0. Since

for $|x|\geq M$ and
$t\leq T$. Since
$t^{-1}e^{-M^2/(16t)}$ is integrable near t = 0, we see that

with some constants C 0, $C^\prime_0$ depending only on M and T. The proof is now complete.
We consider a homogeneous version Eq. (1.13) in $\mathbf{R}\times(0,\infty)$. Namely, we consider

Proposition 3.3. Let w be a bounded solution of Eq. (3.2) with initial data $w_0\in C(\mathbf{R})\cap L^\infty(\mathbf{R})$. Assume that

Then,

for any $\delta,T$ satisfying
$T \gt \delta \gt 0$.
Proof. We may assume $\tau_1=1$. By Duhamel’s formula, w is of the form

Applying propositions 3.1 and 3.2, we obtain desired results.
Proof. Proof of lemma 1.4
For a function f on $(-L,L)$, we decompose it into its odd and even parts, i.e.,

so that $f=f_\mathrm{odd}+f_\mathrm{even}$. By the structure of the equation,
$V_\mathrm{odd}^\varepsilon$ and
$V_\mathrm{even}^\varepsilon$ solve Eq. (1.13) separately.
At first glance, the locally uniform convergence follows from the maximum or comparison principles for a linear parabolic equation [Reference Protter and Weinberger24]. However, a direct application of the maximum principle is impossible since the domains of functions V ɛ and V are different. We first show the convergence where initial data are smooth.
For the odd part, the term $V\partial_y(1_{y \gt 0})$ does not affect since V = 0 at y = 0. Thus, Eq. (1.13) is reduced to

where $V_0^\varepsilon(y)=V^\varepsilon (y,0)$ for
$|y| \lt L/\varepsilon$. Let V ɛ be its solution.
We extend an odd function $V_0^\varepsilon$ to
$\left(L/\varepsilon,3L/\varepsilon\right)$ for x to be ‘even’ with respect to
$L/\varepsilon$, i.e.,

where $\widetilde{V_0^\varepsilon}$ is its extension. We extend
$\widetilde{V_0^\varepsilon}$ outside
$\left(-L/\varepsilon,3L/\varepsilon\right)$ so that the extension
$\overline{V_0^\varepsilon}$ is periodic in R with period
$4L/\varepsilon$. Since
$V_0^\varepsilon$ is even with respect to
$L/\varepsilon$,
$\left(\overline{V_0^\varepsilon}\right)_y \left(\pm L/\varepsilon\right)=0$ if
$V_{0y}^\varepsilon\left(\pm L/\varepsilon\right)=0$ and smooth. Solution V ɛ is the restriction on
$\left(-L/\varepsilon,L/\varepsilon\right)$ of a solution W ɛ of

Although the maximum principle implies

our assumption of the convergence $V_0^\varepsilon\to V_0$ does not guarantee
$\left\lVert \overline{V_0^\varepsilon} - V_0 \right\rVert_{L^\infty(\mathbf{R})}\to0$. We argue differently.
We approximate V 0 by $V_{0\delta}=V_0*\rho_\delta$, where ρδ is a symmetric mollifier. We also approximate
$V_0^\varepsilon$ by

We set $W_\delta^\varepsilon=W^\varepsilon*\rho_\delta$, where W ɛ is the solution of Eq. (3.4). Since Eq. (3.4) is of constant coefficients, this
$W_\delta^\varepsilon$ solves Eq. (3.4)1 with initial data
$V_{0\delta}^\varepsilon$. Let
$V_\delta^\varepsilon$ be the restriction of
$W_\delta^\varepsilon$ on
$(-L/\varepsilon,L/\varepsilon)$. It follows from the parity and periodic condition that this
$V_\delta^\varepsilon$ solves Eq. (3.3). We set
$V_\delta=V*\rho_\delta$ and observe that Vδ is the bounded solution of Eq. (3.4)1 with initial data
$V_{0\delta}$ since Eq. (3.4) is of constant coefficients. For fixed δ > 0, we observe that

for M > 0. Indeed, by the maximum principle

where $\|\cdot\|_{\infty,\varepsilon}$ is the sup norm on
$\left(-L/\varepsilon, L/\varepsilon\right)$. Here, the notation
$\partial_tV_{0\delta}^\varepsilon$ for a function
$V_{0\delta}^\varepsilon$ of y should be interpreted as

Because of the mollifier, the right-hand side is bounded by a constant multiple of $\left\lVert V_0^\varepsilon\right\rVert_{\infty,\varepsilon}$, which is uniformly bounded for ɛ < 1. By the Arzelà–Ascoli theorem and a diagonal argument,
$V_\delta^\varepsilon$ converges (locally uniformly in
$\mathbf{R}\times[0,\infty)$) to a bounded (weak) solution to Eq. (3.4) with initial data
$V_{0\delta}$ by taking a subsequence. Since V is bounded, by the uniqueness of the limit problem, the convergence is now full (without taking a subsequence). Note that we only invoke the locally uniform convergence of
$V_0^\varepsilon$ to V 0 other than the uniform bound on derivatives.
We note that

and observe that

where the norm is $\|\cdot\|$ taken in
$L^\infty(-M,M)$ for M > 0. By the maximum principle,

Since

and $\|\rho_\delta*f\|_\infty\leq\|f\|_\infty$, we see that

Thus,

Fixing δ > 0 and sending ɛ → 0, we observe that

since $V_\delta^\varepsilon\to V_\delta$ in
$L^\infty\left((-M,M)\times[0,T]\right)$. Sending
$\delta\downarrow0$, we obtain

since V 0 is uniformly continuous.
We next study the even part. The general strategy is the same but more involved than the odd part. For the even part, we first note that $V_\mathrm{even}^\varepsilon$ solves

where $V_0^\varepsilon(y)=V^\varepsilon(y,0)$ for
$y\in I_\varepsilon$. We suppress the word ‘even’ from now on. We shall approximate
$V_0^\varepsilon$ by a smoother function
$V_{0\delta}^\varepsilon$ and approximate
$V_0=\lim_{\varepsilon\to0}V_0^\varepsilon$ by a smoother function uniformly. There are many possible ways, and we rather like an abstract way. Let
$\text{BUC}(I_\varepsilon)$ denote the space of all bounded uniformly continuous functions in
$\overline{I_\varepsilon}$. It is a Banach space equipped with the norm

If ɛ = 0, then Iɛ should be interpreted as $(-\infty,0]$. Let A be the operator on
$\text{BUC}(I_0)$ defined by

with

A standard theory [Reference Lunardi18] implies that −A generates an analytic semigroup e −tA in $\text{BUC}(I_0)$. In particular,

with some constant C > 0 independent of time $t\in(0,1)$ and
$f\in BUC(I_0)$. For a function
$h\in BUC(I_\varepsilon)$, we extend it to
$\widetilde{h}$ so that
$\widetilde{h}(x)=h\left(-L/\varepsilon\right)$ for
$x \lt -L/\varepsilon$. For V 0, we set
$V_{0\delta}=e^{-\delta A}V_0$. For
$V_0^\varepsilon$, we tempt to set
$V_{0\delta}^\varepsilon=e^{-\delta A}\widetilde{V_0^\varepsilon}$. However, unfortunately,
$V_{0\delta}^\varepsilon$ does not satisfy the boundary condition at
$-L/\varepsilon$ although it satisfies
$\left(\partial_y V_{0\delta}^\varepsilon+bV_{0\delta}^\varepsilon\right)(0)=0$ and is C 2 (actually smooth). We set
$\sigma(x)=\rho_{\delta^\prime}*\left(1-|x|\right)_+$ for a fixed
$\delta^\prime \gt 0$ so that
$\sigma^\prime\left(1/2\right)=-\kappa$ with κ < 0. We set
$\sigma_{\delta^{\prime\prime}}(y)=\delta^{\prime\prime}\sigma\left(y/\delta^{\prime\prime}\right)$ for small
$\delta^{\prime\prime} \gt 0$. For a given
$h\in C^2(-\infty,0]$, we modify

where we take c so that $c\kappa=h_y\left(-L/\varepsilon\right)$,
$-L/\varepsilon-\nu=\delta^{\prime\prime}/2$. By this modification,

and $h_{\delta^{\prime\prime}}\to h$ in
$L^\infty(-\infty,0)$ as
$\delta^{\prime\prime}\to0$. We set

Since $V_{0\delta}^\varepsilon$ satisfies the boundary condition on the boundary of Iɛ and is smooth, we observe that
$\partial_t V_\delta^\varepsilon$ is continuous up to the boundary of
$I_\varepsilon\times[0,T)$, where
$V_\delta^\varepsilon$ denotes the solution of Eq. (3.5) with initial data
$V_{0\delta}^\varepsilon$. By the maximum principle,

where $\partial_tV_{0\delta}^\varepsilon$ for initial data
$V_{0\delta}^\varepsilon$ should be interpreted as in the proof for the odd part. The term involving b appears because of the Robin type boundary condition. As in the case for Eq. (3.4), by the Arzelà–Ascoli theorem and the uniqueness of the limit equation, we can prove that
$V_\delta^\varepsilon$ converges to Vδ locally uniformly in
$(-\infty,0]\times[0,\infty)$. Note that for a fixed δ > 0, the right-hand sides of Eq. (3.6) are uniformly bounded as ɛ → 0 since
$V_0^\varepsilon$ converges to V 0 uniformly. The comparison principle implies that

Thus,

where the norm $\|\cdot\|$ is taken in
$L^\infty(0,M)$ for M > 0. Taking the supremum in
$t\in(0,T)$ and sending ɛ → 0, we obtain that

since we know $\sup_{0 \lt t \lt T}\|V_\delta-V_\delta^\varepsilon\|(t)\to0$ as ɛ → 0. Sending δ → 0, we conclude that V ɛ converges to V locally uniformly in
$[0,\infty)\times[0,\infty)$.
Since we know that V ɛ is continuous up to y = 0 and t = 0, this gives the local uniform convergence of $V^\varepsilon(0,t)$ in
$[0,\infty)$.
So far, we do not invoke the assumption that $\left(V_0(x)-1\right)|x|\to0$ as
$|x|\to\infty$. We shall prove

We only discuss the case of even initial data since the odd case is easier. Since

it suffices to prove that

By our construction, the function $V_\delta(y,t)$ is

for w solving Eq. (3.2) with $w_0=V_0-1$. By proposition 3.3,

The difference $V_\delta-V_\delta^\varepsilon=u^\varepsilon$ satisfies

Let $\bar{u}^\varepsilon$ be the solution of

By the maximum principle,

Since we know that $\sup_{0 \lt t \lt T}\left|u_b^\varepsilon(x)\right|\to0$ as ɛ → 0, sending ɛ → 0 yields

We set $r=u^\varepsilon-\bar{u}^\varepsilon$ and observe that r satisfies

We set $\tilde{r}=e^{a^2t}r$ and observe that
$\tilde{r}$ satisfies

We take $M^\varepsilon=\sup_{0\leq t\leq T}e^{a^2t}\left|(\partial_y V_\delta)(-L/\varepsilon,0)\right|$ and observe that −yM and yM are the super- and subsolution, respectively. By the comparison principle,

This implies that

By Eq. (3.8), $M^\varepsilon L/\varepsilon\to0$ as ɛ → 0. We thus conclude that

Sending $\delta\downarrow0$, we obtain Eq. (3.7). The statement that
$\sup_{0\leq t\leq T}|V-1|(y)|y|\to0$ as
$|y|\to\infty$ follows from proposition 3.3. The proof is now complete.
If we only assume that the initial data for Eqs. (3.4) and (3.5) are bounded and Lipschitz, we have a similar estimate in Eq. (3.6) up to the first derivative of the solution. However, the estimate for the time derivative should be altered. Since we used such an estimate in lemma 2.2, we state it in the case of ɛ = 0 for the reader’s convenience.
Lemma 3.4. Let V be the bounded solution of Eq. (1.13) in $\mathbf{R}\times(0,\infty)$ with the bounded and Lipschitz continuous initial data V 0. Then, for each T > 0, there is a constant C depending only on a, b, and T such that

Proof. We give direct proof. We may assume that $\tau_1=1$. We set
$w=V-1$ and
$u=e^{a^2t}w$ to get

where we denote by x instead of y. We consider this equation with initial data $w_0=V_0-1$. It suffices by simple scaling
$u_\lambda(x,t)=u(\lambda x,\lambda^2t)$ to prove the desired estimate for some T independent of w 0.
Let $E(x,t)$ be the Gauss kernel as before. Then, the solution can be represented as

where

Since we can approximate a smooth w 0, establishing

with some positive constants C and T independent of w 0 suffices, assuming that $\partial_tu$ exists and is bounded in
$\mathbf{R}\times(0,T)$ for small T. By the maximum principle Eq. (3.6) and the corresponding estimate for the odd part, we know that

with c independent of w 0 and t > 0. We estimate $\partial_tu$ in Eq. (3.9). Since
$\|f\ast_xg\|_{L^\infty(\mathbf{R})}\le\|f\|_{L^1(\mathbf{R})}\|g\|_{L^\infty(\mathbf{R})}$ and

we easily see (cf. [Reference Giga, Giga and Saal5, Chapter 1]) that

with c ʹ independent of w 0. The second term of the right-hand side of Eq. (3.9) is more involved than the first term because h contains u. We observe that

Since $|E|\le(4\pi t)^{-1/2}$, it holds that

for $t\in(0,T)$.
Thus,

with some constants Cj ($j=1,2,3$). By Eqs. (3.9) and (3.12), we now observe that

with C 4 independent of w 0, u, and T. Applying estimate for $\|u\|_\infty$ in Eq. (3.11), we conclude Eq. (3.10) for sufficiently small T by absorbing
$C_2T^{1/2}\sup_{0 \lt t \lt T}\|t^{1/2}\partial_tu\|_{L^\infty(\mathbf{R})}$ on the right-hand side to the left.
4. Dirichlet condition for the total variation flow
In this section, we recall a notion of total variation flow for a given v and prove lemma 1.1. We consider

where $\beta\in C\left(I\times[0,T]\right)$ is a given non-negative function; here,
$I=(p_0,p_1)$ is an open interval and T > 0. If we impose the Dirichlet boundary condition

Eq. (4.1) with Eq. (4.2) should be interpreted as an L 2-gradient flow of a time-dependent total variation type energy

when $\int\beta|u_x|$ is a weighted total variation of β and
$\gamma u$ is a trace of u on
$\partial I$. We consider this energy in
$L^2(I)$ by
$\Phi^t(u)=\infty$ when
$\Phi^t(u)$ is not finite. It is clear that
$\Phi^t$ is convex in
$L^2(I)$. If β is spatially constant, it is well known that
$\Phi^t$ is also lower semicontinuous; see e.g. [Reference Andreu-Vaillo, Caselles and Mazón1]. The solution of Eq. (4.1) with Eq. (4.2) should be interpreted as the gradient flow of form

where $\partial\Phi^t$ denotes the subdifferential of
$\Phi^t$ in
$L^2(I)$, i.e.,

It is standard that Eq. (4.3) is uniquely solvable for given initial data $u_0\in L^2(I)$ if
$\Phi^t$ does not depend on time and is lower semicontinuous and convex on the Hilbert space
$L^2(I)$ (see, for instance, [Reference Brézis2, Reference Kōmura16]). It applies to the total variation flow case when β is a constant. In one-dimensional case,
$u\in BV(I)$ implies
$u\in L^\infty(I)$, so the subdifferential becomes

when $\Phi=\Phi^t$; see [Reference Andreu-Vaillo, Caselles and Mazón1, Proposition 5.10 and Lemma 5.13]Footnote 1. Here,
$(\beta z,u_x)$ denotes the Anzellotti pair [Reference Andreu-Vaillo, Caselles and Mazón1] and

Eq. (4.3) is

with $|z|\leq1$ in I and

with $-(\beta z)(p_i)(-1)^i(g-\gamma u)(p_i)=\beta|\gamma u-g|(p_i)$ for
$i=0,1$. We mimic this notion of the solution. A function
$u\in C\left([0,T),L^2(I)\right)$ is a solution to Eq. (4.1) with Eq. (4.2) if there is
$z\in L^\infty\left(I\times(0,T)\right)$ such that



Under this preparation, we shall prove lemma 1.1.
Proof. Proof of lemma 1.1
We set $p_0=-L$,
$p_1=L$ so that
$I=(-L,L)$. Since
$u_t^b=0$, Eq. (4.4) says that
$\beta z$ is a constant c. The condition Eq. (4.5) is equivalent to saying that
$|c|\leq\min\beta=\beta(0)$. Since

with $g(L)=b$,
$g(-L)=0$, Eq. (4.6) is equivalent to

In other words, c must be $\beta(0)$. Thus, the existence of z satisfying
$|z|\leq1$ is guaranteed if and only if

Eq. (4.6) is fulfilled with $u=u^b$ by taking
$g(-L)=0$ and
$g(L)=b$. Thus, ub is a stationary solution to Eq. (1.4) with Eq. (1.7).
5. Numerical experiment
In this section, we calculate the solution of Eqs. (1.10)–(1.12) with $v_0(x)=1-ce^{-a|x|/\varepsilon}$ and compare its value at x = 0 with an explicit solution of Eq. (2.16) whose explicit form is given in Eq. (2.18).
5.1. Numerical scheme
Since the initial function, v 0, is an even function, the original problem Eqs. (1.10)–(1.12) is reduced to

The computational region $[0,L]$ is divided into uniform mesh partitions:

The points x −1 and $x_{N+1}$ are needed to handle the Neumann boundary conditions.
The approximation of v at $x=x_i$ is written as vi. The central finite difference approximates the Laplace operator, and the time derivative is approximated by the backward difference, yielding the following linear system:

where $\hat{v}_i$ is the value known at the current time and vi is the value to be found at the next time. The Neumann boundary conditions at x = 0 and x = L are approximated as

by the central finite differences.
5.2. Results
Some results are shown for different values of c for parameters

The results of the numerical experiments are summarized in figure 2. In (a), (b), and (c), the horizontal axis represents time, and the vertical axis represents the value at the origin. As ɛ is decreased, the numerical solution converges to the exact solution to the extent that the exact and numerical solutions overlap. Indeed, the table of $L^\infty$-errors for different values of ɛ and c is shown in (d). The errors for
$\varepsilon=2^{-3}$ are of order 10−5, indicating that the solution for small ɛ is an excellent approximation to the solution of the fractional time differential equation obtained as the singular limit.

Figure 2. Results of numerical experiments: (a) c = 0, (b) $c=1/4$, (c) c = 2, (d) table of
$L^\infty$-errors for different ɛ and c values.
Acknowledgements
The authors are grateful to the anonymous reviewer for his/her valuable expository suggestions. The work of the first author was partly supported by JSPS KAKENHI Grant Numbers JP24K00531, JP19H00639, and JP20K20342, and by Arithmer Inc., Daikin Industries, Ltd. and Ebara Corporation through collaborative grants. The work of the third author was partly supported by JSPS KAKENHI Grant Number JP20K20342. The work of the fifth author was partly supported by JSPS KAKENHI Grant Numbers JP22K03425, JP22K18677, and 23H00086.