On the numerical solution of the elastodynamic problem by a boundary integral equation method

TL;DR

This paper introduces a numerical boundary integral method combining Laguerre transformation and quadrature for solving elastodynamic problems in unbounded 2D regions, demonstrated through numerical results.

Contribution

It presents a novel approach that reduces time-dependent elastodynamic problems to stationary boundary value problems using Laguerre transformation, solved via boundary integral equations.

Findings

Effective numerical solution for elastodynamic problems

Reduction of time-dependent problems to stationary boundary problems

Numerical results validate the method's accuracy

Abstract

A numerical method for the Dirichlet initial boundary value problem for the elastic equation in the exterior and unbounded region of a smooth closed simply connected 2-dimensional domain, is proposed and investigated. This method is based on a combination of a Laguerre transformation with respect to the time variable and a boundary integral equation approach in the spatial variables. Using the Laguerre transformation in time reduces the time-depended problem to a sequence of stationary boundary value problems, which are solved by a boundary layer approach resulting to a sequence of boundary integral equations of the first kind. The numerical discretization and solution are obtained by a trigonometrical quadrature method. Numerical results are included.

Tables

Table 1. Table 1: The first components of the computed and the exact solutions of ( 2.5a )–( 2.5c ), for the specific boundary function ( 5.1 ), at the measurement point y = ( 1.5 , 1 ) . 𝑦 1.5 1 y=(1.5,1).

$M$	${\left({\tilde{u}}_{0,M}\right)}_1\left(y\right)$	${\left({\tilde{u}}_{1,M}\right)}_1\left(y\right)$	${\left({\tilde{u}}_{2,M}\right)}_1\left(y\right)$
8	$0.293581559232289$	$-0.084483725080856$	$-0.146079079028772$
16	$0.284988364785089$	$-0.092525310787524$	$-0.155666923858005$
32	$0.285503199323624$	$-0.092138738384510$	$-0.155404881866276$
64	$0.285503741272164$	$-0.092138337605882$	$-0.155404627504594$
	${\left(u_0^{ex}\right)}_1\left(y\right)$	${\left(u_1^{ex}\right)}_1\left(y\right)$	${\left(u_2^{ex}\right)}_1\left(y\right)$
	$0.285503741272020$	$-0.092138337605708$	$-0.155404627504139$

Table 2. Table 2: The second components of the computed and the exact solutions of ( 2.5a )–( 2.5c ), for the specific boundary function ( 5.1 ), at the measurement point y = ( 1.5 , 1 ) . 𝑦 1.5 1 y=(1.5,1).

$M$	${\left({\tilde{u}}_{0,M}\right)}_2\left(y\right)$	${\left({\tilde{u}}_{1,M}\right)}_2\left(y\right)$	${\left({\tilde{u}}_{2,M}\right)}_2\left(y\right)$
8	0.081036497084071	0.028667287783745	$-0.011013685642670$
16	0.071649152048006	0.017647071803846	$-0.021814816353654$
32	0.071756738147012	0.017837149908881	$-0.021648258690479$
64	0.071756880072043	0.017837482038337	$-0.021647898034988$
	${\left(u_0^{ex}\right)}_2\left(y\right)$	${\left(u_1^{ex}\right)}_2\left(y\right)$	${\left(u_2^{ex}\right)}_2\left(y\right)$
	0.071756880072350	0.017837482039221	$-0.021647898033812$

Table 3. Table 3: Numerical values of the components of the computed ( u N , M ) 1 subscript subscript 𝑢 𝑁 𝑀 1 (u_{N,M})_{1} and the exact solution ( u e ​ x ​ ( N ) ) 1 subscript superscript 𝑢 𝑒 𝑥 𝑁 1 (u^{ex}(N))_{1} (rows in grey) of the problem ( 1.1 )–( 1.4 ) for the boundary function ( 5.2 ). Here κ = 1 / 2 , 𝜅 1 2 \kappa=1/2, and y = ( 1 , 1 ) . 𝑦 1 1 y=(1,1).

$t$	$M$	$N=15$	$N=20$	$N=25$
1	32	0.556497018896238	0.589619531491177	0.661666532791282
	64	0.556495142719243	0.589617717103517	0.661663974784660
	${\left(u^{ex}\right)}_1$	0.556495142721411	0.589617717107183	0.661663974772766
