The stability condition for the FDTD of the optical diffusion equations

Abstract

The finite-difference time-domain (FDTD) method, which is often used to solve Maxwell’s equations, has been proposed to solve the optical diffusion equations to simulate light propagation in turbid media such as biological tissue. Finite-difference methods can be numerically unstable and calculation can diverge unless a stability condition is satisfied. Until now, the stability condition for FDTD of the optical diffusion equations has not been derived. In this paper, we derive a stability condition for FDTD of the optical diffusion equations using the von Neumann stability analysis and verify the condition numerically.

Appendices

Appendix 1: Proof of inequalities (10) and (11)

In this appendix, we prove the inequalities (10) and (11). We rewrite these inequalities below:

$$\begin{aligned} | \gamma _\pm |< 1 \quad {\text {as}} \quad \upsilon \varDelta t < 1/\sqrt{P} \\ | \gamma _ -| \ge 1 \quad {\text {as}} \quad \upsilon \varDelta t \ge 1/\sqrt{P}. \end{aligned}$$

We do not need to know the value $| \gamma _ +|$ as $\upsilon \varDelta t > 1/\sqrt{P}$. Whether the value $| \gamma _ +|$ is less than one or not; this FDTD is unstable when $\upsilon \varDelta t > 1/\sqrt{P}$, because the other eigenvalue $| \gamma _ -|$ is larger than one as we prove below. Substituting Eqs. (3) and (6) into Eq. (9):

$$\begin{aligned} {\gamma _ \pm } = {g_1} \pm {g_2}, \end{aligned}$$

(19)

where

$$\begin{aligned} \begin{aligned} g_1&=\frac{4-(\mu _a \mu _{tr} + 8P) \upsilon ^2 \varDelta t^2}{(2+\upsilon \varDelta t \mu _a)(2+\upsilon \varDelta t \mu _{tr})}\\ g_2&= \frac{2\upsilon \varDelta t\sqrt{4P\left( \mu _a\mu _{tr} + 4P \right) {\upsilon ^2}\varDelta t^2 -16P+{(\mu _{tr}-\mu _a)}^2} }{(2+\upsilon \varDelta t \mu _a)(2+\upsilon \varDelta t\mu _{tr})}. \end{aligned} \end{aligned}$$

(20)

Now, substituting $\upsilon \varDelta t = 1/\sqrt{P}$ into Eq. (19), we get:

$$\begin{aligned} \gamma _+|_{\upsilon \varDelta t= 1/\sqrt{P} }&= -\frac{ \left( 2 - P^{-1/2} \mu _a \right) \left( 2 - P^{-1/2} \mu _{tr} \right) }{ \left( 2 + P^{-1/2} \mu _a \right) \left( 2 + P^{-1/2} \mu _{tr} \right) } \end{aligned}$$

(21a)

$$\begin{aligned} \gamma _- |_{\upsilon \varDelta t = 1/\sqrt{P} }&= - 1. \end{aligned}$$

(21b)

Around $\upsilon \varDelta t = 1/\sqrt{P}$, $\gamma _-$ can be larger than one, and hence, $\upsilon \varDelta t = 1/\sqrt{P}$ is the candidate for stability bound. Considering the function $\left| {\frac{{2 - x}}{{2 + x}}} \right| \,\,\,{\text {as}}\,\,x > 0$, we get $\left| {{{\left. {\,\,{\gamma _ + }} \right| }_{\upsilon \varDelta t = 1/\sqrt{P} }}} \right| < 1$.

The square root in $g_2$ is zero if

$$\begin{aligned} \upsilon \varDelta t = \sqrt{\frac{{16P - {{\left( {{\mu _{tr}} - {\mu _a}} \right) }^2}}}{{4P\left( {{\mu _a}{\mu _{tr}} + 4P} \right) }}}. \end{aligned}$$

Then, we define $\beta _r$ by:

$$\begin{aligned} {\beta _r}^2 = \frac{{16P - {{\left( {{\mu _{tr}} - {\mu _a}} \right) }^2}}}{{4P\left( {{\mu _a}{\mu _{tr}} + 4P} \right) }}. \end{aligned}$$

(22)

