Transparent boundary conditions for wave propagation in fractal trees: convolution quadrature approach

Abstract

In this work we propose high-order transparent boundary conditions for the weighted wave equation on a fractal tree, with an application to the modeling of sound propagation in a human lung. This article follows the recent work (Joly et al. in Netw Heterog Media 14(2):205–264, 2019), dedicated to the mathematical analysis of the corresponding problem and the construction of low-order absorbing boundary conditions. The method proposed in this article consists in constructing the exact (transparent) boundary conditions for the semi-discretized problem, in the spirit of the convolution quadrature method developed by Ch. Lubich. We analyze the stability and convergence of the method, and propose an efficient algorithm for its implementation. The exposition is concluded with numerical experiments.

Appendices

Proof of Theorem 2.4

It remains to prove the upper bound on $\varvec{\Lambda }_{{\mathfrak {a}}}(\omega )$. Without loss of generality, we will show it for $\varvec{\Lambda }_{{\mathfrak {n}}}(\omega )$. First, $\varvec{\Lambda }_{{\mathfrak {n}}}(\omega )$ can be defined via the solution of the frequency-domain problem:

$$\begin{aligned} \varvec{\Lambda }_{{\mathfrak {n}}}(\omega )=-\partial _s \lambda (M^*),\quad \omega \in {\mathbb {C}}^+, \end{aligned}$$

(101)

where $\lambda \in H _\mu ^1({\mathcal {T}})$ solves the boundary-value problem:

$$\begin{aligned} \omega ^2\int \limits _{{\mathcal {T}}} \mu (s)\, \lambda \, {\overline{v}}-\int \limits _{{\mathcal {T}}} \mu (s)\, \partial _s\lambda \, \partial _s {\overline{v}}=0, \quad \text {for all }v\in V_{{\mathfrak {n}}}, \quad \lambda (M^*)=1. \end{aligned}$$

(102)

Let us define $\Vert v\Vert _{\omega }:=\int \limits _{{\mathcal {T}}}\mu \left( |\partial _s v|^2+|\omega v|^2\right)$. We proceed as follows:

• first prove the bound $|\varvec{\Lambda }_{{\mathfrak {n}}}(\omega )|^2$ by the energy of the solution (notice that $\lambda (M^*)=1$):

  $$\begin{aligned} \left| \varvec{\Lambda }_{{\mathfrak {n}}}(\omega )\right| ^2 \le |\omega |^2+C_0(1+{\text {Im}}{\omega })\Vert \lambda \Vert _{\omega }^2, \quad C_0>0. \end{aligned}$$

  (103)

• next show that the energy of the solution is bounded by $\frac{1}{2}|\varvec{\Lambda }_{{\mathfrak {n}}}(\omega )|^2$, with $C_0$ as above:

  $$\begin{aligned} \begin{aligned} C_0(1+{\text {Im}}{\omega })\Vert \lambda \Vert _{\omega }^2&\le \frac{1}{2}|\varvec{\Lambda }_{{\mathfrak {n}}}(\omega )|^2+C_1\max (1, ({\text {Im}}{\omega })^{-2})|\omega |^2, \quad C_1>0. \end{aligned} \end{aligned}$$

  (104)

• combine (103) and (104) to obtain the desired bound:

  $$\begin{aligned} \left| \varvec{\Lambda }_{{\mathfrak {n}}}(\omega )\right| ^2 \le C\max (1, ({\text {Im}}{\omega })^{-2})|\omega |^2. \end{aligned}$$

Proof of the bound (103) Let $v_0(s)=\chi (s)\partial _s \lambda$, where $\chi \in C^1({\mathcal {T}}; {\mathbb {R}})$, ${\text {supp}}\chi (s)\subseteq \varSigma _{0,0}$, $\chi (M^*)=1$ and $\chi (M_{0,0})=0$. The weak formulation (102) implies that $\lambda$ satisfies

$$\begin{aligned} \partial _s^2 \lambda +\omega ^2 \lambda =0 \text { on }\varSigma _{0,0}. \end{aligned}$$

Testing the above with $v_0(s)$, we obtain the following identity on the edge $\varSigma _{0,0}$, parametrized by $s\in [0,1]$ (recall that we work with the reference tree, and thus the length of $\varSigma _{0,0}$ is 1):