2	32	0.512733429371186	0.536980500925658	0.478568254585875
	64	0.512733615180447	0.536980917568962	0.478569337106983
	${\left(u^{ex}\right)}_1$	0.512733615179871	0.536980917566782	0.478569337119777
3	32	0.240000095287178	0.141039343209334	0.117399715612992
	64	0.240001133548879	0.141039958261864	0.117400385116931
	${\left(u^{ex}\right)}_1$	0.240001133547631	0.141039958259843	0.117400385110933

Table 4. Table 4: Numerical values of the components of the computed ( u N , M ) 2 subscript subscript 𝑢 𝑁 𝑀 2 (u_{N,M})_{2} and the exact solution ( u e ​ x ​ ( N ) ) 2 subscript superscript 𝑢 𝑒 𝑥 𝑁 2 (u^{ex}(N))_{2} (rows in grey) of the problem ( 1.1 )–( 1.4 ) for the boundary function ( 5.2 ). Here κ = 1 / 2 , 𝜅 1 2 \kappa=1/2, and y = ( 0.5 , − 1.5 ) . 𝑦 0.5 1.5 y=(0.5,-1.5).

$t$	$M$	$N=15$	$N=20$	$N=25$
1	32	$-0.027394214935554$	$-0.029780601479552$	$-0.037360920446497$
	64	$-0.027394199243065$	$-0.029780534307372$	$-0.037360395287884$
	${\left(u^{ex}\right)}_2$	$-0.027394199243094$	$-0.029780534308750$	$-0.037360395289490$
2	32	$-0.006641666911909$	$-0.007677796534028$	$-0.001408341947825$
	64	$-0.006641661641063$	$-0.007677775150468$	$-0.001408727582784$
	${\left(u^{ex}\right)}_2$	$-0.006641661641210$	$-0.007677775149996$	$-0.001408727581196$
3	32	0.008544988585762	0.014762278037365	$0.016882729931766$
	64	0.008544972022448	0.014762136068806	0.016882546539029
	${\left(u^{ex}\right)}_2$	0.008544972022701	0.014762136070815	$0.016882546538771$

Table 5. Table 5: The first component of the numerical solution of the problem ( 1.1 )–( 1.4 ) for the boundary function ( 5.4 ), for κ = 1 / 2 , 𝜅 1 2 \kappa=1/2, at the measurement position y = ( 0.5 , − 1.5 ) . 𝑦 0.5 1.5 y=(0.5,-1.5).

$t$	$M$	$N=10$	$N=15$	$N=20$
1	16	0.062572635051987	0.047833811805535	0.046009612087451
	32	0.062555821083946	0.047821505207762	0.045995757574454
	64	0.062555817944294	0.047821498423333	0.045995752862933
2	16	0.206527137283687	0.230300799019877	0.230584463332204
	32	0.206507930918013	0.230271954401413	0.230552865835900
	64	0.206507936378995	0.230271967614367	0.230552881497197
3	16	0.281765785091472	0.290022752671760	0.293418322035422
	32	0.281753694023185	0.290010686252088	0.293412971755498
	64	0.281753701138142	0.290010693385048	0.293412971498783

Table 6. Table 6: The second component of the numerical solution of the problem ( 1.1 )–( 1.4 ) for the boundary function ( 5.4 ), for κ = 1 / 2 , 𝜅 1 2 \kappa=1/2, at position y = ( 0.5 , − 1.5 ) . 𝑦 0.5 1.5 y=(0.5,-1.5).

$t$	$M$	$N=10$	$N=15$	$N=20$
1	16	0.104385591001347	0.089274523790871	0.089029847496310
	32	0.104414411445436	0.089311568889684	0.089061508274427
	64	0.104414399045771	0.089311556699605	0.089061496285464
2	16	0.242963484605096	0.268615650361266	0.270356959200451
	32	0.242960849558260	0.268601030910481	0.270336784091780
	64	0.242963476811622	0.268601022577740	0.270336775642960
3	16	0.294596074694631	0.301687858830015	0.299978778147128
	32	0.294584218004291	0.301675713615639	0.299984785021505
	64	0.294584215781299	0.301675711488577	0.299984782636057

Equations

$$\frac{{\partial }^{2}u}{\partial t^2}={\Delta }^{\ast }u\text{in}D\times \left(0,\infty \right),$$

