Numerical approximation of a linear elasticity model

Integro-differential equations arise while modeling many problems in physics, biology and other areas which involve non-local diffusion/dispersal mechanisms [2, 10, 11]. A nonlocal elastic wave model can be written as [11]

$$\frac{\partial^{2}}{\partial t^{2}}u(x,t)=\rho\mathcal{L}u(x,t)+g(x,t)\quad(t,x)\in\Omega\times(0,T),$$

(1)

where

$$\mathcal{L}\phi(x)=\int_{\Omega}\mathcal{J}^{\infty}(x-y)\phi(y)dy+\alpha(x)\phi(x),$$

$(0,T)$, $T>0$ is the time interval under consideration, $\alpha(x)$ is a time-independent coefficient and $g(x,t)$ is an inhomogeneity. Here we consider $u(x,0)=u_{0}$, and $\frac{\partial}{\partial t}u(x,0)=v_{0}$, $\Omega\subset\mathbb{R}$, where $u(x,t)$ denotes the displacement field and $\partial_{t}u(x,t)$ represents velocity. Here the $L^{1}(\mathbb{R})$ kernel $\mathcal{J}^{\infty}(x)$ satisfies $\mathcal{J}^{\infty}(x)=\mathcal{J}^{\infty}(-x)$ (symmetry).

We assume that $\alpha(x)$ is a constant given by $\alpha=-\int_{\Omega}\mathcal{J}^{\infty}(x)dx,$ and for simplicity we consider a normalized kernel function ($\int\mathcal{J}^{\infty}(x)dx=1$). Then (1) can be written as

$$\frac{\partial^{2}}{\partial t^{2}}u(x,t)=\rho\mathcal{L}u(x,t)+g(x,t)$$

(2)

where

$$\mathcal{L}\phi(x)=\int_{\Omega}\mathcal{J}^{\infty}(x-y)(\phi(y)-\phi(x))dy.$$

Here we restrict ourself by $\mathcal{J}^{\infty}\geq 0$ (nonnegative) as well.

Now let $y-x=z$, then $dy=dz$. Considering $\Omega=\mathbb{R}$ [3] in (2) we have

$$I_{1}=\varepsilon\int_{\mathbb{R}}\mathcal{J}^{\infty}(x-y)\left(u(y,t)-u(x,t)\right)dy=\varepsilon\int_{\mathbb{R}}\mathcal{J}^{\infty}(z)\left(u(x+z,t)-u(x,t)\right)dz.$$

Expanding $u(x+z)$ in terms of Taylor series around $z=0$ we have

$$I_{1}=\varepsilon\int_{\mathbb{R}}\mathcal{J}^{\infty}(z)\left(u(x,t)+zu_{x}(x,t)+\frac{z^{2}}{2!}u_{xx}(x,t)+\cdots\cdots-u(x,t)\right)dz.$$

Thus we have

$$u_{tt}=\varepsilon u_{x}\int_{\mathbb{R}}z\mathcal{J}^{\infty}(z)dz+\frac{\varepsilon u_{xx}}{2!}\int_{\mathbb{R}}z^{2}\mathcal{J}^{\infty}(z)dz+\cdots\cdots+g(x,t).$$

(3)

As $J(z)$ is symmetric, integrals with odd powers of $z$ are zero. So (3) takes the form $u_{tt}=C_{2}u_{xx}+C_{4}u_{xxxx}+...+g(x,t).$ where $C_{2m}=\frac{\varepsilon}{2m!}\int_{\mathbb{R}}z^{2m}\mathcal{J}^{\infty}(z)dz$; for each $m=1,2,3,\cdots.$ Ignoring the higher order terms we have

$$u_{tt}=C_{2}u_{xx}+g(x,t).$$

(4)

The above equation is the familiar non-homogeneous wave equation, and is a first order approximation of the non-local equation (2). A comparison between the operators $u_{xx}$ and $\mathcal{L}u$ is well presented in [3, 8]. Now following [3] it is well understood that $\|\mathcal{L}\|$ is bounded and $\mathcal{L}$ is negative semi-definite, if the kernel function $\mathcal{J}^{\infty}$ is considered non-negative, symmetric and normalized. These properties are important for the analysis we performed in this study.