Since we know that ${\left( {{\mu _{tr}} - {\mu _a}} \right) ^2} > 0$ and $\mu _a$, $\mu _{tr}>0$:

$$\begin{aligned} {\beta _r}^2 = \frac{{16P - {{\left( {{\mu _{tr}} - {\mu _a}} \right) }^2}}}{{4P\left( {{\mu _a}{\mu _{tr}} + 4P} \right) }} < \frac{{16}}{{{4^2}P}} = {\left( {\frac{1}{{\sqrt{P} }}} \right) ^2}. \end{aligned}$$

(23)

Therefore, the region over which $g_2$ is not real is smaller than $\upsilon \varDelta t = 1/\sqrt{P}$. Figure 3 depicts this relation. Note that if $16P \le {\left( {{\mu _{tr}} - {\mu _a}} \right) ^2}$, then $g_2$ is real for any $\upsilon \varDelta t>0$.

Fig. 3

The schematic behavior of $\gamma _\pm$

If we can prove the states (A1) and (A2) below, we show that $|\gamma _{\pm }|\le 1$ as $\upsilon \varDelta t \le 1/\sqrt{P}$, and (A3) leads $|\gamma _{\pm }|>1$ as $\upsilon \varDelta t > 1/\sqrt{P}$. We can, therefore, prove what we set out to prove in this appendix.

• (A1) $|\gamma _\pm | < 1\quad$ as $\quad 0<\upsilon \varDelta t \le \beta _r$.

• (A2) $| \gamma _+|< 1,| \gamma _-| < 1 \quad$ as $\quad \upsilon \varDelta t \in \left( \beta _r,1/\sqrt{P}\right)$.

• (A3) $\gamma _{-}=-1$ as $\upsilon \varDelta t = 1/\sqrt{P}$, and $\gamma _-$ monotonically decreases as $\upsilon \varDelta t > {\beta _r}$

Proof of (A1)

We can obtain $|\gamma _\pm |^2$ from Eq. (19):

$$\begin{aligned} {\left| \gamma _\pm \right| ^2} = g_1^2 + g_2^2 = \frac{{{{\left( {2 - \upsilon \varDelta t{\mu _a}} \right) }^2}{{\left( {2 - \upsilon \varDelta t{\mu _{tr}}} \right) }^2}}}{{{{\left( {2 + \upsilon \varDelta t{\mu _a}} \right) }^2}{{\left( {2 + \upsilon \varDelta t{\mu _{tr}}} \right) }^2}}}; \end{aligned}$$

(24)

function $\left| {\frac{{2 - x}}{{2 + x}}} \right|$ is less than unity as $x>0$, and therefore, as we have considered in Eq. (21a), we get $\left| {{\gamma _ \pm }} \right| < 1$.

Proof of (A2)

As for ${\gamma _ + } = {g_1} + {g_2}$, $\left| {{\gamma _ + }} \right| = \left| {{g_1} + {g_2}} \right| \le + 1$ is equivalent to two inequalities $- 1 - {g_1} \le {g_2}$ and ${g_2} \le 1 - {g_1}$. First, we show $- 1 - {g_1} \le {g_2}$. Consider the region $\upsilon \varDelta t < 1/\sqrt{P}$, so that $1 - {\upsilon ^2}\varDelta {t^2}P > 0$, and

$$\begin{aligned} - 1 - {g_1} = - \frac{{8\left( {1 - P{\upsilon ^2}\varDelta {t^2}} \right) + 2\upsilon \varDelta t\left( {{\mu _a} + {\mu _{tr}}} \right) }}{{\left( {2 + \upsilon \varDelta t{\mu _a}} \right) \left( {2 + \upsilon \varDelta t{\mu _{tr}}} \right) }} < 0. \end{aligned}$$

(25)

Because apparently $g_2>0$, so that $- 1 - {g_1} \le {g_2}$. Next, we want to show that $- 1 - {g_1} \le {g_2}$, and so we can calculate:

$$\begin{aligned} {\left( {1 - {g_1}} \right) ^2} - {g_2}^2 = \frac{{4{\upsilon ^2}\varDelta {t^2}\left( {4P + {\mu _a}{\mu _{tr}}} \right) }}{{\left( {2 + \upsilon \varDelta t{\mu _a}} \right) \left( {2 + \upsilon \varDelta t{\mu _{tr}}} \right) }} > 0. \end{aligned}$$

(26)

Therefore ${g_2}^2 \le {\left( {1 - {g_1}} \right) ^2}$. $g_2$ is positive and

$$\begin{aligned} 1 - {g_1} = \frac{{2\upsilon \varDelta t\left( {4P\upsilon \varDelta t + {\mu _{tr}} + {\mu _a} + \upsilon \varDelta t{\mu _a}{\mu _{tr}}} \right) }}{{\left( {2 + \upsilon \varDelta t{\mu _a}} \right) \left( {2 + \upsilon \varDelta t{\mu _{tr}}} \right) }} > 0, \end{aligned}$$

(27)

because all the parameters in this equation are positive. Since the square function is an increasing monotonic function in the positive domain, we have ${g_2} \le 1 - {g_1}$. We have shown that $\left| {{\gamma _ + }} \right| \le 1$ as $\upsilon \varDelta t \in \left( {{\beta _r},1/\sqrt{P} } \right)$. Although it is omitted here, the bound on $\gamma _-$ can be shown in the same way.

Proof of (A3)

We have shown that $\gamma _-=-1$ as $\upsilon \varDelta t = 1/ \sqrt{P}$ in Eq. (21b), so it only remains to show that $\gamma _-$ monotonically decreases as for $\upsilon \varDelta t$.

Differentiating ${g_2}^2$ with respect to $\upsilon \varDelta t > 0$:

$$\begin{aligned} \begin{aligned} \frac{\partial {g_2}^2}{{\partial \left( {\upsilon \varDelta t} \right) }}&= \frac{{32P\upsilon \varDelta t\left( {4P + {\mu _a}{\mu _{tr}}} \right) }}{{{{\left( {2 + \upsilon \varDelta t \mu _a} \right) }^3}{{\left( {2 + \upsilon \varDelta t{\mu _{tr}}} \right) }^3}}} \left( {4\left( \upsilon ^2 \varDelta t^2 - {\beta _r}^2 \right) + 4{\upsilon ^2} \varDelta t^2} + 2 \upsilon ^2 \varDelta t^3 \left( {{\mu _{tr}} + \mu _a } \right) { + {\upsilon ^2}\varDelta {t^2}{\beta _r}^2{\mu _{tr}}{\mu _a}} \right) \\&> 0. \end{aligned} \end{aligned}$$

(28)

Therefore, ${g_2}^2$ is monotonically increasing. Because $g_2>0$, $g_2$ is also monotonically increasing.

Differentiating $g_1$:

$$\begin{aligned} \frac{\partial {g_1}}{{\partial \left( {\upsilon \varDelta t} \right) }} = - \frac{{\left( {{\upsilon ^2}\varDelta {t^2}{\mu _a}{\mu _{tr}} + 4 + 8P{\upsilon ^2}\varDelta {t^2}} \right) \left( {{\mu _a} + {\mu _{tr}}} \right) }}{{{{\left( {2 + \upsilon \varDelta t{\mu _a}} \right) }^2}{{\left( {2 + \upsilon \varDelta t{\mu _{tr}}} \right) }^2}/2}} - \frac{{8\upsilon \varDelta t{\mu _{tr}}{\mu _a} + 32P\upsilon \varDelta t}}{{{{\left( {2 + \upsilon \varDelta t{\mu _a}} \right) }^2}{{\left( {2 + \upsilon \varDelta t{\mu _{tr}}} \right) }^2}/2}}<0 \end{aligned}$$

(29)

shows that $g_1$ is monotonically decreasing and $\gamma _-=g_1-g_2$ is also monotonically decreasing.

Appendix 2: The other results of the simulation of the different initial distribution

As shown in Sect. 3, the stability condition depends on the dimensions of the field; $\lambda =\sqrt{3}$ in 1D, and $\lambda =\sqrt{3/2}$ in 2D. We, therefore, include the results of the simulation in these dimensions. We use the same simulation parameters as in Sect. 4 except that the initial conditions, as stated in Eq. (18) for the 3-D case, become Eqs. (30a) and (30b) for 1-D and 2-D, respectively:

$$\begin{aligned} \varPhi _0^{1{\text {D}}}\left( {x,y,z} \right)&= 1 - \cos \left( {2\pi \,x} \over {p_x} \right) \end{aligned}$$

(30a)

$$\begin{aligned} \varPhi _0^{2{\text {D}}}\left( {x,y,z} \right)&= 1 - \frac{1}{2}\cos \left( {{{2\pi \,x} \over {{p_x}}}} \right) - \frac{1}{2}\cos \left( {{2\pi \,y} \over {{p_y}}} \right) , \end{aligned}$$

(30b)

where $p_x$ and $p_y$ are the period of the initial distribution along each axis and we put ${p_x} = {p_y} = 20\varDelta x$. As in the 3D simulation, we choose the Courant numbers to be 1% larger and 1% lower than the stability bound. Figure 4 shows that the results again support the boundary condition expressed by Eq. (17).

Fig. 4

Simulation results around stability bound in 1D and 2D. The peak-to-valley (P–V) of $\varPhi$ in each time-step a with the 1D initial distribution expressed by Eq. (30a) and b with the 2D initial distribution expressed by Eq. (30b). Open square: $\lambda = 0.99\sqrt{3/d}$, diagonal cross: $\lambda = 1.01\sqrt{3/d}$, and dashed line is $C |\gamma _-|^n$. ${\gamma _-} \cong - 1.0872$ when $\lambda =1.01 \sqrt{3}$ in 1D simulation, and ${\gamma _ - } \cong - 1.1162$ when $\lambda =1.01 \sqrt{3/2}$ in 2D-simulation. The dotted line indicates the theoretical amplitude

Appendix 3: The theoretical curve in the single-mode distribution in the first numerical simulation

In this appendix, we calculate the time dependence of the peak-to-valley when the initial distribution is expressed by Eq. (18).

Because the diffusion equations are linear, we only need to consider these terms:

$$\begin{aligned} \begin{aligned} \phi \left( {x,y,z;t} \right) = {b_\phi }\exp \left( {{{{2\pi \mathrm {i}\,x} \over {{p_x}}}}} \right) \\ {{\varvec{J}}}\left( {x,y,z;t} \right) = \left( {\begin{array}{c} {{b_j}\exp \left( {{{{2\pi \mathrm {i}\,x} \over {{p_x}}}}} \right) }\\ 0\\ 0 \end{array}} \right) , \end{aligned} \end{aligned}$$

(31)

where $b_\phi$ and $b_j$ depend only on time. $p_x$ is a spatial frequency.

Substitute the above equations into:

$$\begin{aligned} \begin{aligned} \frac{1}{\upsilon }\frac{{\partial {{\varvec{J}}}\left( {{{\varvec{r}}},t} \right) }}{{\partial t}} + \frac{1}{3}\nabla \phi \left( {{{\varvec{r}}},t} \right) + {\mu _{tr}}{{\varvec{J}}}\left( {{{\varvec{r}}},t} \right)&= 0\\ \frac{1}{\upsilon }\frac{{\partial \phi \left( {{{\varvec{r}}},t} \right) }}{{\partial t}} + \nabla \cdot {{\varvec{J}}}\left( {{{\varvec{r}}},t} \right) + {\mu _a}\phi \left( {{{\varvec{r}}},t} \right)&= 0. \end{aligned} \end{aligned}$$

We then obtain:

$$\begin{aligned} \frac{\partial }{{\partial t}}\left( {\begin{array}{*{20}{c}} {{b_\phi }}\\ {{b_j}} \end{array}} \right) = - \upsilon \left( {\begin{array}{*{20}{c}} {{\mu _a}}&{}{{\textstyle {{2\pi \mathrm {i}} \over {{p_x}}}}}\\ {{\textstyle {{2\pi \mathrm {i}} \over {3{p_x}}}}}&{}{{\mu _{tr}}} \end{array}} \right) \left( {\begin{array}{*{20}{c}} {{b_\phi }}\\ {{b_j}}. \end{array}} \right) . \end{aligned}$$