$$\frac{\partial u}{\partial t}(\cdot ,0)=u(\cdot ,0)=0\text{in}D,$$

$$u=f,\text{on}\Gamma \times (0,\infty ),$$

$$f(x,0)=\frac{\partial f}{\partial t}(x,0)=0,x\in \Gamma .$$

$$c_s=\sqrt{\frac{\mu }{\rho }},c_p=\sqrt{\frac{\lambda +2\mu }{\rho }},$$

$$u(x,t)\to 0,|x|\to \infty ,$$

$$L_n(t)=\sum _{k=0}^n\left(\genfrac{}{}{0ex}{}{n}{k}\right)\frac{(-t{)}^k}{k!},$$

$$(n+1)L_{n+1}(t)=(2n+1-t)L_n(t)-n{L}_{n-1}(t).$$

$$L_{n+1}^{\prime }(t)=L_n^{\prime }(t)-L_n(t),n=0,1,\dots ,$$

$$v(t)=\kappa \sum _{n=0}^{\infty }{v}_n{L}_n(\kappa t)$$

$$v_n=\underset{0}{\overset{\infty }{\int }}{e}^{-\kappa t}{L}_n(\kappa t)v(t)dt,n=0,1,\dots .$$

$$v_n^{\prime }=\kappa \sum _{m=0}^n{v}_m$$

$$v_n^{\prime \prime }={\kappa }^2\sum _{m=0}^n(n-m+1)v_m.$$

$${\Delta }^{\ast }{u}_n-{\kappa }^2{u}_n$$

$$u_n$$

$$u_n(x)$$

$$u(x,t)=\kappa \sum _{n=0}^{\infty }{u}_n(x)L_n(\kappa t),$$

$${\Delta }^{\ast }{E}_n(x,y)-\sum _{m=0}^n{\beta }_{n-m}{E}_m(x,y)=\delta (x-y)I.$$

$$v_n(\gamma ,r)=\sum _{m=0}^{\left[\frac{n}{2}\right]}{a}_{n,2m}(\gamma )r^{2m},w_n(\gamma ,r)=\sum _{m=0}^{\left[\frac{n-1}{2}\right]}{a}_{n,2m+1}(\gamma )r^{2m+1}$$

$$a_{n,0}(\gamma )$$

$$a_{n,n}(\gamma )$$

$$a_{n,m}(\gamma )=\frac{1}{2\gamma m}\left\{4{\left[\frac{m+1}{2}\right]}^2{a}_{n,m+1}(\gamma )-{\gamma }^2\sum _{k=m-1}^{n-1}(n-k+1)a_{k,m-1}(\gamma )\right\},$$

$${\Phi }_n(\gamma ,r)=K_0(\gamma r)v_n(\gamma ,r)+K_1(\gamma r)w_n(\gamma ,r)$$

$${\int }_{r/a}^{\infty }\frac{e^{-\kappa t}{L}_n(\kappa t)}{\sqrt{t^2-(r/a{)}^2}}dt={\Phi }_n(\frac{\kappa }{a},r)$$

$${\int }_{r/a}^{\infty }\frac{e^{-\kappa t}(\kappa t{)}^2{L}_n(\kappa t)}{\sqrt{t^2-(r/a{)}^2}}dt=\sum _{k=-2}^2{\chi }_{k,n}{\Phi }_{n+k}(\frac{\kappa }{a},r),$$

$$G(x,y;t)=\frac{\theta (t-|x-y|/a)}{\sqrt{t^2-(|x-y|/a{)}^2}},$$

$$\frac{1}{a^2}\frac{{\partial }^{2}G(x,y;t)}{\partial t^2}-\Delta G(x,y;t)=\delta (x-y)\delta (t).$$

$$\Delta G_n(x,y)-\frac{{\kappa }^2}{a^2}\sum _{m=0}^n(n-m+1)G_m(x,y)=\delta (x-y),n=0,1,\dots$$

$$G_n(x,y)={\int }_0^{\infty }G(x,y;t)e^{-\kappa t}{L}_n(\kappa t)dt.$$

$$t^2{L}_n(t)=\sum _{k=-2}^2{\chi }_{k,n}{L}_{n+k}(t).$$

$$E_n(x,y)={\Phi }_{1,n}(|x-y|)I+{\Phi }_{2,n}(|x-y|)J(x-y)$$