Emmrich et. al. [11] consider the non-local elastic model (2). They analyze mathematical and numerical solutions of the model. They compare the non-local model with that of a local continuum model. They claim two advantages of using this non-local model. They propose and implement a quadrature based numerical scheme considering infinite spatial domain.

It is well understood from [4, 9] that an integro-differential equation of type (2) defined in the infinite domain can also be defined in a truncated finite domain $[A,B]$ where $A$ and $B$ depend on the decay of the kernel function $J^{\infty}(x)$. It is to note that a closed form formula to find suitable $A$ and $B$ is well presented in [9]. Thus the analysis and the approximation in a periodic $\Omega=[0,1]$ spatial domain can also be applied in any bounded interval $[A,B]$ (and vice-versa). Considering the kernel of the convolution integral as $g_{\delta}(x)=\sqrt{\frac{1}{2\pi\delta}}\exp\left(-\frac{y^{2}}{2\delta^{2}}\right),$ the author in [9] formulate $A=+\sqrt{-2\delta^{2}\log(\delta\varepsilon\sqrt{2\pi})},$ and $B=-A,$ where $\varepsilon>0$ considered so that $g_{\delta}(x)\geq\varepsilon$.

One can consider the model in a spatial periodic domain [3] as well. If we consider a spatially one-periodic initial function $u(x,0)$, then for all $x\in\mathbb{R}$ and $t\in\mathbb{R}_{+}$ $u(x,t)=u(x+L,t).$ Then (2) can be written as

$\displaystyle u_{tt}$

$\displaystyle=$

$\displaystyle(\mathcal{L}u(\cdot,t))(x)+g(x,t)=\int_{\Omega}\mathcal{J}(x-y)\left(u(y,t)-u(x,t)\right)dy+g(x,t)$

(5)

where $\mathcal{J}(x)=\sum_{r=-\infty}^{\infty}\mathcal{J}^{\infty}(x-rL)$ for all $x\in[0,L]$, $L>0$. We are interested to consider both periodic and non-periodic domain for spatial approximations of the model.

Convolution models have been extensively investigated both theoretically and numerically [5, 11], but a lot more yet to be done to get a highly accurate solver. Since an analytic solution can not be obtained always, an efficient and highly accurate numerical scheme has to be developed. Spatial approximation of this type of convolution models are interesting as well as challenging. The unknown function under integral sign and the nonlinearity involved in the model make the approximation more challenging. There are various ways to handle such problems. In most cases scientists use some lower order schemes (with midpoint quadrature rules for integration [3, 6]) to serve their purpose. Thus there is much room for improvement and we find an interest of presenting and analyzing a higher order technique for space integration. To the best of our knowledge, Legendre spectral Galerkin method has not been performed for this model problem yet and we find an interest to investigate such a scheme.

In this article, we start with numerical accuracy analysis of a spectral galerkin approximation in section 2 followed by numerical implementations of the scheme in section 3. We present several computer generated solutions with relevent discussion in section 4. We implement the schemes in MATLAB.

In this section, we look for spectral Galerkin scheme for space approximation and investigate convergence of such a scheme. Here we consider Legendre spectral polynomials. It is our goal to find $u=u(x,t)$, $x\in\Omega$ such that

$$u_{tt}(x,t)=\mathcal{L}u(x,t)+g(x,t),\forall x\in\Omega,$$

and we seek for weak form to find $u\in L_{2}(\Omega)$ such that

$$(u_{tt},v)=(\mathcal{L}u(x,t),v)+(g(x,t),v),\quad\forall\quad v\in L_{2}(\Omega).$$

(6)

We denote $\mathbb{N}$ be the set of all nonnegative integers. For any $N\in\mathbb{N}$, $\mathbb{P}_{N}$ denotes the set of all algebraic polynomials of degree at most $N$ in $\Omega$. We define Legendre spectral Galerkin approximation of (2) as: Find