This differential equation is solved by orthogonalizing the matrix in this equation:

$$\begin{aligned} \frac{\partial }{{\partial t}}\left( {\begin{array}{*{20}{c}} {{b_\phi }}\\ {{b_j}} \end{array}} \right) = - \upsilon {P^{ - 1}}\Lambda P\left( {\begin{array}{*{20}{c}} {{b_\phi }}\\ {{b_j}} \end{array}} \right) , \end{aligned}$$

where

$$\begin{aligned} \Lambda= & {} \left( {\begin{array}{*{20}{c}} {{\lambda _ - }}&{}0\\ 0&{}{{\lambda _ + }} \end{array}} \right) \\ P= & {} \left( {\begin{array}{*{20}{c}} {{p_ + }}&{}1\\ {{p_ - }}&{}1 \end{array}} \right) \\ {\lambda _ \pm }= & {} \frac{1}{2}\left( {{\mu _a} + {\mu _{tr}}} \right) \pm \frac{{\sqrt{9p_x^2{{\left( {{\mu _a} - {\mu _{tr}}} \right) }^2} - 48{\pi ^2}} }}{{6{p_x}}}\\ {p_ \pm }= & {} \frac{i}{{4\pi }}\left[ {3{p_x}\left( {{\mu _{tr}} - {\mu _a}} \right) \pm \sqrt{9p_x^2{{\left( {{\mu _a} - {\mu _{tr}}} \right) }^2} - 48{\pi ^2}} } \right] . \end{aligned}$$

We then obtain the solution with the initial values $b_\phi =1$ and $b_j=b_j^0$ at $t=0$:

$$\begin{aligned} \left( {\begin{array}{*{20}{c}} {{b_\phi }}\\ {{b_j}} \end{array}} \right) = {P^{ - 1}}\exp \left( { - \upsilon \Lambda t} \right) P\left( {\begin{array}{*{20}{c}} {{b_\phi }^0}\\ {{b_j}^0} \end{array}} \right) . \end{aligned}$$

In the simulation described in the article, we can use $b_j^0$, and from the (1,1)-element of ${P^{ - 1}}\exp \left( { - \upsilon \Lambda t} \right) P$, we have the solution of the differential equation:

$$\begin{aligned} \phi \left( {x,y,z;t} \right)&= {b_\phi }^0 \exp \left( \frac{2\pi \mathrm {i}\,x}{p_x} \right) \exp \left[ -\upsilon t\left( \frac{\mu _a + \mu _{tr}}{2} + \frac{\eta }{6 p_x} \right) \right] \\&\quad \times \left\{ \left( \frac{1}{2}-\frac{3 p_x (\mu _a -\mu _{tr})}{2\eta } \right) \exp \left( \frac{\upsilon \eta }{3 p_x} t \right) + \frac{1}{2} + \frac{3{p_x}\left( {\mu _a} - {\mu _{tr}} \right) }{2\eta } \right\} , \end{aligned}$$

where $\eta = \sqrt{9p_x^2{{\left( {{\mu _a} - {\mu _{tr}}} \right) }^2} - 48{\pi ^2}}$

Finally, we obtain the theoretical time dependence of peak-to-valley as:

$$\begin{aligned} 2{b_\phi }^0\exp \left[ {-\upsilon t\left( {\frac{\mu _a + \mu _{tr}}{2} + \frac{\eta }{6 p_x}} \right) } \right] \left\{ {\left( {\frac{1}{2} - \frac{{3{p_x}\left( {{\mu _a} - {\mu _{tr}}} \right) }}{{2\eta }}} \right) \exp \left( {\frac{{\upsilon \eta }}{{3{p_x}}}t} \right) + \frac{1}{2} + \frac{{3{p_x}\left( {{\mu _a} - {\mu _{tr}}} \right) }}{{2\eta }}} \right\} . \end{aligned}$$

When $p_x>0$, this value is exponentially increasing as t increases, and this expression shows that the increasing rate becomes higher as the spatial period $p_x$ becomes smaller. This is the reason why the error forms the checkerboard pattern whose spatial frequency is the highest in the 2D image, as shown in Fig. 1b or Fig. 2.