$${\Phi }_{\ell ,n}(r)$$

Full text

Keywords: elastic equation; initial boundary value problem; Laguerre transformation; fundamental sequence; single and double layer potentials; boundary integral equations of the first kind; trigonometrical quadrature method.

1 Introduction

The problem of numerically solving time-dependent boundary value problems has a long history. The simplest approach consists of using the FDM, which has the obvious limitations of simple domains and interior problems. The most used general scheme, for this kind of problems, first reduces the dimensions of the problem by some semi-discretization approach and then solves the simpler problem by a numerical method. For example, using the Galerkin method with respect to the spatial variables, as a semi-discretization technique, results to a Cauchy problem for a system of ordinary differential equations. Then, single- or multi-steps methods can be applied.

On the other hand, we can apply semi-discretization with respect to the time variable and reduce the given problem to a set of stationary problems. This can be done by an integral transformation or by Rothe’s method. Then, we can apply a numerical method suitable for stationary problems, for example FEM, FDM, BEM and others.

The main aim and result of this paper is to present an integral equation method for the time-dependent elastic equation in an unbounded two-dimensional domain. The possible variants to use an integral equation approach for time-dependent boundary value problems are discussed in [7]. If the given differential equation has a fundamental solution, the problem can be reduced to a time-boundary integral equation by direct or indirect methods (see [1] for the elastic equation). Also by second semi-discretization approach the method of integral equation can be applied for receiving stationary problems.

Our goal is to extend the idea of the Fourier-Laguerre expansion of the solution and consider it as a semi-discretization approach with respect to time [2, 3, 4, 5, 6, 8]. Then, we need to solve a recurrent sequence of stationary problems for Navier equations. The advantage of this semi-discretization approach, compared to Laplace transform, is the simple representation of the numerical solution. In our case, it is a partial sum of the Fourier-Laguerre series. Since the solution domain is unbounded, the use of classical integral equation methods with a volume potential would be ineffective. Instead, we present the solution in terms of a boundary potentials using explicit fundamental solutions of the obtained sequence. This approach gives us the possibility to reduce the stationary differential problems to a sequence of boundary integral equations.

The outline of the paper is as follows: In Section 2, we describe the semi-discretization procedure in time via the Laguerre transformation. Then, the initial boundary value problem for the elastic equation is transformed to a sequence of boundary value problems for Navier equations. In Section 3, we show the reduction of the sequence of stationary problems to integral equations merely involving boundary potentials. To do so, we created fundamental solutions for the sequence of stationary equations. In Section 4, we outline and describe how the well-established numerical methods based on trigonometrical quadratures can be adjusted and applied to solve numerically the derived sequence of boundary integral equations. In the last Section 5, we demonstrate the feasibility of our approach through three numerical examples.

Before closing this section we formulate the problem to be studied. Let $D\subset{\rm I\!R}^{2}$ be a unbounded domain such that its complement is bounded and simply connected and assume that the boundary $\Gamma$ of $D$ is of class $C^{2}$. Consider the initial boundary value problem for the hyperbolic elastic equation

[TABLE]

with the Lamé operator defined by $\Delta^{*}:=c_{s}^{2}\Delta+(c_{p}^{2}-c_{s}^{2})\nabla\nabla\cdot$ , supplied with the homogeneous initial conditions

[TABLE]

and the boundary condition

[TABLE]

where $f$ is a given function, satisfying the compatibility condition

[TABLE]

Here the velocities $c_{s}$ and $c_{p}$ have the following form

[TABLE]

where $\rho$ is the density, and $\lambda$ and $\mu$ are the Lamé constants.

Since we have an unbounded solution domain, we specify that at infinity

[TABLE]

uniformly with respect to all directions $\frac{x}{{|x|}}$ and all $t\in[0,\infty)$. Other types of boundary conditions can potentially be handled as well.

2 Semi-discretization in the time

For the semi-discretization with respect to the time variable in the problem (1.1)–(1.4) we use the Laguerre transformation.

The Laguerre polynomials have the explicit representation

[TABLE]

for $n=0,1,\ldots$. These polynomials are orthogonal with respect to the standard $L^{2}$ inner product in $(0,\infty)$ appended with the weight $e^{-t}$. The following recurrence relation holds

[TABLE]

For our purpose though, we need to note the relation

[TABLE]