$$u^{N}\in\mathbb{P}_{N},\quad\text{such that}$$

$$(u^{N}_{tt},v)=(\mathcal{L}u^{N}(x,t),v)+(g(x,t),v),\quad\forall\quad v\in\mathbb{P}_{N}(\Omega).$$

(7)

We define projection operator $\mathbb{P}_{N}$ such that

$$(u,v)=(P_{N}u,v),\quad\forall\quad v\in\mathbb{P}_{N}.$$

Then [12, 7]

$$\|u-P_{N}u\|\leq CN^{-m}|u|_{H^{m}(\Omega)}\leq CN^{-m}\|u\|_{H^{m}(\Omega)}$$

where

$$\|u\|_{H^{m}(\Omega)}=\left(\sum_{k=0}^{m}\|u^{k}\|\right)^{1/2},$$

and

$$|u|_{H^{m}(\Omega)}=\left(\sum_{k=\min{m,N+1}}^{m}\|u^{k}\|\right)^{1/2},$$

and

$$\|u\|_{m,\infty}=\max_{0\leq k\leq m}\left\{\|u^{k}\|_{\infty}\right\},$$

and $C$ denotes a nonnegative constant which is independent of $N$.

Let $e^{N}=u-u^{N}$ be the error in spectral Galerkin approximation to the Legendre spectral Galerkin solution $u^{N}$ of (6). Now from (6) and (7) we have

$$((u-u^{N})_{tt},v)=(\mathcal{L}(u-u^{N}),v).$$

(8)

We need the following results to prove Theorem 1.

Here we implement spectral Galerkin schemes considering Legendre spectral polynomials. The Legendre polynomials satisfy three term recurrence relation:

$$L_{0}(x)=1,$$

$$L_{1}(x)=x,$$

and

$$(n+1)L_{n+1}(x)=(2n+1)xL_{n(x)}-nL_{n-1}(x),\quad n\geq 1,$$

with

$$\int_{-1}^{1}L_{k}(x)L_{j}(x)dx=\frac{2}{2k+1}\delta_{kj},$$