$$\begin{aligned} I_1+I_2=0, \quad \text { where } I_1=\int \limits _0^1 \partial _s^2 \lambda \, \chi (s)\, \partial _s {\overline{\lambda }}ds, \quad I_2=\omega ^2 \int \limits _0^1 \lambda \, \chi (s)\, \partial _s {\overline{\lambda }}ds. \end{aligned}$$

(105)

Let $\omega =\omega _r+i\omega _i$, $\omega _r\in {\mathbb {R}}$ and $\omega _i>0$. Let us consider the real part of the above:

$$\begin{aligned} {\text {Re}}{I}_1=\frac{1}{2}\int \limits _0^1\frac{d}{ds}|\partial _s \lambda |^2\chi (s)ds=-\frac{1}{2}\left| \varvec{\Lambda }_{{\mathfrak {n}}}(\omega )\right| ^2-\frac{1}{2}\int \limits _0^1\chi '(s)|\partial _s \lambda |^2ds, \end{aligned}$$

(106)

where in the last identity we used $\chi (0)=1$ and $\chi (1)=0$. Combining (105), (106), we deduce

$$\begin{aligned} |\varvec{\Lambda }_{{\mathfrak {n}}}(\omega )|^2\le 2|{\text {Re}}{I}_2|+c_1\int \limits _0^1|\partial _s \lambda |^2 ds, \quad c_1>0. \end{aligned}$$

(107)

Similarly,

$$\begin{aligned} {\text {Re}}{I}_2&=\frac{1}{2}{\text {Re}}{\omega ^2}\int \limits _0^{1}\chi (s) \frac{d}{ds}|\lambda |^2ds-{\text {Im}}\omega ^2 \int \limits _0^1 \chi (s){\text {Im}}(\lambda \partial _s{\overline{\lambda }})ds\\&=-\frac{1}{2}(\omega ^2_r-\omega _i^2)-\frac{1}{2}(\omega ^2_r-\omega _i^2)\int \limits _0^1\chi '(s)|\lambda |^2ds-2\omega _i \omega _r \int \limits _0^1 \chi (s){\text {Im}}(\lambda \partial _s{\overline{\lambda }})ds. \end{aligned}$$

where we used $\chi (0)=1$, $\chi (1)=0$ and $\lambda (0)=1$. Applying to the last integral the Young inequality we obtain the following bound, with $c_2, c_3>0$,

$$\begin{aligned} |{\text {Re}}{I}_2|\le \frac{1}{2}|\omega |^2+c_2|\omega |^2\int \limits _0^1|\lambda |^2ds+c_3\omega _i\left( |\omega _r|^2\int \limits _0^1|\lambda |^2ds+\int \limits _0^1|\partial _s \lambda |^2 ds\right) . \end{aligned}$$

(108)

Inserting (108) into (107) we prove (103).

Proof of the bound (104) Testing the Helmholtz equation corresponding to (102) with $\omega \lambda (s)$ and integrating by parts we obtain the following identity (recall that $\lambda (M^*)=1$):

$$\begin{aligned} {\overline{\omega }}\varvec{\Lambda }_{{\mathfrak {n}}}(\omega )={\overline{\omega }}\int \limits _{{\mathcal {T}}}\mu |\partial _s \lambda |^2-|\omega |^2\omega \int \limits _{{\mathcal {T}}}\mu |\lambda |^2. \end{aligned}$$

Taking the imaginary part of the above results in

$$\begin{aligned} {\text {Im}}\left( {\overline{\omega }}\varvec{\Lambda }_{{\mathfrak {n}}}(\omega )\right) =-\omega _i\left( \int \limits _{{\mathcal {T}}}\mu |\partial _s \lambda |^2+|\omega |^2\int \limits _{{\mathcal {T}}}\mu |\lambda |^2\right) =-\omega _i\Vert \lambda \Vert _{\omega }^2. \end{aligned}$$

Multiplying both sides of the above by $-C_0(1+\omega _i)\omega _i^{-1}$, with $C_0$ is as in (103), we obtain

$$\begin{aligned} -C_0\left( \omega _i^{-1}+1\right) {\text {Im}}\left( {\overline{\omega }}\varvec{\Lambda }_{{\mathfrak {n}}}(\omega )\right) =C_0(1+\omega _i)\Vert \lambda \Vert _{\omega }^2. \end{aligned}$$

(109)

It suffices to notice that the left hand side in the above equality is bounded:

$$\begin{aligned}&\left| -C_0\left( \omega _i^{-1}+1\right) {\text {Im}}\left( {\overline{\omega }}\varvec{\Lambda }_{{\mathfrak {n}}}(\omega )\right) \right| \\&\quad \le C_0\left( \omega _i^{-1}+1\right) |\omega ||\varvec{\Lambda }_{{\mathfrak {n}}}(\omega )|\le \frac{1}{2}|\varvec{\Lambda }_{{\mathfrak {n}}}(\omega )|^2+ \frac{C_0^2}{2}\left( \omega _i^{-1}+1\right) ^2|\omega |^2, \end{aligned}$$

where we used the Young inequality. In the above we bound further $C_0(\omega _i^{-1}+1)\le 2\max (1, \omega _i^{-1})$. Inserting the bound into (109) gives

$$\begin{aligned}&C_0(1+\omega _i)\Vert \lambda \Vert _{\omega }^2\le \frac{1}{2}|\varvec{\Lambda }_{{\mathfrak {n}}}(\omega )|^2+2C_0^2\max \left( 1,\omega _i^{-2}\right) |\omega |^2, \end{aligned}$$

i.e. (104). Combining (103) and (104) proves the statement of the theorem.

Proof of Lemma 4.1

We first show (86). By definition, $\tan \omega =i\frac{1-z}{1+z}$, with $z=\mathrm {e}^{2i\omega }, \, \omega =\omega _r+i\omega _i$. Then,