which is immediate from the explicit representation of $L_{n}$.

A function being square integrable with respect to the above weight $e^{-t}$ over the interval $(0,\infty)$ can then be expanded in a scaled Fourier-Laguerre series

[TABLE]

with Fourier-Laguerre coefficients

[TABLE]

These two relations are interpreted as the inverse and direct Laguerre transformations, respectively. Here, $v$ is the given function and the sequence $\{v_{n}\}$ is the image; $\kappa>0$ is a fixed scaling parameter. In the rest of this work, we use the above scaled convention as the Laguerre transform.

Let $\{v_{n}^{\prime}\}$ be the sequence obtained when applying the Laguerre transformation to the derivative of a sufficiently smooth function $v$ with $v(0)=0$, then

[TABLE]

and for the second derivative of $v$ with $v(0)=v^{\prime}(0)=0$ holds

[TABLE]

Applying the Laguerre transformation to the problem (1.1)–(1.4) with respect to the time variable, together with the relation (2.3), we obtain the following sequence of stationary boundary value problems

[TABLE]

where $n=0,1,\ldots$, $\beta_{n}=\kappa^{2}(n+1)$ and $\{u_{n}\}$ and $\{f_{n}\}$ are the Fourier-Laguerre sequences of coefficients of the functions $u$ and $f$, respectively.

By the maximum principle and induction we have the following uniqueness result.

Theorem 2.1

The sequence of stationary problems (2.5a)–(2.5c)* has at most one solution.*

The function of the form

[TABLE]

with $u_{n}$ solving (2.5a)–(2.5c), is clearly a solution to the initial boundary value problem (1.1)–(1.4). On the other hand, assuming that the solution to (1.1)–(1.4) has the right smoothness properties such that it can be expanded in time in terms of the Laguerre polynomials it follows that the coefficients will form a sequence $\{u_{n}\}$ and satisfy (2.5a)–(2.5c). Thus we state the following theorem.

Theorem 2.2

A sufficiently smooth function (2.6) is the solution of the time-dependent problem (1.1)–(1.4)* if and only if its Fourier-Laguerre coefficients $u_{n}$ for $n=0,1,\ldots$, solve the sequence of stationary problems (2.5a)–(2.5c).*

3 A boundary integral equations approach for the stationary problems

First we determine a sequence of fundamental solutions for equations (2.5a).

Definition 3.1

The sequence of $2\times 2$ matrices $\{E_{n}(x,y)\}\;n=0,1,\ldots$ is called fundamental solutions of the equations (2.5a) if

[TABLE]

Here $I$ is the $2\times 2$ identity matrix, $\delta$ denotes the Dirac function and the differentiation in (3.1) is taken with respect to $x$.

Let’s consider the polynomials

[TABLE]

for $n=0,1,\ldots,N-1$ ($w_{0}=0$), where the coefficients $a_{n,m}$ satisfy the recurrence relations

[TABLE]

and

[TABLE]

for $m=n-1,\ldots,1.$ Here $[r]$ denotes the integer part of $r\geq 0$. Next we introduce the sequence of functions

[TABLE]

where $K_{0}$ and $K_{1}$ are the modified Hankel functions of order zero and one, respectively. Throughout this paper all functions and constants with a negative index number are set equal to zero.

Lemma 3.2

The following formulas hold

[TABLE]

and

[TABLE]

where $\chi_{-2,n}=n(n-1)$, $\chi_{-1,n}=-4n^{2}$, $\chi_{0,n}=2(3n^{2}+3n+1)$, $\chi_{1,n}=-4(n+1)^{2}$ and $\chi_{2,n}=(n+1)(n+2)$.

*Proof. *Let us consider the fundamental solution for the wave equation

[TABLE]

where $\theta$ is the Heaviside function. Clearly it satisfies the equation

[TABLE]

If we apply the Laguerre transformation to this equation with respect to time, we receive the following sequence

[TABLE]

for the Laguerre coefficients

[TABLE]

In [6] it was found with the reducing of (3.5) to ordinary differential equations and using its exact solution that $G_{n}(x,y)=\Phi_{n}(\frac{\kappa}{a},|x-y|)$. Thus the formula (3.3) is proved.

The recurrence relation (2.1) gives us the following representation

[TABLE]

This relation together with (3.3) results to the formula (3.4).   $\Box$