$$\delta_{kj}=\left\{\begin{array}[]{cc}1&\text{if}\quad k=j\\ 0&\text{if}\quad k\neq j.\end{array}\right.$$

Here we consider $\Omega=[-1,1]$ for numerical implementations and one may extend the scheme to any interval $[-L,L]$, $L>0$. We define

$$u^{N}(x,t)=\sum_{k=0}^{N}a_{k}(t)L_{k}(x).$$

Replacing $v\in\mathbb{P}_{N}$ by $L_{k}\in\mathbb{P}_{N}$ in (6) we have

$$(u^{N}_{tt}(x,t),L_{k}(x))=(\mathcal{L}u^{N}(x,t),L_{k}(x))+(g(x,t),L_{k}(x)),\quad\forall\quad L_{k}(x)\in\mathbb{P}_{N}(\Omega).$$

(10)

Now

	$\displaystyle(u^{N}_{tt}(x,t),L_{k}(x))$	$\displaystyle=$	$\displaystyle\sum_{k=0}^{N}\frac{d^{2}}{dt^{2}}a_{j}(t)(L_{j}(x),L_{k}(x))$
		$\displaystyle=$	$\displaystyle\frac{2}{2k+1}\frac{d^{2}}{dt^{2}}a_{k}(t),$

and

$$(g(x,t),L_{k}(x))=\sum_{l=0}^{N}w_{l}g(x_{l})L_{k}(x_{l}).$$

where $w_{l}$ are Gauss weights and $x_{l}$ are Gauss quadrature points. We define

	$\displaystyle\mathcal{L}u(x,t)$	$\displaystyle=$	$\displaystyle\mathcal{L}_{1}u(x,t)+\mathcal{L}_{2}u(x,t)$
		$\displaystyle=$	$\displaystyle\int\mathcal{J}^{\infty}(x-y)u(y,t)dy-u(x,t)\int\mathcal{J}^{\infty}(x-y)dy.$

Now

$$(\mathcal{L}u^{N}(x,t),L_{k}(x))=(\mathcal{L}_{1}u^{N}(x,t),L_{k}(x))+(\mathcal{L}_{2}u^{N}(x,t),L_{k}(x))$$

where

$$(\mathcal{L}_{1}u^{N}(x,t),L_{k}(x))=\sum_{j}a_{j}(t)\int\int\mathcal{J}^{\infty}(x-y)L_{k}(x)L_{j}(y)dydx$$

and

$$(\mathcal{L}_{2}u^{N}(x,t),L_{k}(x))=\sum_{j}a_{j}(t)\int L_{k}(x)L_{j}(x)\int\mathcal{J}^{\infty}(x-y)dydx.$$

Since we consider $\mathcal{J}^{\infty}$ such that

$$\int_{\Omega}\mathcal{J}^{\infty}(x)=1,\quad\forall\quad x\in\Omega,$$

we have

$$(\mathcal{L}_{2}u^{N}(x,t),L_{k}(x))=\sum_{j}a_{j}(t)\int L_{k}(x)L_{j}(x)dx=\frac{2}{2k+1}a_{k}(t).$$

Adding all these integrals (10) can be written as

	$\displaystyle\frac{2}{2k+1}\frac{d^{2}}{dt^{2}}a_{k}(t)$	$\displaystyle=$	$\displaystyle\sum_{j=0}^{N}a_{j}(t)\int_{\Omega}\int_{\Omega}\mathcal{J}^{\infty}(x-y)L_{k}(x)L_{j}(y)dydx$		(11)
			$\displaystyle\quad-\frac{2}{2k+1}a_{k}(t)+\int g(x,t)L_{k}(x)dx,$

for all $k=0,1,\cdots,N$. Thus (11) can be written as

$$M\frac{d^{2}}{dt^{2}}a(t)=Aa(t)+b(t),$$

with $a(0)=M^{-1}u_{0}$, and $a^{\prime}(0)=M^{-1}v_{0}$, $A=A_{1}+A_{2}$ and elements of $A_{1}$ are defined from the continuous operator $\mathcal{L}_{1}$ as well as the elements of $A_{2}$ are defined from the continuous operator of $\mathcal{L}_{2}$. Since $M$ is a nonsingular diagonal matrix, the above second order ordinary differential equation can be written as

$$\frac{d^{2}}{dt^{2}}a(t)=M^{-1}Aa(t)+M^{-1}b(t).$$

(12)

We approximate the integrals by Gaussian quadrature formula

	$\displaystyle\int_{\Omega}\int_{\Omega}\mathcal{J}^{\infty}(x-y)L_{k}(x)L_{j}(y)dydx$	$\displaystyle\approx\sum_{l=0}^{M}w_{l}L_{k}(x_{l})\int_{\Omega}\mathcal{J}^{\infty}(x_{l}-y)L_{j}(y)dy$
		$\displaystyle\approx\sum_{l=0}^{M}\sum_{m=0}^{M}w_{l}L_{k}(x_{l})w_{m}\mathcal{J}^{\infty}(x_{l}-x_{m})L_{j}(x_{m}),$

and

$\displaystyle\int_{\Omega}g(x,t)L_{k}(x)dx$

$\displaystyle\approx\sum_{l=0}^{M}w_{l}g(x_{l},t)L_{j}(x_{l}),$

where $x_{l}$ are quadrature points, and $w_{l}$ are quadrature weights. (11) is a system of second order time dependent ordinary differential equations which can be solved using a simple finite difference scheme to approximate the coefficients $a_{k}(t)$.

Depending on the nature of the kernel function, it is also wise to divide the integral domain into pieces. Then one can use small number of quadrature points to compute the integrals, since over each subinterval/subdomain $J(x)$ will be flatter (if $J(x)$ is a smooth function over $\Omega$.) We discuss the quadrature approximations in Section A.

In this section, we exhibit several numerical results to illustrate the scheme. Here we use a simple scheme for time integration. We use MATLAB to serve our purpose. We approximate the time dependent system of equations (12) using a central difference scheme ¹¹1Mathematica and MATLAB builtin functions can also be used to solve the time dependent system of equations (12).. For simplicity we consider $\Omega=[-1,1]$ and the scheme can be extended to $[-L,L]$, for any $L>0$. We approximate (12) by

$$\frac{a^{j+1}-2a^{j}+a^{j-1}}{dt^{2}}=M^{-1}Aa^{j+1}+M^{-1}b(t_{j})$$

(13)

where $a^{j}=a(t_{j}).$ This scheme is accurate [1] of order $\mathcal{O}(\Delta t^{2}+N^{-m})$.

We know that $\rho(A)\leq\|A\|$ for any matrix $A$ where $\rho(A)$ stands for condition number of $A$. Also $\|\frac{1}{I-\Delta t^{2}M^{-1}A}\|\leq\frac{1}{1-\|dt^{2}M^{-1}A\|}.$ Now

	$\displaystyle|A_{jk}^{1}|=|\int\int\mathcal{J}^{\infty}(x-y)L_{k}(x)L_{j}(y)dydx|$	$\displaystyle\leq\max_{x\in\Omega}\mathcal{J}^{\infty}(x)|\int\int L_{k}(x)L_{j}(y)dydx|$
		$\displaystyle\leq|\int L_{k}(x)dx||\int L_{j}(y)dy|\leq 4K,$

since

$$\int_{\Omega}L_{j}(x)dx=\left\{\begin{array}[]{cc}2&j=0\\ 0&j\neq 0,\end{array}\right.,$$

$\mathcal{J}^{\infty}$ is a normalized function and $K=\max_{x\in\Omega}|\mathcal{J}^{\infty}(x)|$. We have

$$\|A_{1}\|=\sqrt{\sum_{i,j=1}^{N+1}|A_{ij}^{1}|^{2}}\leq\sqrt{\sum_{i,j=1}^{N+1}16K}=\sqrt{16(N+1)^{2}K^{2}}=4(N+1)K,$$

and

$$\|A_{2}\|=\sqrt{\sum_{j,j=1}^{N+1}|A^{2}_{jj}|^{2}}\leq 2\sqrt{\sum_{j=1}^{N+1}\left|\frac{1}{2j-1}\right|^{2}}<4.$$

Thus

$$\|A\|\leq\|A_{1}\|+\|A_{2}\|\leq 4\left((N+1)K+1\right),$$

and $1\leq\|M\|\leq 4$, shows that $\|M^{-1}A\|\leq 4\left((N+1)K+1\right)$ which is a bounded quantity for any finite $N$. So $\Delta t^{2}\|M^{-1}A\|\longrightarrow 0$ as $\Delta t\longrightarrow 0$. Hence

$$\rho\left(\frac{1}{I-\Delta t^{2}M^{-1}A}\right)\leq\|\frac{1}{I-\Delta t^{2}M^{-1}A}\|\leq\frac{1}{1-\|\Delta t^{2}M^{-1}A\|}\leq 1,$$

and thus by the definition of condition number of a matrix

$$\rho\left(\frac{1}{I-\Delta t^{2}M^{-1}A}\right)=1,\ \text{as}\ \Delta t\rightarrow 0\ \text{and}\ N<\infty.$$

Thus the scheme (13) is unconditionally stable.

In the following Figures 1-2, we present the numerical solutions for various choices of initial functions, $\rho$, and $g(x,t)$ to demonstrate the scheme considering $\mathcal{J}^{\infty}(x)=\sqrt{\frac{400}{\pi}}e^{-400x^{2}}$. Here we notice that the wave form is similar to the solutions of local linear wave model and the numerical solutions behaves well for any finite time preserving wave phenomena.

Here in this study we have shown theoretically that the error in spectral approximation decays exponentially. More precisely, it is shown that if the kernel function and the solution underlying the model problem is smooth, then the errors obtained by the Legendre spectral scheme decays exponentially, which one desires from spectral method.

There is a drawback of this approximation. Because of the presence of the convolution operator, the discrete operator is a dense matrix and as a result numerical computations become slower, needs huge storage for computations in multidimensional domain.