$$\begin{aligned} -{\text {Im}}\left( \tan \omega \right) ^{-1}&={\text {Re}}\frac{1+z}{1-z}={\text {Re}}\frac{(1+z)(1-{\overline{z}})}{|1-z|^2}=\frac{1-|z|^2}{1+|z|^2-2{\text {Re}}z}\ge \frac{1-|z|^2}{(1+|z|)^2}\nonumber \\&=\frac{1-|z|}{1+|z|}=\frac{1-\mathrm {e}^{-2\omega _i}}{1+\mathrm {e}^{-2\omega _i}}\ge \displaystyle \left\{ \begin{array}{ll} \frac{\mathrm {e}^{2\omega _i}-1}{\mathrm {e}^{2\omega _i}+1}\ge \frac{2\omega _i}{\mathrm {e}^2+1}, &{} \text { if } 0<\omega _i\le 1,\\ &{}\\ \frac{1-\mathrm {e}^{-2}}{1+\mathrm {e}^{-2}}, &{} \text { if }\omega _i>1, \end{array} \right. \end{aligned}$$

(110)

hence the bound (86). Let us show (87). After straightforward computations,

$$\begin{aligned} \left| 1-i\left( \tan \omega \right) ^{-1}\right|&=\frac{2|z|}{|1-z|}\le \frac{2|z|}{|1-|z||}= \frac{2\mathrm {e}^{-2\omega _i}}{1-\mathrm {e}^{-\omega _i}}\le C\max \left( 1, \omega _i^{-1}\right) \mathrm {e}^{-2\omega _i}, \end{aligned}$$

where the last bound follows by noticing that, for $\omega _i>0$,

$$\begin{aligned} 1-\mathrm {e}^{-\omega _i}\ge \left\{ \begin{array}{ll} 1-\mathrm {e}^{-1}, &{} \text {if }\omega _i\ge 1, \\ \mathrm {e}^{-1}\omega _i, &{} \text { if }\omega _i <1 \end{array}\right. \quad \ge c\min \left( 1, \omega _i\right) , \quad c>0. \end{aligned}$$

(111)

Proof of Proposition 4.3

To prove Proposition 4.3, we need the following auxiliary result.

Lemma C.1

Let $0<\rho <1$, $\varepsilon >0$, and $\lambda _{s,n}^{\varDelta t, \varepsilon }, \, n=0, \ldots , N_t$ be given by (93), with $N\ge N_t+1$, where $\max \limits _{k}\left| \varvec{\Lambda }^{s,\varepsilon }\left( \omega _k\right) -\varvec{\Lambda }^s\left( \omega _k\right) \right| <\varepsilon$. Then

$$\begin{aligned} \begin{aligned}&\max \limits _{n=0, \ldots , N_t}\left| \lambda _{s,n}^{\varDelta t, \varepsilon }-\lambda _{s,n}^{\varDelta t}\right| <\rho ^{-N_t}\varepsilon +\rho ^NC_N(\rho ), \\&C_N(\rho )=(1-\rho ^N)^{-1}\left( 1+N_t\varDelta t+N\varDelta t(1-\rho ^N)^{-1}\right) . \end{aligned} \end{aligned}$$

(112)

Proof

For all $n=0, \ldots , N_t$,

$$\begin{aligned} \left| \lambda _{s,n}^{\varDelta t, \varepsilon }-\lambda _{s,n}^{\varDelta t}\right|&\le S_1+S_2, \nonumber \\ S_1&=\left| \frac{\rho ^{-n}}{N}\sum \limits _{k=0}^{N-1}\mathrm {e}^{-i\frac{2\pi kn}{N}}\left( \varvec{\Lambda }^{s,\varepsilon }\left( \omega _k\right) -\varvec{\Lambda }^{s}\left( \omega _k\right) \right) \right| , \nonumber \\ S_2&=\left| \frac{\rho ^{-n}}{N}\sum \limits _{k=0}^{N-1}\mathrm {e}^{-i\frac{2\pi kn}{N}}\varvec{\Lambda }^{s}\left( \omega _k\right) -\lambda _{s,n}^{\varDelta t}\right| . \end{aligned}$$

(113)

An upper bound for $S_1$ follows from the triangle inequality and the assumption of the proposition: $S_1\le \rho ^{-n}\varepsilon \le \rho ^{-N_t}\varepsilon$ (because $\rho <1$).

As for $S_2$, it suffices to replace $\varvec{\Lambda }^{s}\left( \omega _k\right)$ in the above sum by $\sum \nolimits _{\ell =0}^{\infty }\lambda _{s,\ell }^{\varDelta t}\rho ^{\ell }\mathrm {e}^{i\frac{2\pi \ell k}{N}}$, cf. (92), and use the aliasing argument. In particular,

$$\begin{aligned} \frac{\rho ^{-n}}{N}\sum \limits _{k=0}^{N-1}\mathrm {e}^{-i\frac{2\pi kn}{N}}\varvec{\Lambda }^{s}\left( \omega _k\right) =\frac{\rho ^{-n}}{ N}\sum \limits _{k=0}^{N-1}\sum \limits _{\ell =0}^{\infty }\lambda _{s,\ell }^{\varDelta t} \, \rho ^{\ell } \, \mathrm {e}^{i\frac{2\pi k(\ell -n)}{N}}. \end{aligned}$$

Since $N^{-1}\sum \nolimits _{k=0}^{N-1}\mathrm {e}^{i\frac{2\pi k(\ell -n)}{N}}=1$ when $\ell -n$ is a multiple of N and vanishes otherwise, and $n\le N_t\le N-1$, the above can be rewritten as follows:

$$\begin{aligned}&\frac{\rho ^{-n}}{N}\sum \limits _{k=0}^{N-1}\mathrm {e}^{-i\frac{2\pi kn}{N}}\varvec{\Lambda }^{s}\left( \omega _k\right) =\lambda _{s, n}^{\varDelta t}+\rho ^{-n}\sum \limits _{k=1}^{\infty }\lambda _{s,kN+n}^{\varDelta t}\rho ^{kN+n},\quad \text { and }\\&S_2\le \rho ^{-n}\sum \limits _{k=1}^{\infty }\left| \lambda _{s,n+kN}^{\varDelta t}\right| \rho ^{n+kN}\overset{(95)}{\le }C\sum \limits _{k=1}^{\infty }\max (1, (n+kN)\varDelta t)\rho ^{kN}. \end{aligned}$$

The above sum is then bounded:

$$\begin{aligned} S_2&\le \sum \limits _{k=1}^{\infty }\rho ^{kN} +\sum \limits _{k=1}^{\infty }\rho ^{kN}(n+kN)\varDelta t\\&\le \rho ^N(1-\rho ^N)^{-1}(1+n\varDelta t)+N\varDelta t\,\rho ^N (1-\rho ^N)^{-2}. \end{aligned}$$

The result follows by bounding in the above $n\varDelta t$ by $N_t\varDelta t$ and combining bounds for $S_1$ and $S_2$ into (113). $\square$

The bound of Lemma C.1 allows us to quantify the choice of $\rho$, N in (93).

Proof of Proposition 4.3

The desired bound follows by applying the result of Lemma C.1. In particular, $C_N(\rho )$ can be estimated by providing an adequate estimate on $1-\rho ^N=1-\varepsilon ^{\frac{N}{N+N_t-1}}$. Because $N_t\ge 1$, the function $N\mapsto 1-\varepsilon ^{\frac{N}{N+N_t-1}}=1-\varepsilon \varepsilon ^{\frac{1-N_t}{N+N_t-1}}$ grows in N. Since, additionally, $N\ge N_t+1$, we have

$$\begin{aligned} 1-\rho ^N\ge 1-\varepsilon ^{\frac{N_t+1}{2N_t}}>1-\sqrt{\varepsilon }>1-\sqrt{\frac{1}{2}},\text { for all } 0<\varepsilon <\frac{1}{2}. \end{aligned}$$

Plugging in this bound into (112) yields $C_N(\rho )\le C(1+(N+N_t)\varDelta t)$ and

$$\begin{aligned} \max \limits _{n=0, \ldots , N_t}|\lambda _{s,n}^{\varDelta t, \varepsilon }-\lambda _{s,n}^{\varDelta t}|&<\varepsilon ^{\frac{N-1}{N+N_t-1}}+C\varepsilon ^{\frac{N}{N+N_t-1}}(1+(N+N_t)\varDelta t), \end{aligned}$$

from which the desired bound is obtained immediately. $\square$

Proof of Lemma 4.2

Let us show (a), which basically follows from Section 5.2.1 in [3]. The frequencies $\omega _k$ defined in (74), namely,

$$\begin{aligned} \omega _k=i\frac{\delta \left( \rho \mathrm {e}^{i\frac{2\pi k}{N}}\right) }{\varDelta t}=\frac{2i}{\varDelta t}\frac{1-\rho \mathrm {e}^{i\varphi _k}}{1+\rho \mathrm {e}^{i\varphi _k}}, \quad \varphi _k = \mathrm {e}^{i\frac{2\pi k}{N}}, \end{aligned}$$

lie on the circle centered at $c_{\rho , \varDelta t}$ of radius $R_{\rho , \varDelta t}$ (this follows from the fact that $z\mapsto \frac{1-z}{1+z}$ is a homography), with

$$\begin{aligned} c_{\rho ,\varDelta t}=\frac{2i}{\varDelta t}\frac{1+\rho ^2}{1-\rho ^2}, \quad R_{\rho ,\varDelta t}=\frac{2}{\varDelta t}\frac{2\rho }{1-\rho ^2}, \end{aligned}$$

(114)

i.e. $\omega _k=c_{\rho ,\varDelta t}+R_{\rho ,\varDelta t}\mathrm {e}^{i\psi _k}$, for some $\psi _k\in [0, 2\pi )$. Hence

$$\begin{aligned} {\text {Im}}\omega _k\ge \inf \limits _{\tiny 0\le \varphi <2\pi } {\text {Im}}\left( i\frac{\delta (\rho \mathrm {e}^{i\varphi })}{\varDelta t}\right) =\frac{2}{\varDelta t}\frac{1+\rho ^2-2\rho }{1-\rho ^2}=\frac{2}{\varDelta t}\frac{1-\rho }{1+\rho }, \end{aligned}$$

(115)

and, as $\rho <1$, ${\text {Im}}\omega _k>\frac{1-\rho }{\varDelta t}$. For $\rho$ defined in (97),

$$\begin{aligned} 1-\rho =1-\varepsilon ^{\frac{1}{2N_t}}=1-\exp \left( -\frac{\log \varepsilon ^{-1}}{2N_t}\right)>c_0\min \left( 1, \frac{\log \varepsilon ^{-1}}{N_t}\right) , \quad c_0>0, \end{aligned}$$

(116)

where the last bound follows from (111). Therefore, as $\varDelta t<1$, and $\varepsilon <\frac{1}{2}$,

$$\begin{aligned} {\text {Im}}\omega _k>c\min \left( 1, \frac{1}{N_t\varDelta t}\right) ,\quad c>0. \end{aligned}$$

To show (b), we use the same property (114), which results in

$$\begin{aligned} |\omega _k|\le \frac{2}{\varDelta t}\left( \frac{1+\rho ^2}{1-\rho ^2}+\frac{2\rho }{1-\rho ^2}\right) \le \frac{2(1+\rho )}{\varDelta t(1-\rho )}<\frac{4}{\varDelta t(1-\rho )}. \end{aligned}$$

Using (116), and then $\varepsilon <\frac{1}{2}$, we deduce the following inequality, for some $C, C'>0$,

$$\begin{aligned} |\omega _k|<\frac{C}{\varDelta t}\max \left( 1, N_t \left( \log \varepsilon ^{-1}\right) ^{-1}\right) \le \frac{C'}{\varDelta t} \max \left( 1, N_t\right) \le \frac{C'}{N_t\varDelta t}N_t^2. \end{aligned}$$

$\square$