Note here, that from (3.6) it follows $\sum_{k=-2}^{2}\chi_{k,n}=0$. Let us introduce the notation $J(x)=\frac{xx^{\top}}{|x|^{2}}$ for $x\in{\rm I\!R}^{2}\backslash\{0\}$.

Theorem 3.3

The sequence of matrices

[TABLE]

*are fundamental solutions of (2.5a).

Here*

[TABLE]

for $\ell=1,2$.

*Proof. *We consider the fundamental solution of the time-dependent elastodynamic equation (1.1) (see [1])

[TABLE]

where $r=|x-y|$. From the definition 3.1 it is clear that

[TABLE]

Thus the statement of the theorem follows from (3.8) with the use of formulas (3.3) and (3.4).   $\Box$Note that the fundamental matrix $E_{0}$ from (3.7) coincides with the fundamental matrix for the harmonic elastodynamic equation (see [1]).

Noe we can analyze the singularities in the fundamental matrix. The modified Hankel functions have the following series representations

[TABLE]

with

[TABLE]

and

[TABLE]

Here, we set $\psi(0)=0$,

[TABLE]

and let $C=0.57721\ldots$ denote Euler’s constant. Thus we can rewrite the functions $\Phi_{n}$ as follows

[TABLE]

where

[TABLE]

and

[TABLE]

Clearly we have the following asymptotic behavior with respect to $r$

[TABLE]

with

[TABLE]

and

[TABLE]

[TABLE]

Then we have the following representation for the functions in (3.7)

[TABLE]

with

[TABLE]

and

[TABLE]

Taking into account the definition of the coefficients $\chi_{k,n}$ and $a_{n,m}(\gamma)$ and following (3.10) we get the asymptotic expansion

[TABLE]

and

[TABLE]

Thus we are convinced that our fundamental sequence has only the logarithmic singularity.

We shall then construct a solution to the sequence of problems (2.5a)–(2.5c). Let $\{U_{n}\}$ be a sequence of single-layer potentials

[TABLE]

and $\{V_{n}\}$ be a sequence of double-layer potentials

[TABLE]

$n=0,1,\ldots$, where ${q_{m}}\in C(\Gamma)$ are unknown densities, $\{E_{n}\}$ is the fundamental sequence (3.7) and $T$ is a tracing operator

[TABLE]

with the unitary matrix

[TABLE]

As follows from the representation of fundamental matrices (3.7) and the expansion (3.11), the classical jump and regularity properties of the logarithmic potentials (see [9]) can be applied also to the present situation. Hence we have the following transformations into sequences of boundary integral equations.

Theorem 3.4

The sequence of single-layer potentials (3.12) is a solution of the sequence of boundary value problems (2.5a)–(2.5c)* provided that their densities satisfy the following sequence of boundary integral equations of the first kind*

[TABLE]

*for $n=0,1,\ldots$ .

The sequence of double-layer potentials (3.13) is a solution of the sequence of boundary value problems (2.5a)–(2.5c) provided that their densities satisfy the following sequence of boundary integral equations of the second kind*

[TABLE]

for $n=0,1,\ldots$ .

We proceed by investigating integral equations of the first kind (3.15). The case of the sequence of integral equations of the second kind doesn’t contain any principal different.

Theorem 3.5

For any sequence $f_{n}$ in $C^{1,\alpha}(\Gamma)$ the system (3.15) possesses a unique solution $q_{n}$ in $C^{0,\alpha}(\Gamma)$.

*Proof. *By standard arguments (see [9] for the case of the Laplace equation) it can be seen that the integral equation with logarithmic singularity

[TABLE]

has a unique solution $q_{0}\in C^{0,\alpha}(\Gamma)$ for any $f_{0}$ in $C^{1,\alpha}(\Gamma)$. Then the statement of the theorem follows by induction.   $\Box$

4 A quadrature method for full discretization

We assume that the boundary curve $\Gamma$ is given through

[TABLE]

where $x:{\rm I\!R}\to{\rm I\!R}^{2}$ is $C^{1}$ and $2\pi$–periodic with $|x^{\prime}(s)|>0$ for all $s$, such that the orientation of $\Gamma$ is counter-clockwise. Then we transform (3.15) into the parametric form

[TABLE]

where we have set $\psi_{n}(s):=|x^{\prime}(s)|\,q_{n}(x(s))$, $g_{n}(s):=f_{n}(x(s))$ and where the kernels are given by

[TABLE]

for $s\neq\tau$ and $n=0,1,\ldots.$

The kernels $H_{n}$ have logarithmic singularities and can be written in the form

[TABLE]

where

[TABLE]

and

[TABLE]

with the diagonal terms

[TABLE]

and

[TABLE]

Here we used the diagonal values for the matrix $J,$

[TABLE]

We choose $M\in{\rm I\!N}$ and an equidistant mesh by setting $s_{k}:=\frac{k\pi}{M}\,,\,k=0,\ldots,2M-1,$ and use the following quadrature rules

[TABLE]

and

[TABLE]

with the weights

[TABLE]

These quadratures are obtained by replacing the integrand $f$ by its trigonometric interpolation polynomial of degree $M$ with respect to the grid points $s_{k}$, $k=0,\ldots,2M-1.$

We use the quadrature rules (4.2)–(4.3) to approximate the integrals in the integral equations (4.1) and collocate at the nodal points to obtain the sequence of linear systems

[TABLE]

which we have to solve for the nodal values $\psi_{n,M}(s_{j})$. For the right hand sides we have

[TABLE]

For a more detailed description of this numerical solution method and an error and convergence analysis based on interpreting the above method as a fully discrete projection method in a Hölder space setting and in Sobolev space setting we refer to [9]. In particular, this error analysis implies exponential convergence

[TABLE]

for some positive constants $C_{n}$ and $\sigma,$ provided that the boundary values are also analytic. Of course due to the accumulation of the errors the constants $C_{n}$ will increase with $n$.

Given the approximate solution $\psi_{n,M}$ of the integral equation (4.1), the approximate solution of the initial boundary value problem is obtained by first evaluating the parametrized form of the potential (3.12) using the trapezoidal rule, that is, by

[TABLE]

and then summing up

[TABLE]

according to the series (2.6).

5 Numerical results

For the numerical examples we consider a kite-shaped boundary with parametrization

[TABLE]

In all examples we choose the Lamé parameters to be $\lambda=2,\mu=1$ and the density $\rho=1.$

In the first example, we set $\kappa=1,$ and we test the feasibility of the stationary problems (2.5a)–(2.5c). We choose two arbitrary points: a source point $z\in{\rm I\!R}^{2}\setminus D,$ and a measurement point $y\in D.$ We define the vector-valued boundary function

[TABLE]

where $[\cdot]_{1},$ denotes the first column of the tensor. Then, the field

[TABLE]

is clearly a solution of (2.5a)–(2.5c) for the boundary function defined above. We consider the points $z=(0.2,0.5)$ and $y=(1.5,1),$ see Figure 1. We present in Tables 1 and 2, the numerical values $\tilde{u}_{n,M}(y),$ see (4.5), and compare them with the exact solutions $u_{n}^{ex}(y),$ for $n=0,1,2$ and varying $M.$ The exponential convergence with respect to the spatial discretization is clearly exhibited, as we can we see also in Figure 2 where we plot the $L^{2}$ norm of the difference in logarithmic scale.

In the second example we consider the time-dependent problem. We set in (1.3) as boundary function

[TABLE]

Then, the exact solution is given by $u^{ex}(x,t)=[E(x,z;t)]_{1},\,x\in D,$ where its Fourier-Laguerre coefficients satisfy (2.5a)–(2.5c) for a boundary function $f_{n}$ as in the first example. We consider different source point $z=(0.4,0.2).$ We compare the computed solution $u_{N,M}(y,t),$ formula (4.6), with the exact, considering the truncated form

[TABLE]

In Table 3, we see the results of the first components, for $\kappa=1/2,$ at the position $y=(1,1),$ for different values of $M$ and $N$ at various time positions. The values of the second components, at a different position $y=(0.5,-1.5)$ are presented in Table 4.

In the third example, we consider the spatial independent boundary function

[TABLE]

which admits the expansion

[TABLE]

The numerical solution of the problem (1.1)–(1.4) is presented in Table 5 (the first component) and in Table 6 (the second component). Here, we don’t know the exact solution but we observe the convergence with respect to the discretization. We set $\kappa=1/2,$ and we compute the solution at the measurement point $y=(0.5,-1.5)$. Again we see the exponential convergence with respect to $M$ and the convergence with respect to $N$.