A roadmap for Generalized Plane Waves and their interpolation properties

TL;DR

This paper develops a systematic approach for constructing and analyzing Generalized Plane Waves (GPWs) as approximate solutions to variable coefficient PDEs, enhancing wave-based numerical methods like Discontinuous Galerkin.

Contribution

It introduces a comprehensive roadmap and algorithm for building GPWs with high-order interpolation properties for PDEs beyond the Helmholtz equation.

Findings

Established necessary conditions for high-order interpolation

Provided an explicit construction algorithm for GPWs

Extended the framework to new classes of PDEs

Abstract

This work focuses on the study of partial differential equation (PDE) based basis function for Discontinuous Galerkin methods to solve numerically wave-related boundary value problems with variable coefficients. To tackle problems with constant coefficients, wave-based methods have been widely studied in the literature: they rely on the concept of Trefftz functions, i.e. local solutions to the governing PDE, using oscillating basis functions rather than polynomial functions to represent the numerical solution. Generalized Plane Waves (GPWs) are an alternative developed to tackle problems with variable coefficients, in which case Trefftz functions are not available. In a similar way, they incorporate information on the PDE, however they are only approximate Trefftz functions since they don't solve the governing PDE exactly, but only an approximated PDE. Considering a new set of PDEs beyond the Helmholtz equation, we propose to set a roadmap for the construction and study of local interpolation properties of GPWs. Identifying carefully the various steps of the process, we provide an algorithm to summarize the construction of these functions, and establish necessary conditions to obtain high order interpolation properties of the corresponding basis.

Equations

$$-\Delta u+\beta (x,y))u=0,$$

$$\left({\partial }_x^2+(d+\overline{d}){\partial }_x{\partial }_y+|d{|}^2{\partial }_y^2\right)F+(d-\overline{d})x{\partial }_{y}F-\left(1+\frac{1}{\mu }+x(x+y)\right)F=0.$$

$${\mathcal{L}}_{M,\alpha }:=\sum _{\ell =0}^M\sum _{k=0}^{\ell }{\alpha }_{k,\ell -k}\left(x,y\right){\partial }_x^k{\partial }_y^{\ell -k}.$$

$${\mathcal{L}}_{M,\alpha }u=0\text{ on a domain }\Omega \subset {\mathbb{R}}^2,$$

$$\left\{\begin{array}[]{l}\phi(x,y)=\exp i\kappa(\cos\theta x+\sin\theta y)\\ \left[-\Delta-\kappa^{2}\right]\phi=0\end{array}\right.\text{ versus }\left\{\begin{array}[]{l}\varphi(x,y)=\exp(i\kappa(\mathbf{\cos}\theta x+\sin\theta y)+HOT)\\ \mathcal{L}_{M,\alpha}\varphi\approx 0\end{array}\right.$$

$$P(x,y)=\sum _{0\le i+j\le dP}{\lambda }_{ij}(x-x_0{)}^i(y-y_0{)}^j,$$

$$\sum _{k=0}^M{\alpha }_{k,M-k}(x_0,y_0)X^k{Y}^{M-k}=({\gamma }_1{X}^2+{\gamma }_{2}XY+{\gamma }_3{Y}^2{)}^{\frac{M}{2}}$$

$$\sum _{k=0}^M{\alpha }_{k,M-k}(x_0,y_0)X^k{Y}^{M-k}={\left({\mu }_1(A_{11}X+A_{12}Y{)}^2+{\mu }_2(A_{21}X+A_{22}Y{)}^2\right)}^{\frac{M}{2}}.$$

$$\nabla \cdot (A\nabla u)+\mathbf{d}\cdot \nabla u+k^{2}mu=0,$$

$${\mathcal{L}}_{2,\beta }{e}^{P(x,y)}=\left(-{\partial }_x^{2}P-({\partial }_{x}P{)}^2-{\partial }_y^{2}P-({\partial }_{y}P{)}^2+{\beta }_{0,0}(x,y)\right)e^{P(x,y)},$$

$$\text{Find }P\in \mathbb{C}[x,y]\text{ such that }{\partial }_x^{2}P(x,y)+({\partial }_{x}P{)}^2(x,y)+{\partial }_y^{2}P(x,y)+({\partial }_{y}P{)}^2(x,y)={\beta }_{0,0}(x,y).$$

$${\mathcal{L}}_{M,\alpha }{e}^{P(x,y)}=O(\Vert (x-x_0,y-y_0){\Vert }^q)$$

$${\mathcal{D}}^{(i,j)}:f\in {\mathcal{C}}^{i+j}\mapsto \frac{1}{i!j!}{\partial }_x^i{\partial }_y^{j}f.$$

$${\mathcal{L}}_{M,\alpha }{e}^{P(x,y)}=$$

$$\displaystyle+\sum_{\ell=1}^{M}\sum_{k=0}^{\ell}\alpha_{k,\ell-k}\left(x,y\right)k!{(\ell-k)}!\sum_{1\leq\mu\leq\ell}\sum_{s=1}^{\ell}\sum_{p_{s}((k,\ell-k),\mu)}\prod_{m=1}^{s}\frac{1}{k_{m}!}\left(\mathcal{D}^{(i_{m},j_{m})}P(x,y)\right)^{k_{m}}\Bigg{)},$$

$$\forall(\mu,\nu)\in(\mathbb{N}^{2})^{2},\mu\prec\nu\Leftrightarrow\begin{array}[]{l}1.\ \mu_{1}+\mu_{2}<\nu_{1}+\nu_{2};\text{ or }\\ 2.\ \mu_{1}+\mu_{2}=\nu_{1}+\nu_{2}\text{ and }\mu_{1}<\mu_{2},\end{array}$$

$$\begin{array}[]{l}\Bigg{\{}(k_{1},\cdots,k_{s};(i_{1},j_{1}),\cdots,(i_{s},j_{s})):k_{i}>0,0\prec(i_{1},j_{1})\prec\cdots\prec(i_{s},j_{s}),\\ \displaystyle\phantom{\Bigg{\{}}\sum_{l=1}^{s}k_{l}=\mu,\sum_{l=1}^{s}k_{l}i_{l}=i,\sum_{l=1}^{s}k_{l}j_{l}=j\Bigg{\}}.\end{array}$$

$${\mathfrak{L}}_{2,\gamma }{e}^P=$$

$$\displaystyle\phantom{e^{P}\Bigg{(}}+\gamma_{1,0}\partial_{x}P+\gamma_{0,1}\partial_{y}P+\gamma_{0,0}\Bigg{)}.$$

$${\mathcal{L}}_{M,\alpha }^A=\sum _{\ell =1}^M\sum _{k=0}^{\ell }k!(\ell -k)!{\alpha }_{k,\ell -k}\sum _{1\le \mu \le \ell }\sum _{s=1}^{\ell }\sum _{p_s((k,\ell -k),\mu )}\prod _{m=1}^s\frac{1}{k_m!}{\left({\mathcal{D}}^{(i_m,j_m)}(\cdot )\right)}^{k_m},$$

$${\mathcal{L}}_{M,\alpha }^A:f\in {\mathcal{C}}^M\mapsto \frac{{\mathcal{L}}_{M,\alpha }{e}^f}{e^f}-{\alpha }_{0,0}.$$

$${\mathfrak{L}}_{2,\gamma }^{A}P=-{\partial }_x^{2}P+{\gamma }_{1,1}{\partial }_x{\partial }_{y}P+{\gamma }_{0,2}{\partial }_y^{2}P-({\partial }_{x}P{)}^2+{\gamma }_{1,1}{\partial }_{x}P{\partial }_{y}P+{\gamma }_{0,2}({\partial }_{y}P{)}^2+{\gamma }_{1,0}{\partial }_{x}P+{\gamma }_{0,1}{\partial }_{y}P.$$

$${\mathcal{L}}_{M,\alpha }^{A}P(x,y)=-{\alpha }_{0,0}(x,y)+O(|(x-x_0,y-y_0){|}^q).$$

$$\text{Find }P\in \mathbb{C}[x,y],\text{ s.t. }\forall (I,J)\in {\mathbb{N}}^2,0\le I+J<q,{\mathcal{D}}^{(I,J)}{\mathcal{L}}_{M,\alpha }^{A}P(x_0,y_0)=-{\mathcal{D}}^{(I,J)}{\alpha }_{0,0}(x_0,y_0).$$

$${\mathcal{L}}_{M,\alpha }^L=\sum _{\ell =1}^M\sum _{k=0}^{\ell }{\alpha }_{k,\ell -k}{\partial }_x^k{\partial }_y^{\ell -k},$$

$${\mathcal{L}}_{M,\alpha }^N=\sum _{\ell =1}^M\sum _{k=0}^{\ell }k!(\ell -k)!{\alpha }_{k,\ell -k}\sum _{2\le \mu \le \ell }\sum _{s=1}^{\ell }\sum _{p_s((k,\ell -k),\mu )}\prod _{m=1}^s\frac{1}{k_m!}{\left({\mathcal{D}}^{(i_m,j_m)}(\cdot )\right)}^{k_m}.$$

$$\left\{\begin{array}[]{l}\displaystyle\mathfrak{L}_{2,\gamma}^{L}P=-\partial_{x}^{2}P+\gamma_{1,1}\partial_{x}\partial_{y}P+\gamma_{0,2}\partial_{y}^{2}P+\gamma_{1,0}\partial_{x}P+\gamma_{0,1}\partial_{y}P,\\ \displaystyle\mathfrak{L}_{2,\gamma}^{N}P=-(\partial_{x}P)^{2}+\gamma_{1,1}\partial_{x}P\partial_{y}P+\gamma_{0,2}(\partial_{y}P)^{2}.\end{array}\right.$$

$${\mathcal{D}}^{(I,J)}\left[{\mathcal{L}}_{M,\alpha }^{L}P\right](x_0,y_0)$$

$$\begin{array}[]{l}\displaystyle\mathcal{D}^{(I,J)}\left[\mathcal{L}_{M,\alpha}^{L}P\right](x_{0},y_{0})\\ \displaystyle\phantom{=}=\frac{1}{I!J!}\sum_{k=0}^{M}\alpha_{k,{M}-k}(x_{0},y_{0})\partial_{x}^{k+I}\partial_{y}^{{M}-k+J}P(x_{0},y_{0})\\ \phantom{==}\displaystyle+\sum_{k=0}^{M}\sum_{\tilde{i}=0}^{I-1}\sum_{\tilde{j}=0}^{J-1}\frac{1}{\tilde{i}!\tilde{j}!}\mathcal{D}^{(I-\tilde{i},J-\tilde{j})}\alpha_{k,{M}-k}(x_{0},y_{0})\partial_{x}^{k+\tilde{i}}\partial_{y}^{{M}-k+\tilde{j}}P(x_{0},y_{0})\\ \phantom{==}\displaystyle+\sum_{\ell=1}^{M-1}\sum_{k=0}^{\ell}\sum_{\tilde{i}=0}^{I}\sum_{\tilde{j}=0}^{J}\frac{1}{\tilde{i}!\tilde{j}!}\mathcal{D}^{(I-\tilde{i},J-\tilde{j})}\alpha_{k,\ell-k}(x_{0},y_{0})\partial_{x}^{k+\tilde{i}}\partial_{y}^{\ell-k+\tilde{j}}P(x_{0},y_{0}).\end{array}$$

$$\forall(i,j)\in\mathbb{N}^{2},\mathcal{D}^{(i,j)}P(x_{0},y_{0})=\left\{\begin{array}[]{ll}\lambda_{i,j}&\text{ if }i+j\leq dP,\\ 0&\text{otherwise},\end{array}\right.$$

$$\begin{array}[]{c}\displaystyle\text{Find }\{\lambda_{i,j},(i,j)\in\mathbb{N}^{2},0\leq i+j\leq M+q-1\}\text{ such that }\hfill\\ \displaystyle P(x,y):=\sum_{i=0}^{M+q-1}\sum_{j=0}^{M+q-1-i}\lambda_{i,j}(x-x_{0})^{i}(y-y_{0})^{j}\in\mathbb{C}[x,y],\text{ satisfies }\\ \displaystyle\forall(I,J)\in\mathbb{N}^{2},0\leq I+J<q,\mathcal{D}^{(I,J)}\mathcal{L}_{M,\alpha}^{A}P(x_{0},y_{0})=-\mathcal{D}^{(I,J)}\alpha_{0,0}(x_{0},y_{0}).\end{array}$$

Full text

1 Introduction

Trefftz methods are Galerkin type of methods that rely on function spaces of local solutions to the governing partial differential equations (PDEs). They were initially introduced in [34, 26], and the original idea was to use trial functions which satisfy the governing PDE to derive error bounds. They have been widely used in the engineering community [18] since the 60’s, for instance for Laplace’s equation [29], to the biharmonic equation [32] and to elasticity [21]. Later the general idea of taking advantage of analytical knowledge about the problem to build a good approximation space was used to develop numerical methods: in the presence of corner and interface singularities [9, 33], boundary layers, rough coefficients, elastic interactions [28, 1, 27, 2], wave propagation [1, 8]. In the context of boundary value problems (BVPs) for time-harmonic wave propagation, several methods have been proposed following the idea of functions that solves the governing PDE, [22], relying on incorporating oscillating functions in the function spaces to derive and discretize a weak formulation. Wave-based numerical methods have received attention from several research groups around the world, from the theoretical [22] and computational [12] point of view, and the pollution effect of plane wave Discontinuous Galerkin (DG) methods was studied in [10]. Such methods have also been implemented in industry codes111http://www.waveller.com/Waveller_Acoustics/waveller_acoustics.shtml, for acoustic applications. The use of Plane Wave (PW) basis functions has been the most popular choice, while an attempt to use Bessel functions was reported in [24]. In [23], the authors present an interesting comparison of performance between high order polynomial and wave-based methods. More recently, application to space-time problems have been studied in [31, 3, 20, 30, 19].

In this context, numerical methods rely on discretizing a weak formulation via a set of exact solutions of the governing PDE. When no exact solutions to the governing PDE are available, there is no natural choice of basis functions to discretize the weak formulation. This is in particular the case for variable coefficient problems. In order to take advantage of Trefftz-type methods for problems with variable coefficients, Generalized Plane Waves (GPWs) were introduced in [15], as basis functions that are local approximate solutions - rather than exact solutions - to the governing PDE. GPWs were designed adding higher order terms in the phase of classical PWs, choosing these higher order terms to ensure the desired approximation of the governing PDE. In [13], the construction and interpolation properties of GPWs were studied for the Helmholtz equation

[TABLE]

with a particular interest for the case of a sign-changing coefficient $\beta$, including propagating solutions ($\beta<0$), evanescent solutions ($\beta>0$), smooth transition between them ($\beta=0$) called cut-offs in the field of plasma waves. The interpolation properties of a set $\mathbb{V}$ spanned by resulting basis functions, namely $\|(I-P_{\mathbb{V}})u\|$ where $P_{\mathbb{V}}$ is the orthogonal projector on $\mathbb{V}$ while $u$ is the solution to the original problem, play a crucial role in the error estimation of the corresponding numerical method [4]. For this same equation the error analysis of a modified Trefftz method discretized with GPWs was presented in [16]. In [17], Generalized Plane Waves (GPWs) were used for the numerical simulation of mode conversion modeled by the following equation:

[TABLE]

In the present work, we answer questions related to extending the work on GPW developed in [13] - the construction of GPWs on the one hand, and their interpolation properties on the other hand - from the Helmholtz operator $-\Delta+\beta$ to a wide range of partial differential operators. A construction process valid for some operators of order two or higher is presented, while a proof of interpolation properties is limited to some operators of order two. We propose a road map to identify crucial steps in our work:

Construction of GPWs $\varphi$ such that $\mathcal{L}\varphi\approx 0$

1. (a)

Choose an ansatz for $\varphi$ (Section 2). 2. (b)

Identify the corresponding $N_{dof}$ degrees of freedom and $N_{eqn}$ constraints (Subsection 2.1). 3. (c)

Choose the number of degrees of freedom adequately $N_{dof}\geq N_{eqn}$ (Subsection 2.1). 4. (d)

Study the structure of the resulting system and identify $N_{dof}-N_{eqn}$ additional constraints (Subsections 2.2 and 2.3). 5. (e)

Compute the remaining $N_{eqn}$ degrees of freedom at minimal computational cost (Subsection 2.4). 2. 2.

Interpolation properties

1. (a)

Study the properties of the remaining $N_{eqn}$ degrees of freedom with respect to the $N_{dof}-N_{eqn}$ additional constraints 2. (b)

Identify a simple reference case depending only on the $N_{dof}-N_{eqn}$ additional constraints (Section 3). 3. (c)

Study the interpolation properties of this reference case (Subsection 4.1). 4. (d)

Relate the general case to the reference case (Subsections 3.1 and 3.2). 5. (e)

Prove the interpolation properties of the GPWs from those of the reference case (Subsection 4.2).

We will consider linear partial differential operators with variable coefficients, defined as follows.

Definition 1**.**

A linear partial differential operator of order $M\geq 2$, in two dimensions, with a given set of complex-valued coefficients $\alpha=\{\alpha_{k,\ell-k},(k,\ell)\in\mathbb{N}^{2},0\leq k\leq\ell\leq M\}$ will be denoted hereafter as

[TABLE]

Our goal is to build a basis of functions well suited to approximate locally any solution $u$ to a given homogeneous variable-coefficient partial differential equation

[TABLE]

where by locally we mean piecewise on a mesh $\mathcal{T}_{h}$ of $\Omega$. Such interpolation properties are a building block for the convergence proof of Galerkin methods. For a constant coefficient operator, it is natural to use the same basis on each element $K\in\mathcal{T}_{h}$. However, with variable coefficients, it cannot be optimal to expect a single basis to have good approximation properties on the whole domain $\Omega\subset\mathbb{R}^{2}$. For instance, for the Helmholtz equation with a sign-changing coefficient, it can not be optimal to look for a single basis that would give a good approximation of solutions both in the propagating region and in the evanescent region. Therefore it is natural to think of local bases defined on each $K\in\mathcal{T}_{h}$: with GPWs we focus on local properties around a given point $(x_{0},y_{0})\in\mathbb{R}^{2}$ rather than on a given domain $\Omega$. A simple idea would then be freezing the coefficients of the operator, that is to say studying, instead of $\mathcal{L}_{M,\alpha}$, the constant coefficient operator $\mathcal{L}_{M,\bar{\alpha}}$ with constant coefficients $\bar{\alpha}=\{\alpha_{k,l}(x_{0},y_{0}),0\leq k+l\leq M\}$. However, as observed in [13, 14], this leads to low order approximation properties, while we are interested in high order approximation properties. This is why new functions are needed to handle variable coefficients. This work will focus on two aspects: the construction and the interpolation properties of GPWs.

We follow the GPW design proposed in [13, 15]. Retaining the oscillating feature while aiming for higher order approximation, GPW were designed with Higher Order Terms ($HOT$) in the phase function of a plane wave. These higher order terms are to be defined to ensure that a GPW function $\varphi$ is an approximate solution to the PDE:

[TABLE]

In Section 2, the construction of a GPW $\varphi(x,y)=e^{P(x,y)}$ will be described in detail, then a precise definition of GPW will be provided under the following hypothesis:

Hypothesis 1**.**

Consider a given point $(x_{0},y_{0})\in\mathbb{R}^{2}$, a given approximation parameter $q\in\mathbb{N}$, $q\geq 1$, a given $M\in\mathbb{N}$, $M\geq 2$, and a partial differential operator $\mathcal{L}_{M,\alpha}$ defined by a given set of complex-valued coefficients $\alpha=\{\alpha_{k,l},0\leq k+l\leq M\}$, defined in a neighborhood of $(x_{0},y_{0})$, satisfying

• •

$\alpha_{k,l}$* is $\mathcal{C}^{q-1}$ at $(x_{0},y_{0})$ for all $(k,l)$ such that $0\leq k+l\leq M$,*

• •

$\alpha_{M,0}(x_{0},y_{0})\neq 0$.

This construction is equivalent to the construction of the bi-variate polynomial

[TABLE]

and is performed by choosing the degree $dP$, and providing an explicit formula for the set of complex coefficients $\{\lambda_{ij}\}_{\{(i,j)\in\mathbb{N}^{2},0\leq i+j\leq dP\}}$, in order for $\varphi$ to satisfy $\mathcal{L}_{M,\alpha}\varphi(x,y)=O\left(\|(x,y)-(x_{0},y_{0})\|^{q}\right)$. An algorithm to construct a GPW is provided. In Section 3 properties of the $\lambda_{ij}$s are studied, while the interpolation properties of the corresponding set of basis functions are studied for the case $M=2$ in Section 4, under the following hypothesis:

Hypothesis 2**.**

Under Hypothesis 1 we consider only operators $\mathcal{L}_{M,\alpha}$ such that $M$ is even and the terms of order $M$ satisfy

[TABLE]

for some complex numbers $(\gamma_{1},\gamma_{2},\gamma_{3})$ such that there exists $(\mu_{1},\mu_{2})\in\mathbb{C}^{2}$, $\mu_{1}\mu_{2}\neq 0$, a non-singular matrix $A\in\mathbb{C}^{2\times 2}$ satisfying $\Gamma=A^{t}DA$ where $\Gamma=\begin{pmatrix}\gamma_{1}&\gamma_{2}/2\\ \gamma_{2}/2&\gamma_{3}\end{pmatrix}$ and $D=\begin{pmatrix}\mu_{1}&0\\ 0&\mu_{2}\end{pmatrix}$ , and therefore

[TABLE]

For instance, these matrices are $\Gamma=D=Id$ for $\mathcal{L}_{H}:=-\Delta-\kappa^{2}(x,y)$ or $\mathcal{L}_{B}:=\Delta\mathcal{L}_{H}$, and $\Gamma=D=c(x_{0},y_{0})Id$ for $\mathcal{L}_{C}:=-\nabla\cdot(c(x,y)\nabla)-\kappa^{2}(x,y)$. Note that if $\Gamma$ is real, this is simply saying that its eigenvalues are non-zero. Finally, corresponding numerical results are then provided, for various operators $\mathcal{L}_{M,\alpha}$ of order $M=2$ in Section 5.

Our previous work was limited to the Helmholtz equation (1) for propagating and evanescent regions, transition between the two, absorbing regions, as well as caustics. The interpolation properties presented here cover more general second order equations, in particular equations that can be written as

[TABLE]

with variable coefficients $A$ matrix-valued, real and symmetric with non-zero eigenvalues, $\mathbf{d}$ vector-valued and $m$ scalar-valued. It includes for instance

• •

Helmholtz equation with absorption corresponding to $A=I$ with $\Re(m)>0$ and $\Im(m)\neq 0$ ;

• •

the mild-slop equation [7] modeling the amplitude of the free-surface water waves corresponding to $m=c_{p}c_{g}$ being the product of $c_{p}$ the phase speed of the waves and $c_{g}$ the group speed of the waves with $A=mId$ ;

• •

if $\mu$ is the permeability and $\epsilon$ the permittivity, then the transverse-magnetic mode of Maxwell’s equations for $A=\frac{1}{\mu}I$ and $m=\epsilon$, while the transverse-electric mode of Maxwell’s equations for $A=\frac{1}{\epsilon}I$ and $m=\mu$.

Throughout this article, we will denote by $\mathbb{N}$ the set of non-negative integers, by $\mathbb{N}^{*}$ the set of positive integers, by $\mathbb{R}^{+}=[0;+\infty)$ the set of non-negative real numbers, and by $\mathbb{C}[z_{1},z_{2}]$ the space of complex polynomials with respect to the two variables $z_{1}$ and $z_{2}$. As the first part of this work is dedicated to finding the coefficients $\lambda_{ij}$, we will reserve the word unknown to refer to the $\lambda_{i,j}$s. The length of the multi-index $(i,j)\in\mathbb{N}^{2}$ of an unknown $\lambda_{ij}$, $|(i,j)|=i+j$, will play a crucial role in what follows.

2 Construction of a GPW

The task of constructing a GPW is attached to a homogeneous PDE, it is not global on $\mathbb{R}^{2}$ but it is local as it is expressed in terms of a Taylor expansion. It consists in finding a polynomial $P\in\mathbb{C}[x,y]$ such that the corresponding GPW, namely $\varphi:=e^{P}$, is locally an approximate solution to the PDE.

Consider $M=2$, $\beta=\{\beta_{0,0},\beta_{0,1}=\beta_{1,0}=\beta_{1,1}=0,\beta_{0,2}\equiv-1,\beta_{2,0}\equiv-1\}$, and the corresponding the operator $\mathcal{L}_{2,\beta}=-\partial_{x}^{2}-\partial_{y}^{2}+\beta_{0,0}(x)$. Then for any polynomial $P\in\mathbb{C}[x,y]$:

[TABLE]

so the construction of an exact solution to the PDE would be equivalent to the following problem:

[TABLE]

Consider then the following examples.

• •

If $\beta_{0,0}(x,y)$ is constant, then it is straightforward to find a polynomial of degree one satisfying Problem (5); $\beta_{0,0}$ being negative this would correspond to a classical plane wave.

• •

If $\beta_{0,0}(x,y)=x$, then there is no solution to (5), since the total degree of $\partial_{x}^{2}P+(\partial_{x}P)^{2}+\partial_{y}^{2}P+(\partial_{y}P)^{2}$ is always even.

• •

If $\beta_{0,0}(x,y)$ is not a polynomial function, it is also straightforward to see that no polynomial $P$ can satisfy Problem (5).

From these trivial examples we see that in general there is no such function, $\varphi(x,y)=e^{P(x,y)}$, $P$ being a complex polynomial, solution to a variable coefficient partial differential equation exactly. It could seem that the restriction for $P$ to be a polynomial is very strong. However since we are interested in approximation and smooth coefficients, rather than looking for a more general phase function we restrict the identity $\mathcal{L}\varphi=0$ on $\Omega$ into an approximation on a neighborhood of $(x_{0},y_{0})\in\mathbb{R}^{2}$ in the following sense. We replace the too restrictive cancellation of $\mathcal{L}_{M,\alpha}e^{P(x,y)}$ by the cancellation of the lowest terms of its Taylor expansion around $(x_{0},y_{0})$. So this section is dedicated to the construction of a polynomial $P\in\mathbb{C}[x,y]$, under Hypothesis 1, to ensure that the following local approximation property

[TABLE]

is satisfied. The parameter $q$ will denote throughout this work the order of approximation of the equation to which the GPW is designed. In summary, the construction is performed:

• •

for a partial differential operator $\mathcal{L}_{M,\alpha}$ of order $M$ defined by a set of smooth coefficients $\alpha$,

• •

at a point $(x_{0},y_{0})\in\mathbb{R}^{2}$,

• •

at order $q\in\mathbb{N}^{*}$,

• •

to ensure that $\mathcal{L}_{M,\alpha}e^{P(x,y)}=O(|(x-x_{0},y-y_{0})|^{q}).$

Even though the construction of a GPW will involve a non-linear system we propose to take advantage of the structure of this system to construct a solution via an explicit formula. In this way, even though a GPW $\varphi:=e^{P}$ is a PDE-based function, the polynomial $P$ can be constructed in practice from this formula, and therefore the function can be constructed without solving numerically any non-linear - or even linear - system. This remark is of great interest with respect to the use of such functions in a Discontinuous Galerkin method to solve numerically boundary value problems.

In order to illustrate the general formulas that will appear in this section, we will use the specific case $\mathfrak{L}_{2,\gamma}$ where $\gamma=\{\gamma_{0,0},\gamma_{1,0},\gamma_{0,1},\gamma_{2,0}\equiv-1,\gamma_{1,1},\gamma_{0,2}\}$, for which we can write explicitly many formulas is a compact form. In order to simplify certain expressions that will follow we propose the following definition.

Definition 2**.**

Assume $(i,j)\in\mathbb{N}^{2}$ and $(x_{0},y_{0})\in\mathbb{R}^{2}$. We define the linear partial differential operator $D^{(i,j)}$ by

[TABLE]

A precise definition of GPW will be provided at the end of this section.

2.1 From the Taylor expansion to a non-linear system

We are seeking a polynomial $\displaystyle P(x,y)=\sum_{0\leq i+j\leq dP}\lambda_{ij}(x-x_{0})^{i}(y-y_{0})^{j}$ satisfying the Taylor expansion (6). Defining such a polynomial is equivalent to defining the set $\{\lambda_{ij};(i,j)\in\mathbb{N}^{2},0\leq i+j\leq dP\}$, and therefore we will refer to the $\lambda_{ij}$s as the unknowns throughout this construction process. The goal of this subsection is to identify the set of equations to be satisfied by these unknowns to ensure that $P$ satisfies the Taylor expansion (6), and in particular to choose the degree of $P$ so as to guarantee the presence of linear terms in each equation of the system.

According to the Faa di Bruno formula, the action of the partial differential operator $\mathcal{L}_{M,\alpha}$ on a function $\varphi(x,y)=e^{P(x,y)}$ is given by

[TABLE]

where the linear order $\prec$ on $\mathbb{N}^{2}$ is defined by

[TABLE]

and where $p_{s}((i,j),\mu)$ is equal to

[TABLE]

For the operator $\mathfrak{L}_{2,\gamma}$ the Faa di Bruno formula becomes

[TABLE]

In order to single out the terms depending on $P$ in the right hand side, this leads to the following definition.

Definition 3**.**

Consider a given $M\in\mathbb{N}$, $M\geq 2$, a given set of complex-valued functions $\alpha=\{\alpha_{k,l},0\leq k+l\leq M\}$, and the corresponding partial differential operator $\mathcal{L}_{M,\alpha}$. We define the partial differential operator $\mathcal{L}_{M,\alpha}^{A}$ associated to $\mathcal{L}_{M,\alpha}$ as

[TABLE]

or equivalently, since the exponential of a bounded quantity is bounded away from zero:

[TABLE]

For the operator $\mathfrak{L}_{2,\gamma}$ this gives

[TABLE]

Since, for any polynomial $P$, the function $e^{P}$ is locally bounded, and since {\mathcal{L}_{M,\alpha}[e^{P}]}=\big{(}\mathcal{L}_{M,\alpha}^{A}{e^{P}}+\alpha_{0,0}\big{)}e^{P}, then for a polynomial $P$ to satisfy the approximation property (6), it is sufficient to satisfy

[TABLE]

Therefore, the problem to be solved is now:

[TABLE]

In order to define a polynomial $\displaystyle P(x,y)=\sum_{0\leq i+j\leq dP}\lambda_{ij}(x-x_{0})^{i}(y-y_{0})^{j}$, the degree $dP$ of the polynomial determines the number of unknowns: there are $N_{dof}=\frac{(dP+1)(dP+2)}{2}$ unknowns to be defined, namely the $\{\lambda_{i,j}\}_{\{(i,j)\in\mathbb{N},0\leq i+j\leq dP}$. In order to design a polynomial $P$ satisfying Equation (7), the parameter $q$ determines the number of equations to be solved: there are $N_{eqn}=\frac{q(q+1)}{2}$ terms to be canceled from the Taylor expansion. The first step toward the construction of a GPW is to define the value of $dP$ for a given value of $q$.

At this point it is clear that if $dP\leq q-1$, then the resulting system is over-determined. Our choice for the polynomial degree $dP$ relies on a careful examination of the linear terms in $\mathcal{L}_{M,\alpha}^{A}P$. We can already notice that, under Hypothesis 1, in $\mathcal{L}_{M,\alpha}^{A}P$ there is at least one non-zero linear term, namely $\alpha_{M,0}(x_{0},y_{0})\partial_{x}^{M}P$, and there is at least one non-zero non-linear term, namely $\alpha_{M,0}(x_{0},y_{0})(\partial_{x}P)^{M}$. This non-linear term corresponds to the following parameters from the Faa di Bruno formula: $\mu=M$, $s=1$, $(k_{1},(i_{1},j_{1}))=(M,(1,0))$. The linear terms can only correspond to $s=1$, $\mu=1$ and $p_{1}((k,\ell-k),1)=\{(1,(k,\ell-k))\}$, see Definition 3. We can then split $\mathcal{L}_{M,\alpha}^{A}$ into its linear and non-linear parts.

Definition 4**.**

Consider a given $M\in\mathbb{N}$, $M\geq 2$, a given set of complex-valued functions $\alpha=\{\alpha_{k,l},0\leq k+l\leq M\}$, and the corresponding partial differential operator $\mathcal{L}_{M,\alpha}$. The linear part of the partial differential operator $\mathcal{L}_{M,\alpha}^{A}$ is defined by $\mathcal{L}_{M,\alpha}^{L}:=\mathcal{L}_{M,\alpha}-\alpha_{0,0}\partial_{x}^{0}\partial_{y}^{0}$, or equivalently

[TABLE]

and its non-linear part $\mathcal{L}_{M,\alpha}^{N}:=\mathcal{L}_{M,\alpha}^{A}-\mathcal{L}_{M,\alpha}^{L}$ can equivalently be defined by

[TABLE]

For the operator $\mathfrak{L}_{2,\gamma}$ this gives respectively

[TABLE]

Consider the $(I,J)$ coefficients of the Taylor expansion of $\mathcal{L}_{M,\alpha}^{L}P$ for $(I,J)\in\mathbb{N}^{2}$ and $0\leq I+J<q$:

[TABLE]

so that in order to isolate the derivatives of highest order, i.e. of order $M+I+J$, we can write

[TABLE]

Back to Problem (8), the $(I,J)$ terms (9) a priori depend on the unknowns $\{\lambda_{i,j},(i,j)\in\mathbb{N}^{2},0\leq i+j\leq dP\}$. Since

[TABLE]

then under Hypothesis 1 any $(I,J)$ term in System (8) has at least one non-zero linear term, as long as $I+J\leq dP-M$, namely $\frac{(M+I)!}{I!}\alpha_{M,0}(x_{0},y_{0})\lambda_{M+I,J}$, while it does not necessarily have any linear term as soon as $I+J>dP-M$. Avoiding equations with no linear terms is natural, and it will be crucial for the construction process described hereafter.

Choosing the polynomial degree to be $dP=M+q-1$ therefore guarantees the existence of at least one linear term in every equation of System (8). Therefore, from now on the polynomial $P$ will be of degree $dP=M+q-1$ and the new problem to be solved is

[TABLE]

As a consequence the number of unknowns is $N_{dof}=\frac{(M+q)(M+q+1)}{2}$, and the system is under-determined : $N_{dof}-N_{eqn}=Mq+\frac{M(M+1)}{2}$. See Figure 1 for an illustration of the equation and unknown count.

Note that this system is always non-linear. Indeed, under Hypothesis 1, the $(0,0)$ equation of the system always includes the non-zero non-linear term $\alpha_{M,0}(x_{0},y_{0})(\lambda_{1,0})^{M}$, corresponding to the following parameters from the Faa di Bruno formula: $\mu=M$, $s=1$, $(k_{1},(i_{1},j_{1}))=(M,(1,0))$.

The key to the construction procedure proposed next is a meticulous gathering of unknowns $\lambda_{i,j}$ with respect the length of their multi-index $i+j$. As we will now see, this will lead to splitting the system into a hierarchy of simple linear sub-systems.

2.2 From a non-linear system to linear sub-systems

The different unknowns appearing in each equation of System (10) can now be studied. A careful inspection of the linear and non-linear terms will reveal the underlying structure of the system, and will lead to identify a hierarchy of simple linear subsystems.

The inspection of the linear terms is very straightforward thanks to Equation (9). The description of the unknowns in the linear terms is summarized here.

Lemma 1**.**

Consider a point $(x_{0},y_{0})\in\mathbb{R}^{2}$, a given $q\in\mathbb{N}^{*}$, a given $M\in\mathbb{N}$, $M\geq 2$, a given set of complex-valued functions $\alpha=\{\alpha_{k,l}\in\mathcal{C}^{q-1}\text{ at }(x_{0},y_{0}),0\leq k+l\leq M\}$, and the corresponding partial differential operator $\mathcal{L}_{M,\alpha}$. In each equation $(I,J)$ of System (10), the linear terms can be split as follows:

• •

a set of unknowns with length of the multi-index equal to $M+I+J$, corresponding to $\ell=M$ and $(\tilde{i},\tilde{j})=(I,J)$,

• •

a set of unknowns with length of the multi-index at most equal to $M+I+J-1$.

Under Hypothesis 1, both sets are never empty.

Proof.

In terms of unknowns $\{\lambda_{i,j},(i,j)\in\mathbb{N}^{2},0\leq i+j\leq M+q-1\}$, Equation (9) reads :

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The result is immediate for $I=J=0$ from (11). The following comments are valid for the right hand sides of (12), (13), and (14): the third term only contains unknowns with a length of the multi-index equal to $\ell+\tilde{i}+\tilde{j}\leq M-1+I+J$, while the second term only contains unknowns with a length of the multi-index equal to $M+\tilde{i}+\tilde{j}\leq M+I+J-2$ ; as to the first term, it only contains unknowns with a length of the multi-index equal to $M+I+J$. This proves the claim. ∎

We then focus on the inspection of the non-linear terms. Each non-linear term in $\mathcal{L}_{M,\alpha}^{A}P$ reads from the definition of $\mathcal{L}_{M,\alpha}^{N}$

[TABLE]

and yields a sum of non-linear terms with respect to the unknowns $\{\lambda_{ij}\}_{\{(i,j),0\leq i+j\leq M+q-1\}}$, implicitly given by the following formula:

[TABLE]

Therefore coming from the term (15), only a restricted number of unknowns contribute to the $(I,J)$ equation of Problem (10).

In order to identify the unknowns contributing to (16), here are two simple yet important reminders are provided in Appendix C.

It is now straightforward to describe the unknowns $\lambda_{i,j}$ appearing in the non-linear terms of the equation $(I,J)$ of System (10), unwinding formula (16).

Lemma 2**.**

Consider a point $(x_{0},y_{0})\in\mathbb{R}^{2}$, a given $q\in\mathbb{N}^{*}$, a given $M\in\mathbb{N}$, $M\geq 2$, a given set of complex-valued functions $\alpha=\{\alpha_{k,l}\in\mathcal{C}^{q-1}\text{ at }(x_{0},y_{0}),0\leq k+l\leq M\}$, and the corresponding partial differential operator $\mathcal{L}_{M,\alpha}$. In each equation $(I,J)$ of System (10), the unknowns $\lambda_{i,j}$ appearing in the non-linear terms have a length of the multi-index $i+j<M+I+J$.

Proof.

Each term $\partial_{x}^{\tilde{i}}\partial_{y}^{\tilde{j}}\left[\prod_{m=1}^{s}\left(\partial_{x}^{i_{m}}\partial_{y}^{j_{m}}P\right)^{k_{m}}\right]$ in $\mathcal{L}_{M,\alpha}^{A}P$ is a polynomial, and its constant coefficient contains coefficients of the polynomial $\prod_{m=1}^{s}\left(\partial_{x}^{i_{m}}\partial_{y}^{j_{m}}P\right)^{k_{m}}$ with a length of the multi-index length of the multi-index at most equal to $\tilde{i}+\tilde{j}$, that is to say coefficients of the polynomials $\partial_{x}^{i_{m}}\partial_{y}^{j_{m}}P$ with a length of the multi-index length of the multi-index at most equal to $\tilde{i}+\tilde{j}$ for every $(i_{m},j_{m})$ from the Faa di Bruno’s formula, so coefficients $\lambda_{i,j}$ of the polynomial $P$ with a length of the multi-index at most equal to $\tilde{i}+\tilde{j}+i_{m}+j_{m}$. Since the indices are such that $\tilde{i}\leq I$, $\tilde{j}\leq J$, and $i_{m}+j_{m}\leq\ell<M$, the unknowns $\lambda_{i,j}$ appearing in each term $\partial_{x}^{\tilde{i}}\partial_{y}^{\tilde{j}}\left[\prod_{m=1}^{s}\left(\partial_{x}^{i_{m}}\partial_{y}^{j_{m}}P\right)^{k_{m}}\right](x_{0},y_{0})$ have a length of the multi-index at most equal to $M+I+J-1$. It is therefore true for any linear combination such as (16). ∎

From the two previous Lemmas, we see that, in each equation $(I,J)$ of System (10), unknowns with a length of the multi-index equal to $M+I+J$ appear only in linear terms, namely in

[TABLE]

whereas all the remaining unknowns have a length of the multi-index at most equal to $M+I+J-1$. It is consequently natural to subdivide the set of unknowns with respect to the length of their multi-index $M+\mathfrak{L}$, for $\mathfrak{L}$ between [math] and $q-1$ in order to take advantage of this linear structure.

2.3 Hierarchy of triangular linear systems

Our goal is now to construct a solution to the non-linear system (10), and our understanding of its linear part will lead to an explicit construction of such a solution without any need for any approximation.

The crucial point of our construction process is to take advantage of the underlying layer structure with respect to the length of the multi-index: it is only natural now to gather into subsystems all equations $(I,\mathfrak{L}-I)$ for $I$ between [math] and $\mathfrak{L}$, while gathering similarly all unknowns with length of the multi-index equal to $M+\mathfrak{L}$. In the subsystem of layer $\mathfrak{L}$, we know that the unknowns with a length of the multi-index equal to $M+I+J$ only appear in linear terms, and we rewrite each equation $(I,J)$ as

[TABLE]

For the sake of clarity, the resulting right-hand side terms can defined as follows.

Definition 5**.**

Consider a point $(x_{0},y_{0})\in\mathbb{R}^{2}$, a given $q\in\mathbb{N}^{*}$, a given $M\in\mathbb{N}$, $M\geq 2$, a given set of complex-valued functions $\alpha=\{\alpha_{k,l}\in\mathcal{C}^{q-1}\text{ at }(x_{0},y_{0}),0\leq k+l\leq M\}$, and the corresponding partial differential operator $\mathcal{L}_{M,\alpha}$. We define the quantity $N_{I,J}$ from Equation $(I,J)$ from (10) as

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[EX] For the operator $\mathfrak{L}_{2,\gamma}$ the non-linear terms in $N_{0,0}$, $N_{1,0}$ and $N_{0,1}$ are respectively

[TABLE]

[TABLE]

[TABLE]

We now consider the following subsystems for given $\mathfrak{L}$ between [math] and $q-1$:

[TABLE]

The layer structure follows from our understanding of the non-linearity of the original system:

Corollary 1**.**

Consider a point $(x_{0},y_{0})\in\mathbb{R}^{2}$, a given $q\in\mathbb{N}^{*}$, a given $M\in\mathbb{N}$, $M\geq 2$, a given set of complex-valued functions $\alpha=\{\alpha_{k,l}\in\mathcal{C}^{q-1}\text{ at }(x_{0},y_{0}),0\leq k+l\leq M\}$, and the corresponding partial differential operator $\mathcal{L}_{M,\alpha}$. For any $(I,J)\in\mathbb{N}^{2}$ such that $I+J<q$, the quantity $N_{I,J}$ only depends on unknowns $\lambda_{i,j}$ with length of the multi-index at most equal to $M+I+J-1$.

Proof.

The result is straightforward from Lemmas 1 and 2. ∎

Assuming that all unknowns $\lambda_{i,j}$ with length of the multi-index at most equal to $M+I+J-1$ are known, then (21) is a well-defined linear under-determined system with

• •

$\mathfrak{L}$ linear equations, namely the $(I,J)=(I,\mathfrak{L}-I)$ equations from System (10) for $I$ between [math] and $\mathfrak{L}$;

• •

$M+\mathfrak{L}+1$ unknowns, namely the $\lambda_{ij}$ for $i+j=M+\mathfrak{L}$.

Therefore, if all unknowns $\lambda_{i,j}$ with length of the multi-index at most equal to $M+I+J-1$ are known, we expect to be able to compute a solution to the subsystem $\mathfrak{L}$ ; this is the layer structure of our original problem (10). Figure 2 highlights the link between the layers of unknowns and equations of the initial nonlinear system on the one hand, and the layers unknowns and equations of the linear subsystems on the other hand.

At this stage, we have identified a hierarchy of under-determined linear subsystems, for increasing values of $\mathfrak{L}$ from [math] to $q-1$, and we are now going to propose one procedure to build a solution to each subsystem. There is no unique way to do so, however if either $\alpha_{M,0}(x_{0},y_{0})\neq 0$ or $\alpha_{0,M}(x_{0},y_{0})\neq 0$ it provides a natural way to proceed. Indeed, the unknowns involved in an equation $(I,J)=(I,\mathfrak{L}-I)$ are $\{\lambda_{i,M+\mathfrak{L}-i};i\in\mathbb{N},I\leq i\leq I+M\}$ ; and the coefficient of the unknown $\lambda_{I+M,\mathfrak{L}-I}$ is proportional to $\alpha_{M,0}(x_{0},y_{0})$, which is non-zero under Hypothesis 1. Figure 3 provides two examples, in the $(i,j)$ plane, of the indices of one equation’s unknowns: for each equation, the coefficient of the term corresponding to the rightmost marker is non-zero.

By adding $M$ constraints corresponding to fixing the values of $\lambda_{i,M+\mathfrak{L}-i}$ for $0\leq i<M$, that is the unknowns corresponding in the $(i,j)$ plane to first $M$ markers on the left at level $M+\mathfrak{L}$, we therefore guarantee that for increasing values of $I$ from [math] to $\mathfrak{L}$ we can compute successively $\lambda_{I+M,\mathfrak{L}-I}$.

We can easily recast this in terms of matrices. At each level $\mathfrak{L}$, numbering the equations with increasing values of $I$ and the unknowns with increasing values of $i$ highlights the band-limited structure of each subsystem, while the entries of the $M$th super diagonal are all proportional to $\alpha_{M,0}(x_{0},y_{0})$, and therefore non-zero under Hypothesis 1. The matrix of the square linear system at level $\mathfrak{L}$ is then constructed from the first $M$ lines of the identity, corresponding to the additional $M$ constraints, placed on top of the matrix of the subsystem.

Definition 6**.**

Consider a point $(x_{0},y_{0})\in\mathbb{R}^{2}$, a given $q\in\mathbb{N}^{*}$, a given $M\in\mathbb{N}$, $M\geq 2$, a given set of complex-valued functions $\alpha=\{\alpha_{k,l}\in\mathcal{C}^{q-1}\text{ at }(x_{0},y_{0}),0\leq k+l\leq M\}$, and the corresponding partial differential operator $\mathcal{L}_{M,\alpha}$. For a given level $\mathfrak{L}\in\mathbb{N}$ with $\mathfrak{L}<q$, we define the matrix of the square system of level $\mathfrak{L}$, $\mathsf{T}^{\mathfrak{L}}\in\mathbb{C}^{(M+\mathfrak{L}+1)\times(M+\mathfrak{L}+1)}$, as

[TABLE]

or equivalently

[TABLE]

Assuming that all unknowns $\lambda_{i,j}$ with length of the multi-index at most equal to $M+I+J-1$ are known, then, as expected, a solution to the linear under-determined system (21) can be computed as follows.

Proposition 1**.**

Consider a point $(x_{0},y_{0})\in\mathbb{R}^{2}$, a given $q\in\mathbb{N}^{*}$, a given $M\in\mathbb{N}$, $M\geq 2$, a given set of complex-valued functions $\alpha=\{\alpha_{k,l}\in\mathcal{C}^{q-1}\text{ at }(x_{0},y_{0}),0\leq k+l\leq M\}$, and the corresponding partial differential operator $\mathcal{L}_{M,\alpha}$. For a given level $\mathfrak{L}\in\mathbb{N}$ with $\mathfrak{L}<q$, under Hypothesis 1, the matrix $\mathsf{T}^{\mathfrak{L}}\in\mathbb{C}^{(M+\mathfrak{L}+1)\times(M+\mathfrak{L}+1)}$ is non-singular.

We now assume that the unknowns $\{\lambda_{i,j},(i,j)\in\mathbb{N}^{2},i+j<M+\mathfrak{L}\}$ are known, so that the terms $N_{I,\mathfrak{L}-I}$ for $I$ from [math] to $\mathfrak{L}$ can be computed. Consider any vector $\mathsf{B}^{\mathfrak{L}}\in\mathbb{C}^{M+\mathfrak{L}+1}$ satisfying

[TABLE]

Then independently of the first $M$ components of $\mathsf{B}^{\mathfrak{L}}$, solving the linear system

[TABLE]

by forward substitution provides a solution to (21) for

[TABLE]

Proof.

The matrix $\mathsf{T}^{\mathfrak{L}}$ is lower triangular, therefore its determinant is

[TABLE]

which can not be zero under Hypothesis 1. The second part of the claim derives directly from the definition of $\mathsf{T}^{\mathfrak{L}}$ and $\mathsf{B}^{\mathfrak{L}}$ and the fact that the system is lower triangular, and can be illustrated as follows:

[TABLE]

∎

To summarize, we have defined for increasing values of $\mathfrak{L}$ a hierarchy of linear systems, each of which has the following characteristics:

• •

its unknowns are $\{\lambda_{i,M+\mathfrak{L}-i};\ \forall i\in\mathbb{N}\text{ such that }0\leq i\leq M+\mathfrak{L}\}$;

• •

its matrix $\mathsf{T}^{\mathfrak{L}}\in\mathbb{C}^{(M+\mathfrak{L}+1)\times(M+\mathfrak{L}+1)}$ is a square, non-singular, and triangular ;

• •

its right-hand side depends both on $\{\lambda_{i,j};\ \forall(i,j)\in\mathbb{N}^{2}\text{ such that }0\leq i+j<M+\mathfrak{L}\}$ and on $M$ additional parameters.

At each level $\mathfrak{L}$, assuming that the unknowns of inferior levels are known and provided $M$ given values for $\lambda_{i,M+\mathfrak{L}-i}$ for $0\leq i<M$, Proposition 1 provides an explicit formula to compute $\lambda_{i,M+\mathfrak{L}-i}$ for $M\leq i\leq M+\mathfrak{L}$.

2.4 Algorithm

The non-linear system (10) had $N_{dof}^{\eqref{thesyst}}=\frac{(M+q)(M+q+1)}{2}$ unknowns and $N_{eqn}^{\eqref{thesyst}}=\frac{q(q+1)}{2}$ equations, whereas each linear triangular system introduced in the previous subsection has $N_{dof}^{T}=M+\mathfrak{L}+1$ unknowns and $N_{eqn}^{T}=M+\mathfrak{L}+1$ equations for each level $\mathfrak{L}$ such that $0\leq\mathfrak{L}\leq q-1$. Therefore the hierarchy of triangular systems has a total of $N_{dof}^{H}=(M+1)q+\frac{q(q-1)}{2}$ unknowns and $N_{eqn}^{H}=N_{eqn}^{\eqref{thesyst}}+Mq=Mq+\frac{q(q+1)}{2}$ equations, including the $\frac{q(q+1)}{2}$ equations of the initial non-linear system (10).

The remaining $N_{dof}^{\eqref{thesyst}}-N_{dof}^{T}=\frac{M(M+1)}{2}$ unknowns, which are unknowns of none of the triangular systems but appear only on the right hand side of these systems, are the $\{\lambda_{i,j},(i,j)\in\mathbb{N}^{2},0\leq i+j<M\}$. These are the unknowns with length of the multi-index at most equal to $M-1$, and the corresponding indices $(i,j)$ are the only ones that are not marked on the right panel of Figure 2. It is therefore natural to add $\frac{M(M+1)}{2}$ constraints corresponding to fixing the values of the remaining unknowns $\{\lambda_{i,j},(i,j)\in\mathbb{N}^{2},0\leq i+j<M\}$. The final system we consider consists of these $\frac{M(M+1)}{2}$ constraints, guaranteeing that the unknowns $\{\lambda_{i,j},(i,j)\in\mathbb{N}^{2},0\leq i+j<M\}$ are known, together with the hierarchy of triangular systems (22) for increasing values of $\mathfrak{L}$ from [math] to $q-1$; it has $N_{dof}^{F}=\frac{(M+q)(M+q+1)}{2}$ unknowns, namely the unknowns of the original system (10), and $N_{eqn}^{F}=\frac{(M+q)(M+q+1)}{2}$ equations, namely the equations of the original system split into linear subsytems together with a total of $\frac{M(M+1)}{2}+qM$ additional constraints. A counting summary is presented here:

[TABLE]

Thanks to the $\frac{M(M+1)}{2}$ constraints, for increasing values of $\mathfrak{L}$ from [math] to $q-1$, the hypothesis of Proposition 1 is satisfied, the right hand side $\mathsf{B}^{\mathfrak{L}}$ can be evaluated and the triangular system (22) can be solved. So the unknowns $\{\lambda_{i,M+\mathfrak{L}-i};\ \forall i\in\mathbb{N}\text{ such that }0\leq i\leq M+\mathfrak{L}\}$ can be computed by induction on $\mathfrak{L}$, constructing a solution to the initial non-linear system (10) by induction on $\mathfrak{L}$.

The following algorithm requires the value of $\frac{M(M+1)}{2}+qM$ parameters, to fix initially the set of unknowns $\{\lambda_{i,j},(i,j)\in\mathbb{N}^{2},0\leq i+j<M\}$ and then at each level $\mathfrak{L}$ the set of unknowns $\{\lambda_{i,M+\mathfrak{L}-i},i\in\mathbb{N},0\leq i<M\}$. Under Hypothesis 1, the algorithm presents a sequence of steps to construct explicitly a solution to Problem (10) and requires no approximation process.

From the definitions of $\mathsf{T}^{\mathfrak{L}}$ and $\mathsf{B}^{\mathfrak{L}}$ we immediately see that the step $5$ boils down to

[TABLE]

If the set of unknowns $\{\lambda_{i,j},(i,j)\in\mathbb{N}^{2},0\leq i+j<M+q-1\}$ is computed from Algorithm 1, then the polynomial $\displaystyle P(x,y):=\sum_{0\leq i+j\leq q+M-1}\lambda_{i,j}(x-x_{0})^{i}(y-y_{0})^{j}$ is a solution to Problem (10), and therefore the function $\displaystyle\varphi(\mathbf{x}):=\exp P(\mathbf{x})$ satisfies (6). This is true independently of the values fixed in lines 1.1 and 1.3 of the algorithm.

Remark 1**.**

It is interesting to notice that the algorithm applies to a wide range of partial differential operators, including type changing operators such as Keldysh operators, $L_{K}=\partial_{x}^{2}+y^{2m+1}\partial_{y}^{2}+$ lower order terms, or Tricomi operators, $L_{T}=\partial_{x}^{2}+x^{2m+1}\partial_{y}^{2}+$ lower order terms, that change from elliptic to hyperbolic type along a smooth parabolic curve.

To conclude this section, we provide a formal definition of a GPW associated to an partial differential operator at a given point.

Definition 7**.**

Consider a point $(x_{0},y_{0})\in\mathbb{R}^{2}$, a given $q\in\mathbb{N}^{*}$, a given $M\in\mathbb{N}$, $M\geq 2$, a given set of complex-valued functions $\alpha=\{\alpha_{k,l}\in\mathcal{C}^{q-1}\text{ at }(x_{0},y_{0}),0\leq k+l\leq M\}$, and the corresponding partial differential operator $\mathcal{L}_{M,\alpha}$. A Generalized Plane Wave (GPW) associated to the differential operator $\mathcal{L}_{M,\alpha}$ at the point $(x_{0},y_{0})$ is a function $\varphi$ satisfying

[TABLE]

Under Hypothesis 1, a Generalized Plane Wave (GPW) can be constructed as function $\varphi(x,y)=\exp P(x,y)$, where the coefficients of the polynomial $P$ are computed by Algorithm 1, independently of the values fixed in the algorithm.

The crucial feature of the construction process is the exact solution provided in the algorithm: in practice, a solution to the initial non-linear rectangular system is computed without numerical resolution of any system, with an explicit formula.

The choice of the fixed values in Algorithm 1 will be discussed in the next paragraph. Even though these values does not affect the construction process, and the fact that the corresponding $\varphi(x,y)=\exp P(x,y)$ is a GPW, it will be key to prove the interpolation properties of the corresponding set of GPWs.

Remark 2**.**

Under the hypothesis $\alpha_{0,M}(x_{0},y_{0})\neq 0$ it would be natural to fix the values of $\{\lambda_{i,j},0\leq j\leq M-1,0\leq i\leq q+M-1-j\}$ instead of those of $\{\lambda_{i,j},0\leq i\leq M-1,0\leq j\leq q+M-1-i\}$, and an algorithm very similar to Algorithm 1, exchanging the roles of $i$ and $j$ would construct the polynomial coefficients of a GPW.

3 Normalization

We will refer to normalization as the choice of imposed values in Algorithm 1. The discussion presented in this section will be summarized in Definition 8.

Within the construction process presented in the previous section, only the design of the function $\varphi$ as the exponential of a polynomial is related to wave propagation, while Algorithm 1 works for partial differential operators not necessarily related to wave propagation. In particular, the property $\mathcal{L}_{M,\alpha}\varphi(x,y)=O\left(\|(x,y)-(x_{0},y_{0})\|^{q}\right)$ of GPWs is independent of the choice of $(\lambda_{1,0},\lambda_{0,1})$. However, the normalization process described here carries on the idea of adding higher order terms to the phase function of a plane wave, see (3), as was proposed in [13].

We will now restrict our attention to a smaller set of partial differential operators that include several interesting operators related to wave propagation, thanks to an additional hypothesis on the highest order derivatives in $\mathcal{L}_{M,\alpha}$, namely Hypothesis 2. Under this hypothesis we will be able to study the interpolation properties of associated GPWs in a unified framework. As we will see in this section, choosing only two non-zero fixed values in Algorithm 1 is sufficient to generate a set of linearly independent GPWs. It is then natural to study how the rest of the $\lambda_{ij}$s depend on those two values, and the related consequences of Hypothesis 2. These rely on Hypothesis 2 extending the fact that for classical PWs $(i\kappa\cos\theta)^{2}+(i\kappa\sin\theta)^{2}=-\kappa^{2}$ is independent of $\theta$.

3.1 For every GPWs

In Algorithm 1, the number of prescribed coefficients is $\frac{M(M+1)}{2}+Mq$, and the set of coefficients to be prescribed is the set $\{\lambda_{i,j},0\leq i\leq M-1,0\leq j\leq q+M-1-i\}$.

For the sake of simplicity, it is natural to choose most of these values to be zero. Since the unknown $\lambda_{0,0}$ never appears in the non-linear system, there is nothing more natural than setting it to zero: this ensures that any GPW $\varphi$ will satisfy $\varphi(x_{0},y_{0})=1$. Concerning the subset of $Mq$ unknowns corresponding to step 1.3 in Algorithm 1, setting these values to zero simply reduces the amount of computation involved in step 1.5 in the algorithm: indeed for $I=0$ then $\displaystyle\sum_{k=0}^{M-1}\mathsf{T}^{\mathfrak{L}}_{I+M+1,I+k+1}\lambda_{I+k,M+\mathfrak{L}-I-k}=0$, while for $0<I<M$ then

[TABLE]

As for the unknowns $\lambda_{1,0}$ and $\lambda_{0,1}$, they will be non-zero to mimic the classical plane wave case, and their precise choice will be discussed in the next subsection. For the remaining unknowns to be fixed, that is to say the set $\{\lambda_{i,j},2\leq i+j\leq M-1\}$, their values are set to zero, here again in order to reduce the amount of computation in computing the right hand side entries $\mathrm{B}^{\mathfrak{L}}_{M+1+I}$ and in applying 1.5.

For the operator $\mathfrak{L}_{2,\gamma}$ the non-linear terms in $N_{1,0}$ and $N_{0,1}$ respectively become with this normalization

[TABLE]

[TABLE]

Since all but two of the unknowns to be fixed in Algorithm 1 are set to zero, it is now natural to express the $\lambda_{i,j}$ unknowns computed from 1.5 in the algorithm as functions of the two non-zero prescribed unknowns, $\lambda_{1,0}$ and $\lambda_{0,1}$.

Lemma 3**.**

Consider a point $(x_{0},y_{0})\in\mathbb{R}^{2}$, a given $q\in\mathbb{N}^{*}$, a given $M\in\mathbb{N}$, $M\geq 2$, a given set of complex-valued functions $\alpha=\{\alpha_{k,l}\in\mathcal{C}^{q-1}\text{ at }(x_{0},y_{0}),0\leq k+l\leq M\}$, and the corresponding partial differential operator $\mathcal{L}_{M,\alpha}$. Under Hypothesis 1 consider a solution to Problem (10) constructed thanks to Algorithm 1 with all the prescribed values $\lambda_{i,j}$ such that $i<M$ and $i+j\neq 1$ set to zero. Each $\lambda_{i+M,j}$ can be expressed as an element of $\mathbb{C}[\lambda_{1,0},\lambda_{0,1}]$.

Proof.

The fact that $\lambda_{i+M,j}$ can be expressed as a polynomial in two variables with respect to $\lambda_{1,0}$ and $\lambda_{0,1}$ is a direct consequence from the explicit formula in step 1.5 in Algorithm 1 combining with setting $\lambda_{i,j}$ such that $i<M$ and $i+j\neq 1$ to zero. ∎

Since unknowns are expressed as elements of $\mathbb{C}[\lambda_{1,0},\lambda_{0,1}]$, we will now study the degree of various terms from Algorithm 1 as polynomials with respect to $\lambda_{1,0}$ and $\lambda_{0,1}$. To do so, we will start by inspecting the product terms appearing in Faa di Bruno’s formula.

Lemma 4**.**

Consider a point $(x_{0},y_{0})\in\mathbb{R}^{2}$, a given $q\in\mathbb{N}^{*}$, a given $M\in\mathbb{N}$, $M\geq 2$, a given set of complex-valued functions $\alpha=\{\alpha_{k,l}\in\mathcal{C}^{q-1}\text{ at }(x_{0},y_{0}),0\leq k+l\leq M\}$, and the corresponding partial differential operator $\mathcal{L}_{M,\alpha}$. Consider a given polynomial $P\in\mathbb{C}[x,y]$. The non-linear terms $\mathcal{L}_{M,\alpha}^{N}P$, expressed as linear combinations of products of derivatives of $P$, namely $\prod_{m=1}^{s}\left(\partial_{x}^{i_{m}}\partial_{y}^{j_{m}}P\right)^{k_{m}}$, contain products of up to $M$ derivatives of $P$, namely $\partial_{x}^{i_{m}}\partial_{y}^{j_{m}}P$, counting repetitions. The only products that have exactly $M$ terms are $(\partial_{x}P)^{k}(\partial_{y}P)^{M-k}$ for $0\leq k\leq M$, whereas all the other products have less than $M$ terms.

Proof.

Since the operator $\mathcal{L}_{M,\alpha}^{N}$ is defined via Faa di Bruno’s formula, we will proceed by careful examination of the summation and product indices in the latter.

The number of terms in the product term is $s$, with possible repetitions counted thanks to the $k_{m}$s, and the total number of terms counting repetitions is $\mu=\sum_{m=1}^{s}k_{m}$. Since in $\mathcal{L}_{M,\alpha}^{N}$ the indices are such that $1\leq\mu\leq\ell\leq M$, there cannot be more than $M$ terms counting repetitions in any of the $\prod_{m=1}^{s}\left(\partial_{x}^{i_{m}}\partial_{y}^{j_{m}}P\right)^{k_{m}}$.

For $s=1$, in the set $p_{1}((k,\ell-k),\mu)$, $(i_{1},j_{1})\in\mathbb{N}^{2}$ are such that $i_{1}+j_{1}\geq 1$ and $k_{1}\in\mathbb{N}$ is such that $k_{1}(i_{1}+j_{1})=\ell$. Since $\ell\leq M$, such a term appears in Faa di Bruno’s formula as a product of $\mu=M$ terms if and only if $\ell=M$, $k_{1}=M$, and therefore $i_{1}+j_{1}=1$. There are then only two possibilities: either $(i_{1},j_{1})=(1,0)$ corresponding to the term $(\partial_{x}P)^{M}$, or $(i_{1},j_{1})=(0,1)$ corresponding to the term $(\partial_{y}P)^{M}$.

For $s=2$, in the set $p_{2}((k,\ell-k),\mu)$, $(i_{1},j_{1},i_{2},j_{2})\in\mathbb{N}^{4}$ are such that $i_{1}+j_{1}\geq 1$, $i_{2}+j_{2}\geq 1$, $(i_{1},j_{1})\prec(i_{2},j_{2})$, and $(k_{1},k_{2})\in\mathbb{N}^{2}$ is such that $\mu=k_{1}+k_{2}$ and $k_{1}(i_{1}+j_{1})+k_{2}(i_{2}+j_{2})=\ell$. Since $\ell\leq M$ and $\ell=k_{1}(i_{1}+j_{1})+k_{2}(i_{2}+j_{2})\geq k_{1}+k_{2}=\mu$, such a term appears in Faa di Bruno’s formula as a product of $\mu=M$ terms if and only if $\ell=M$ and $k_{1}+k_{2}=M$. There are then two possible cases: either $i_{2}+j_{2}>1$, then $M=k_{1}(i_{1}+j_{1})+k_{2}(i_{2}+j_{2})>k_{1}+k_{2}=M$, so there is no such term in the sum, or $i_{2}+j_{2}=1$, then necessarily $(i_{1},j_{1})=(0,1)$ and $(i_{2},j_{2})=(1,0)$, corresponding to the terms $(\partial_{x}P)^{k}(\partial_{y}P)^{M-k}$ for any $k$ from [math] to $M$.

For $s\geq 3$, in the set $p_{s}((k,\ell-k),\mu)$, for all $m\in\mathbb{N}$ such that $1\leq m\leq s$, $(i_{m},j_{m})\in\mathbb{N}^{2}$ and $k_{m}\in\mathbb{N}$ are such that $i_{m}+j_{m}\geq 1$, $\sum_{m=1}^{s}k_{m}(i_{m}+j_{m})=\ell$, $\mu=\sum_{m=1}^{s}k_{m}$ and $(i_{1},j_{1})\prec(i_{2},j_{2})\prec(i_{3},j_{3})$. Because of this last condition, it is clear that $i_{3}+j_{3}>1$. Since $\ell\leq M$ and $\ell=\sum_{m=1}^{s}k_{m}(i_{m}+j_{m})\geq\sum_{m=1}^{s}k_{m}=\mu$, such a term appears in Faa di Bruno’s formula as a product of $\mu=M$ terms if and only if $\ell=M$ and $\sum_{m=1}^{s}k_{m}=M$. But then $M=\sum_{m=1}^{s}k_{m}(i_{m}+j_{m})>\sum_{m=1}^{s}k_{m}=M$, so there is no such term in the sum.

The claim is proved. ∎

Lemma 5**.**

Consider a point $(x_{0},y_{0})\in\mathbb{R}^{2}$, a given $q\in\mathbb{N}^{*}$, a given $M\in\mathbb{N}$, $M\geq 2$, a given set of complex-valued functions $\alpha=\{\alpha_{k,l}\in\mathcal{C}^{q-1}\text{ at }(x_{0},y_{0}),0\leq k+l\leq M\}$, and the corresponding partial differential operator $\mathcal{L}_{M,\alpha}$. Consider a given polynomial $P\in\mathbb{C}[x,y]$. The quantity $\partial_{x}^{I_{0}}\partial_{y}^{J_{0}}\mathcal{L}_{M,\alpha}^{N}P$ is a linear combination of terms $\partial_{x}^{I_{0}}\partial_{y}^{J_{0}}\left(\prod_{m=1}^{s}\left(\partial_{x}^{i_{m}}\partial_{y}^{j_{m}}P\right)^{k_{m}}\right)$, where the indices come from Faa di Bruno’s formula. Each of these $\partial_{x}^{I_{0}}\partial_{y}^{J_{0}}\left(\prod_{m=1}^{s}\left(\partial_{x}^{i_{m}}\partial_{y}^{j_{m}}P\right)^{k_{m}}\right)$ can be expressed as a linear combination of products $\prod_{m=1}^{t}(\partial_{x}^{a_{m}}\partial_{y}^{b_{m}}P)^{c_{m}}$ where the indices satisfy $\sum_{m=1}^{t}c_{m}(a_{m}+b_{m})\leq{I_{0}}+{J_{0}}+M$.

Proof.

Thanks to the product rule, the derivative $\partial_{x}^{I_{0}}\partial_{y}^{J_{0}}\left(\prod_{m=1}^{s}\left(\partial_{x}^{i_{m}}\partial_{y}^{j_{m}}P\right)^{k_{m}}\right)$ can be expressed as a linear combination of several terms $\prod_{m=1}^{s}\partial_{x}^{I_{m}}\partial_{y}^{J_{m}}\left[\left(\partial_{x}^{i_{m}}\partial_{y}^{j_{m}}P\right)^{k_{m}}\right]$, where $\sum_{m=1}^{t}I_{m}=I_{0}$ and $\sum_{m=1}^{t}J_{m}=J_{0}$.

We can prove by induction on $k$ that $\partial_{x}^{I}\partial_{y}^{J}\left[\left(\partial_{x}^{i}\partial_{y}^{j}P\right)^{k}\right]$ can be expressed, for all $(i,j,I,J)\in\mathbb{N}^{4}$, as a linear combination of products $\prod_{m=1}^{M}(\partial_{x}^{a_{m}}\partial_{y}^{b_{m}}P)^{c_{m}}$ where the indices satisfy $\sum_{m=1}^{M}c_{m}(a_{m}+b_{m})\leq I+J+k(i+j)$:

it is evidently true for $k=1$; 2. 2.

suppose that it is true for $k_{0}\geq 1$, then for any $(i,j,I,J)\in\mathbb{N}^{4}$ the product rule applied to $\partial_{x}^{i}\partial_{y}^{j}P\times\left(\partial_{x}^{i}\partial_{y}^{j}P\right)^{k_{0}}$ yields

[TABLE]

where by hypothesis each $\partial_{x}^{\tilde{i}}\partial_{y}^{\tilde{j}}\left[\left(\partial_{x}^{i}\partial_{y}^{j}P\right)^{k_{0}}\right]$can be expressed as a linear combination of products $\prod_{m=1}^{M}(\partial_{x}^{a_{m}}\partial_{y}^{b_{m}}P)^{c_{m}}$ with $\sum_{m=1}^{M}c_{m}(a_{m}+b_{m})\leq\tilde{i}+\tilde{j}+k_{0}(i+j)$, so that each term in the double sum can be expressed as a linear combination of products $\prod_{m=1}^{M+1}(\partial_{x}^{a_{m}}\partial_{y}^{b_{m}}P)^{c_{m}}$ where $a_{M+1}:=i+I-\tilde{i}$, $b_{M+1}:=j+J-\tilde{j}$ and $c_{M+1}:=1$, which yields $\sum_{m=1}^{M+1}c_{m}(a_{m}+b_{m})=\sum_{m=1}^{M}c_{m}(a_{m}+b_{m})+(i+I-\tilde{i}+j+J-\tilde{j})$ and therefore $\sum_{m=1}^{M+1}c_{m}(a_{m}+b_{m})\leq k_{0}(i+j)+(i+I+j+J)$. This concludes the proof by induction.

Finally the derivative $\partial_{x}^{I_{0}}\partial_{y}^{J_{0}}\left(\prod_{m=1}^{s}\left(\partial_{x}^{i_{m}}\partial_{y}^{j_{m}}P\right)^{k_{m}}\right)$ can be expressed as a linear combination of several terms $\prod_{m=1}^{s}\prod_{\tilde{m}=1}^{M}(\partial_{x}^{a_{\tilde{m}}}\partial_{y}^{b_{\tilde{m}}}P)^{c_{\tilde{m}}}$, with $\sum_{\tilde{m}=1}^{M}c_{\tilde{m}}(a_{\tilde{m}}+b_{\tilde{m}})\leq I_{m}+J_{m}+k_{m}(i_{m}+j_{m})$, in other words it can be expressed as a linear combination of several terms $\prod_{m=1}^{Ms}(\partial_{x}^{a_{m}}\partial_{y}^{b_{m}}P)^{c_{m}}$, with $\sum_{m=1}^{Ms}c_{m}(a_{m}+b_{m})\leq\sum_{m=1}^{s}I_{m}+J_{m}+k_{m}(i_{m}+j_{m})=I_{0}+J_{0}+\sum_{m=1}^{s}k_{m}(i_{m}+j_{m})$. For any $\partial_{x}^{I_{0}}\partial_{y}^{J_{0}}\left(\prod_{m=1}^{s}\left(\partial_{x}^{i_{m}}\partial_{y}^{j_{m}}P\right)^{k_{m}}\right)$ coming from $\partial_{x}^{I_{0}}\partial_{y}^{J_{0}}\mathcal{L}_{M,\alpha}^{N}P$, the summation indices from Faa di Bruno’s formula satisfy $\sum_{m=1}^{s}k_{m}(i_{m}+j_{m})=\ell$, so the products $\prod_{m=1}^{Ms}(\partial_{x}^{a_{m}}\partial_{y}^{b_{m}}P)^{c_{m}}$ are such that $\sum_{m=1}^{Ms}c_{m}(a_{m}+b_{m})\leq I_{0}+J_{0}+M$. ∎

The two following results now turn to $\lambda_{i+M,j}$ computed in Algorithm 1.

Proposition 2**.**

Consider a point $(x_{0},y_{0})\in\mathbb{R}^{2}$, a given $q\in\mathbb{N}^{*}$, a given $M\in\mathbb{N}$, $M\geq 2$, a given set of complex-valued functions $\alpha=\{\alpha_{k,l}\in\mathcal{C}^{q-1}\text{ at }(x_{0},y_{0}),0\leq k+l\leq M\}$, and the corresponding partial differential operator $\mathcal{L}_{M,\alpha}$. Under Hypothesis 1 consider a solution to Problem (10) constructed thanks to Algorithm 1 with all the fixed values $\lambda_{i,j}$ such that $i<M$ and $i+j\neq 1$ set to zero. As an element of $\mathbb{C}[\lambda_{1,0},\lambda_{0,1}]$, $\lambda_{M,0}$ is of degree equal to $M$.

Proof.

The formula to compute $\lambda_{M,0}$ in Algorithm 1 comes from the $(I,J)=(0,0)$ equation in System (10), that is to say $\mathcal{L}_{M,\alpha}^{A}P(x_{0},y_{0})=-\alpha_{0,0}(x_{0},y_{0})$. It reads

[TABLE]

and the sum is actually zero since the $\lambda_{k,M-k}$ unknowns are prescribed to zero for $k<M$. The definitions of $\mathsf{B}^{0}$ and $\mathsf{L}^{0}$ then give

[TABLE]

Since the $\lambda_{k,\ell-k}$ unknowns are prescribed to zero for all $1<\ell<M-1$ and all $k$, the double sum term reduces to $\alpha_{0,1}(x_{0},y_{0})\lambda_{0,1}+\alpha_{1,0}(x_{0},y_{0})\lambda_{1,0}$. The non-linear terms from $\mathcal{L}_{M,\alpha}^{N}P$, namely $\prod_{m=1}^{s}(\partial_{x}^{i_{m}}\partial_{y}^{j_{m}}P)^{k_{m}}$, are products of at most $M$ terms, counting repetitions, according to Lemma 4. So $\mathcal{L}_{M,\alpha}^{N}P(x_{0},y_{0})$ is a linear combination of product terms reading $\prod_{m=1}^{s}(\lambda_{i_{m},j_{m}})^{k_{m}}$ with at most $M$ factors. Moreover, since $P$ is constructed thanks to Algorithm 1, from Corollary 1 we know that these $\lambda_{i_{m},j_{m}}$s have a length of the multi-index at most equal to $M-1$, so they are either $\lambda_{1,0}$ or $\lambda_{0,1}$ or prescribed to zero. This means that in $\mathbb{C}[\lambda_{1,0},\lambda_{0,1}]$ each one of these $\lambda_{i_{m},j_{m}}$ is at most of degree one. So in $\mathbb{C}[\lambda_{1,0},\lambda_{0,1}]$ each $\prod_{m=1}^{s}(\lambda_{i_{m},j_{m}})^{k_{m}}$ is a product of at most $M$ factors each of them of degree at most one, the product is therefore of degree at most $M$. As a result

[TABLE]

as an element of $\mathbb{C}[\lambda_{1,0},\lambda_{0,1}]$ is of degree at most $M$.

Finally, the term $(\partial_{x}P)^{M}$ from $\mathcal{L}_{M,\alpha}^{N}P$ identified in Lemma 4 corresponds to a term $\alpha_{M,0}(x_{0},y_{0})(\lambda_{1,0})^{M}$ in the expression of $\lambda_{M,0}$, and this term is non-zero under Hypothesis 1. As a conclusion $\lambda_{M,0}$ as an element of $\mathbb{C}[\lambda_{1,0},\lambda_{0,1}]$ is of degree equal to $M$. ∎

Proposition 3**.**

Consider a point $(x_{0},y_{0})\in\mathbb{R}^{2}$, a given $q\in\mathbb{N}^{*}$, a given $M\in\mathbb{N}$, $M\geq 2$, a given set of complex-valued functions $\alpha=\{\alpha_{k,l}\in\mathcal{C}^{q-1}\text{ at }(x_{0},y_{0}),0\leq k+l\leq M\}$, and the corresponding partial differential operator $\mathcal{L}_{M,\alpha}$. Under Hypothesis 1 consider a solution to Problem (10) constructed thanks to Algorithm 1 with all the fixed values $\lambda_{i,j}$ such that $i<M$ and $i+j\neq 1$ set to zero. As an element of $\mathbb{C}[\lambda_{1,0},\lambda_{0,1}]$, each $\lambda_{i+M,j}$ has a total degree at most equal to the length of its multi-index $i+j+M$.

Proof.

The formula to compute $\lambda_{I+M,\mathfrak{L}-I}$ in Algorithm 1 comes from the $(I,J)=(I,\mathfrak{L}-I)$ equation in System (10), that is to say $\partial_{x}^{I}\partial_{y}^{\mathfrak{L}-I}\mathcal{L}_{M,\alpha}^{A}P(x_{0},y_{0})=-\partial_{x}^{I}\partial_{y}^{\mathfrak{L}-I}\alpha_{0,0}(x_{0},y_{0})$. It reads

[TABLE]

We will proceed by induction on $\mathfrak{L}$:

the result has been proved to be true for $\mathfrak{L}=0$ in Proposition 2 ; 2. 2.

suppose the result is true for $\mathfrak{L}\in\mathbb{N}$ as well as for all $\tilde{\mathfrak{L}}\in\mathbb{N}$ such that $\tilde{\mathfrak{L}}\leq\mathfrak{L}$, then all the linear terms in $N_{I,\mathfrak{L}+1-I}$ have a length of the multi-index at most equal to $M+\mathfrak{L}$ so by hypothesis their degree as elements of $\mathbb{C}[\lambda_{1,0},\lambda_{0,1}]$ is at most equal to $M+\mathfrak{L}$, and thanks to Lemma 5 all the non-linear terms in $N_{I,\mathfrak{L}+1-I}$ can be expressed as a linear combination of products $\prod_{m=1}^{t}(\lambda_{a_{m},b_{m}})^{c_{m}}$ where the indices satisfy $\sum_{m=1}^{t}c_{m}(a_{m}+b_{m})\leq\mathfrak{L}+1+M$ so by hypothesis their degree as elements of $\mathbb{C}[\lambda_{1,0},\lambda_{0,1}]$ is at most equal to $M+\mathfrak{L}+1$ ; the last step is to prove that the $\lambda_{I+k,M+\mathfrak{L}+1-I-k}$ are also of degree at most equal to $M+\mathfrak{L}+1$, and we will proceed by induction on $I$:

1. (a)

for $I=0$, all $\lambda_{I+k,M+\mathfrak{L}+1-I-k}$ for $0\leq k\leq M-1$ satisfy the two conditions $I+k<M$ and $I+k+M+\mathfrak{L}+1-I-k=M+\mathfrak{L}+1\neq 1$ so they are all prescribed to zero and their degree as element of $\mathbb{C}[\lambda_{1,0},\lambda_{0,1}]$ is at most equal to $M+\mathfrak{L}+1$ that ; 2. (b)

suppose that, for a given $I\in\mathbb{N}$, the $\lambda_{\tilde{I}+k,M+\mathfrak{L}+1-\tilde{I}-k}$ for all $\tilde{I}\in\mathbb{N}$ such that $\tilde{I}\leq I$ are also of degree at most equal to $M+\mathfrak{L}+1$ then it is clear from Equation (24) that $\lambda_{I+1+M,\mathfrak{L}-I-1}$ is also of degree at most equal to $M+\mathfrak{L}+1$.

This concludes the proof.

∎

As explained from an algebraic viewpoint in section 3.2 in [13], the degree of $\lambda_{i+M,j}$ as an element of $\mathbb{C}[\lambda_{1,0},\lambda_{0,1}]$ will be affected by the choice of the last two prescribed values, namely $\lambda_{1,0}$ and $\lambda_{0,1}$. Indeed if $\lambda_{1,0}$ and $\lambda_{0,1}$ satisfy a polynomial identity $P_{l}(\lambda_{1,0},\lambda_{0,1})=0$, then we can consider the quotient ring $\mathbb{C}[\lambda_{1,0},\lambda_{0,1}]/(P_{l})$.

Note that choosing to set $\{\lambda_{i,j},1<i+j\leq M-1\}$ to values different from zero may be useful to treat operators that do not satisfy Hypothesis 2 but this is not our goal here.

3.2 For each GPW

In order to obtain a set of linearly independent GPWs, the values of $\lambda_{1,0}$ and $\lambda_{0,1}$ will be chosen different for each GPW. However the values of $\lambda_{1,0}$ and $\lambda_{0,1}$ will satisfy a common property for every GPWs. Very much as the coefficients of any plane wave of wavenumber $\kappa$ satisfy $(\lambda_{1,0})^{2}+(\lambda_{0,1})^{2}=-\kappa^{2}$, independently of the direction of propagation $\theta$ since $\lambda_{1,0}=\imath\kappa\cos\theta$ and $\lambda_{0,1}=\imath\kappa\sin\theta$, under Hypothesis 2 the coefficients of each GPW will be chosen for the quantity

[TABLE]

to be identical for every GPWs, as we will see in the following proposition and theorem.

This will be crucial to prove interpolation properties of the corresponding set of functions, which will result from the consequence of this common property on the degree of each $\lambda_{i+M,j}$ as an element of $\mathbb{C}[\lambda_{1,0},\lambda_{0,1}]$. As the plane wave case suggests, we will see that $\lambda_{i+M,j}$ can be expressed as a polynomial of lower degree thanks to a judicious choice for $\lambda_{1,0}$ and $\lambda_{0,1}$.

We first need an intermediate result concerning the polynomial $\mathcal{L}_{M,\alpha}^{N}P$.

Lemma 6**.**

Consider a point $(x_{0},y_{0})\in\mathbb{R}^{2}$, a given $q\in\mathbb{N}^{*}$, a given $M\in\mathbb{N}$, $M\geq 2$, a given set of complex-valued functions $\alpha=\{\alpha_{k,l}\in\mathcal{C}^{q-1}\text{ at }(x_{0},y_{0}),0\leq k+l\leq M\}$, and the corresponding partial differential operator $\mathcal{L}_{M,\alpha}$. Consider a given polynomial $P\in\mathbb{C}[x,y]$. For any $\mathfrak{L}\in\mathbb{N}$ and any $I\in\mathbb{N}$ such that $I\leq\mathfrak{L}+1$, the quantity $\partial_{x}^{I}\partial_{y}^{\mathfrak{L}+1-I}\left[\mathcal{L}_{M,\alpha}^{N}P\right]$ can be expressed as a linear combination of products $\prod_{t=1}^{\mu}\partial_{x}^{i_{t}+I_{t}}\partial_{y}^{j_{t}+J_{t}}P$, with $\sum_{t=1}^{\mu}I_{t}=I$, $\sum_{t=1}^{\mu}J_{t}=\mathfrak{L}+1-I$, $\sum_{t=1}^{\mu}i_{t}=k$, and $\sum_{t=1}^{\mu}j_{t}=\ell-k$. Moreover, for each product term, there exists $t_{0}\in\mathbb{N}$, $1\leq t_{0}\leq\mu$ such that $I_{t_{0}}\neq 0$ or $J_{t_{0}}\neq 0$.

Proof.

The quantity $\mathcal{L}_{M,\alpha}^{N}P$ can be expressed, from Faa di Bruno’s formula, as a linear combination of products $\prod_{m=1}^{s}\left(\partial_{x}^{i_{m}}\partial_{y}^{j_{m}}P\right)^{k_{m}}$, with $(i_{m_{1}},j_{m_{1}})\neq(i_{m_{2}},j_{m_{2}})$ for all $m_{1}\neq m_{2}$, $\sum_{m=1}^{s}k_{m}=\mu$, $\sum_{m=1}^{s}k_{m}i_{m}=k$, and $\sum_{m=1}^{s}k_{m}j_{m}=\ell-k$. Therefore $\mathcal{L}_{M,\alpha}^{N}P$ can also be expressed, repeating terms, as a linear combination of products $\prod_{t=1}^{\mu}\partial_{x}^{i_{t}}\partial_{y}^{j_{t}}P$, with possibly $(i_{m_{1}},j_{m_{1}})=(i_{m_{2}},j_{m_{2}})$ for $m_{1}\neq m_{2}$, $\sum_{t=1}^{\mu}i_{t}=k$, and $\sum_{t=1}^{\mu}j_{t}=\ell-k$. So the quantity $\partial_{x}^{I}\partial_{y}^{\mathfrak{L}+1-I}\left[\mathcal{L}_{M,\alpha}^{N}P\right]$ can be expressed, from Leibniz’s rule, as a linear combination of products $\prod_{t=1}^{\mu}\partial_{x}^{i_{t}+I_{t}}\partial_{y}^{j_{t}+J_{t}}P$, with $\sum_{t=1}^{\mu}I_{t}=I$ and $\sum_{t=1}^{\mu}J_{t}=\mathfrak{L}+1-I$.

Consider such a given product term $\prod_{t=1}^{\mu}\partial_{x}^{i_{t}+I_{t}}\partial_{y}^{j_{t}+J_{t}}P$, and suppose that for all $t$ $I_{t}=J_{t}=0$. Then $I=\sum_{t=1}^{\mu}I_{t}=0$ and $\mathfrak{L}+1-I=\sum_{t=1}^{\mu}J_{t}=0$, which is impossible since $\mathfrak{L}+1>0$. ∎

The two following results gather the consequences of this choice on $\lambda_{i+M,j}$s computed in Algorithm 1.

Proposition 4**.**

Consider a point $(x_{0},y_{0})\in\mathbb{R}^{2}$, a given $q\in\mathbb{N}^{*}$, a given $M\in\mathbb{N}$, $M\geq 2$, a given set of complex-valued functions $\alpha=\{\alpha_{k,l}\in\mathcal{C}^{q-1}\text{ at }(x_{0},y_{0}),0\leq k+l\leq M\}$, and the corresponding partial differential operator $\mathcal{L}_{M,\alpha}$. Under Hypotheses 1 and 2 consider a solution to Problem (10) constructed thanks to Algorithm 1 with all the prescribed values $\lambda_{i,j}$ such that $i<M$ and $i+j\neq 1$ set to zero, and

[TABLE]

for some $\theta\in\mathbb{R}$ and $\kappa\in\mathbb{C}^{*}$. As an element of $\mathbb{C}[\lambda_{1,0},\lambda_{0,1}]$, $\lambda_{M,0}$ can be expressed as a polynomial of degree at most equal to $M-1$, and its coefficients are independent of $\theta$.

Note that once we impose this condition on $\lambda_{1,0},\lambda_{0,1}$ any element of $\mathbb{C}[\lambda_{1,0},\lambda_{0,1}]$ can be expressed by different polynomials, possibly with different degrees, simply because under Hypothesis 2 and (25) we have

[TABLE]

See paragraph 3.2 in [13] for an algebraic view point on this comment.

Proof.

Since

[TABLE]

again the term to investigate is $\mathcal{L}_{M,\alpha}^{N}P(x_{0},y_{0})$. Lemma 4 identifies products of $M$ terms in $\mathcal{L}_{M,\alpha}^{N}P$, and from the definition of $\mathcal{L}_{M,\alpha}^{N}$ they appear in the following linear combination

[TABLE]

Back to the expression of $\lambda_{M,0}$, and thanks to Hypothesis 2, the only possible terms of degree $M$ therefore appear in the following linear combination:

[TABLE]

Finally thanks to (25), the only terms of degree $M$ in (26) can be expressed as a polynomial of degree at most equal $M-1$. ∎

Proposition 5**.**

Consider a point $(x_{0},y_{0})\in\mathbb{R}^{2}$, a given $q\in\mathbb{N}^{*}$, a given $M\in\mathbb{N}$, $M\geq 2$, a given set of complex-valued functions $\alpha=\{\alpha_{k,l}\in\mathcal{C}^{q-1}\text{ at }(x_{0},y_{0}),0\leq k+l\leq M\}$, and the corresponding partial differential operator $\mathcal{L}_{M,\alpha}$. Under Hypotheses 1 and 2 consider a solution to Problem (10) constructed thanks to Algorithm 1 with all the fixed values $\lambda_{i,j}$ such that $i<M$ and $i+j\neq 1$ set to zero, and

[TABLE]

for some $\theta\in\mathbb{R}$ and $\kappa\in\mathbb{C}^{*}$. As an element of $\mathbb{C}[\lambda_{1,0},\lambda_{0,1}]$, each $\lambda_{i+M,j}$ can be expressed as a polynomial of degree at most equal to $i+j+M-1$, and its coefficients are independent of $\theta$.

Proof.

From Algorithm 1 the expression of $\lambda_{I+M,\mathfrak{L}-I}$ reads

[TABLE]

We will proceed again by induction on $\mathfrak{L}$:

the result has been proved to be true for $\mathfrak{L}=0$ in Proposition 4 ; 2. 2.

suppose the result is true for $\mathfrak{L}\in\mathbb{N}$ as well as for all $\tilde{\mathfrak{L}}\in\mathbb{N}$ such that $\tilde{\mathfrak{L}}\leq\mathfrak{L}$, then we focus on $N_{I,\mathfrak{L}+1-I}$, given by

[TABLE]

[TABLE]

all the linear terms in $N_{I,\mathfrak{L}+1-I}$, as elements of $\mathbb{C}[\lambda_{1,0},\lambda_{0,1}]$, by hypothesis have degree at most equal to $(I+M)+(\mathfrak{L}+1-I)-1=M+\mathfrak{L}$, and thanks to Lemma 6 all the non-linear terms in $N_{I,\mathfrak{L}+1-I}$ can be expressed as a linear combination of products $\prod_{t=1}^{\mu}\lambda_{a_{t},b_{t}}$ where the indices satisfy $\sum_{t=1}^{\mu}(a_{t}+b_{t})\leq\mathfrak{L}+1+M$ ; in each such product, as element of $\mathbb{C}[\lambda_{1,0},\lambda_{0,1}]$, each $\lambda_{a_{t},b_{t}}$ is either of degree $a_{t}+b_{t}=1$ if $(a_{t},b_{t})\in\{(0,1),(1,0)\}$, or of degree at most equal to $a_{t}+b_{t}-1$ otherwise by hypothesis ; from Lemma 6 there is at least one $t_{0}$ such that $(a_{t_{0}},b_{t_{0}})\notin\{(0,1),(1,0)\}$, therefore each product $\prod_{t=1}^{\mu}\lambda_{a_{t},b_{t}}$, as element of $\mathbb{C}[\lambda_{1,0},\lambda_{0,1}]$, can be expressed as a polynomial of degree at most $\left(\sum_{t=1}^{\mu}(a_{t}+b_{t})\right)-1\leq\mathfrak{L}+M$ ; so all terms in $N_{I,\mathfrak{L}+1-I}$, as elements of $\mathbb{C}[\lambda_{1,0},\lambda_{0,1}]$, have degree at most equal to $M+\mathfrak{L}$ ; the last step is to prove that the $\lambda_{I+k,M+\mathfrak{L}+1-I-k}$ are also of degree at most equal to $M+\mathfrak{L}$, and we will proceed by induction on $I$:

1. (a)

for $I=0$, all $\lambda_{I+k,M+\mathfrak{L}+1-I-k}$ for $0\leq k\leq M-1$ satisfy the two conditions $I+k<M$ and $I+k+M+\mathfrak{L}+1-I-k=M+\mathfrak{L}+1\neq 1$ so they are all prescribed to zero and their degree as element of $\mathbb{C}[\lambda_{1,0},\lambda_{0,1}]$ is at most equal to $M+\mathfrak{L}$ that ; 2. (b)

suppose that, for a given $I\in\mathbb{N}$, the $\lambda_{\tilde{I}+k,M+\mathfrak{L}+1-\tilde{I}-k}$ for all $\tilde{I}\in\mathbb{N}$ such that $\tilde{I}\leq I$ are also of degree at most equal to $M+\mathfrak{L}$ then it is clear from Equation (27) that $\lambda_{I+1+M,\mathfrak{L}-I-1}$ is also of degree at most equal to $M+\mathfrak{L}$.

This concludes the proof.

∎

Finally, since we are interested in the local approximation properties of GPWs, it is natural to study their Taylor expansion coefficients, and how they can be expressed as elements of $\mathbb{C}[\lambda_{1,0},\lambda_{0,1}]$. In particular we will find what is the link between the Taylor expansion coefficients of a GPW, $\partial_{x}^{i}\partial_{y}^{j}\varphi\left(x_{0},y_{0}\right)/(i!j!)$, and that of the corresponding PW, $(\lambda_{0,1})^{j}(\lambda_{1,0})^{i}/(i!j!)$.

Proposition 6**.**

Consider a point $(x_{0},y_{0})\in\mathbb{R}^{2}$, a given $q\in\mathbb{N}^{*}$, a given $M\in\mathbb{N}$, $M\geq 2$, a given set of complex-valued functions $\alpha=\{\alpha_{k,l}\in\mathcal{C}^{q-1}\text{ at }(x_{0},y_{0}),0\leq k+l\leq M\}$, and the corresponding partial differential operator $\mathcal{L}_{M,\alpha}$. Under Hypotheses 1 and 2 consider a solution to Problem (10) constructed thanks to Algorithm 1 with all the fixed values $\lambda_{i,j}$ such that $i<M$ and $i+j\neq 1$ set to zero, and

[TABLE]

for some $\theta\in\mathbb{R}$ and $\kappa\in\mathbb{C}^{*}$, and the corresponding $\displaystyle\varphi(x,y)=\exp\sum_{0\leq i+j\leq q+1}\lambda_{ij}(x-x_{0})^{i}(y-y_{0})^{j}$. Then for all $(i,j)\in\mathbb{N}^{2}$ such that $i+j\leq q+1$ the difference

[TABLE]

can be expressed as an element of $\mathbb{C}[\lambda_{1,0},\lambda_{0,1}]$ such that

• $\bullet$

its total degree satisfies ${\rm d}R_{i,j}\leq i+j-1$,

• $\bullet$

its coefficients only depend on $i$, $j$, and on the derivatives of the PDE coefficients $\alpha$ evaluated at $(x_{0},y_{0})$ but do not depend on $\theta$.

Proof.

Applying the chain rule introduced in Appendix A.2 to the GPW $\varphi$ one gets for all $(i,j)\in\mathbb{N}^{2}$,

[TABLE]

where $p_{s}((i,j),\mu)$ is the set of partitions of $(i,j)$ with length $\mu$:

[TABLE]

For each partition $(k_{l},(i_{l},j_{l}))_{l\in[\![1,s]\!]}$ of $(i,j)$, the corresponding product term, considered as an element of $\mathbb{C}[\lambda_{1,0},\lambda_{0,1}]$, has degree $\displaystyle Deg\ \prod_{l=1}^{s}(\lambda_{i_{l},j_{l}})^{k_{l}}=\sum_{l=1}^{s}k_{l}Deg\ \lambda_{i_{l},j_{l}}$. Combining Proposition 5 and the fact that $\lambda_{i,j}=0$ for all $(i,j)$ such that $1<i+j<M$, we can conclude that this degree is at most equal to

[TABLE]

The partition with two terms $(i,j)=j(0,1)+i(1,0)$ corresponds to the term $(\lambda_{0,1})^{j}(\lambda_{1,0})^{i}$, which is the leading term in $\partial_{x}^{i}\partial_{y}^{j}\varphi\left(x_{0},y_{0}\right)$. Indeed, any other partition will include at least one term such that $i_{l}+j_{l}>1$, and the degree corresponding to this term within the product is either $k_{l}\cdot 0$ or $k_{l}(i_{l}+j_{l}-1)$, and in both case it is at most equal to $k_{l}(i_{l}+j_{l})-1$. As a result, the degree of the product term in (29) is necessarily less than $\displaystyle\sum_{l=1}^{s}k_{l}(i_{l}+j_{l})=i+j$. So $R_{i,j}$, which is defined as the difference between $\partial_{x}^{i}\partial_{y}^{j}\varphi\left(x_{0},y_{0}\right)$ and its leading term $(\lambda_{0,1})^{j}(\lambda_{1,0})^{i}$, is as expected of degree less than $i+j$.

Finally, the coefficients of $R_{i,j}$ share the same property as the coefficients of $\lambda_{ij}$s from Propositions 5. ∎

Remark 3**.**

As mentioned in Remark 2, under the hypothesis $\alpha_{0,M}(x_{0},y_{0})\neq 0$, an algorithm very similar to Algorithm 1 would construct the polynomial coefficients of a GPW, fixing the values of $\{\lambda_{i,j},0\leq j\leq M-1,0\leq i\leq q+M-1-j\}$. The corresponding version of Proposition 6 could then be proved essentially by exchanging the roles of $i$ and $j$ in all the proofs.

3.3 Local set of GPWs

At this point for a given value of $\theta\in\mathbb{R}$ we can construct a GPW as a function $\varphi=\exp P$ where the polynomial $P$ is a solution to Problem (10) constructed thanks to Algorithm 1 with all the fixed values $\lambda_{i,j}$ such that $i<M$ and $i+j\neq 1$ set to zero, and

[TABLE]

This parameter $\theta$ is then equivalent to the direction a classical plane wave, while $|\kappa|$ is equivalent to the wave number of a classical plane wave, and $\theta$ will now be used to construct a set of GPWs. Under Hypotheses 1 and 2, by choosing $p$ different angles $\{\theta_{l},l\in\mathbb{N}^{*},l\leq p\}\in\mathbb{R}^{p}$, we can consider $p$ solutions to Problem (10) to construct $p$ GPWs.

Definition 8**.**

Consider a point $(x_{0},y_{0})\in\mathbb{R}^{2}$, a given $q\in\mathbb{N}^{*}$, a given $M\in\mathbb{N}$, $M\geq 2$, a given set of complex-valued functions $\alpha=\{\alpha_{k,l}\in\mathcal{C}^{q-1}\text{ at }(x_{0},y_{0}),0\leq k+l\leq M\}$, and the corresponding partial differential operator $\mathcal{L}_{M,\alpha}$. Let $p\in\mathbb{N}^{*}$ be the number of desired basis functions. Under Hypotheses 1 and 2, consider the normalization $\lambda_{i,j}$ such that $i<M$ and $i+j\neq 1$ set to zero, and

[TABLE]

The set of corresponding GPWs contructed from Algorithm 1 will be denoted hereafter by

[TABLE]

4 Interpolation properties

This section is restricted to operators of order $M=2$.

We now have built tools to turn to the interpolation properties of GPWs. In particular, since the GPWs are constructed locally, and will be defined separately on each mesh element, we focus on local interpolation properties. Given a partial differential operator $\mathcal{L}$, a point $(x_{0},y_{0})\in\mathbb{R}^{2}$ and an integer $n\in\mathbb{N}$, the question is whether we can find a finite dimensional space $\mathbb{V}_{h}\subset\mathcal{C}^{\infty}$, with the following property:

[TABLE]

that is to say there exists an element of $\mathbb{V}_{h}$ whose Taylor expansion at $(x_{0},y_{0})$ matches the Taylor expansion of $u$ at $(x_{0},y_{0})$ up to order $n$, for any solution $u$ of the PDE $\mathcal{L}u=0$. If $\{f_{i},i\in\mathbb{N}^{*},i\leq p\}$ is a basis of $\mathbb{V}_{h}$, this can be expressed in terms of linear algebra. Consider the vector space $\mathbb{F}$ and the matrix $\mathsf{M}\in\mathbb{C}^{(n+1)(n+2)/2\times p}$ defined as follows:

[TABLE]

[TABLE]

Then (30) is equivalent to

[TABLE]

and the choice of $p$, the number of basis functions, will be crucial to our study.

Our previous work on GPWs was focused on the Helmholtz equation, i.e. $\mathcal{L}=-\Delta+\beta(x,y)$, and in that case the classical PWs are exact solutions to the PDE if the coefficient is constant $\beta(x,y)=-\kappa^{2}$. However, even though the proof of the interpolation properties of GPWs relies strongly on that of classical PWs, it is not required, in order to obtain the GPW result, for classical PW to be solutions of the constant coefficient equation [13]. Indeed, what will be central to the proof that follows is the rank of the matrix $\mathsf{M}$ associated to a set of reference functions - not necessarily classical PWs - that are not required to satisfy any PDE. For the Helmholtz equation, the reference functions used in [13] were classical PWs if $\beta(x_{0},y_{0})<0$ and real exponentials if $\beta(x_{0},y_{0})>0$, and the structure of the proof provides useful guidelines for what follows.

4.1 Comments on a standard reference case

Interpolation properties of classical plane waves were already presented for instance in [13], and in [4], however the link between desired order of approximation $n$ and number $p$ of basis functions was simply provided as $p=2n+1$. We present here a new perspective, focusing on properties of trigonometric functions, to justify this choice. The corresponding set of trigonometric functions will constitute the reference case at the heart of the GPWs interpolation properties.

Definition 9**.**

Consider a given $n\in\mathbb{N}^{*}$ and a given $p\in\mathbb{N}^{*}$. Considering for some $\kappa\in\mathbb{R}^{*}$ a space $\mathbb{V}_{h}^{\kappa}=Span\{\exp\mathrm{i}\kappa(\cos\theta_{l}(x-x_{0})+\sin\theta_{l}(y-y_{0})),1\leq l\leq p,\theta_{l}\in[0,2\pi),\theta_{l_{1}}\neq\theta_{l_{2}}\ \forall l_{1}\neq l_{2}\}$ of classical PWs, we define the corresponding matrix (31) for the plane wave functions spanning $\mathbb{V}_{h}^{\kappa}$, denoted $\mathsf{M}^{C}$, as well as the reference matrix $\mathsf{M}^{R}$, by

[TABLE]

If we denote by $\mathsf{D}_{n}^{RC}=diag(d_{k}^{RC},k\text{ from }1\text{ to }n+1)$ the block diagonal matrix with blocks of increasing size $d_{k}^{RC}=(\mathrm{i}\kappa)^{k-1}I_{k}\in\mathbb{C}^{k\times k}$, it is evident that $\mathsf{M}^{C}_{n}=\mathsf{D}_{n}^{RC}\mathsf{M}^{R}_{n}$, therefore trigonometric functions are closely related to interpolation properties of PWs.

Consider the two sets of functions

[TABLE]

The first one, $\mathcal{F}_{n}$, is a set of $(n+1)(n+2)/2$ functions. The second one, $\mathcal{G}_{n}$, is a set of $2n+1$ linearly independent functions: indeed, any null linear combination of these functions $\sum_{-n\leq k\leq n}\nu_{k}\exp(ik\theta)$ would define a function $f(x)=\sum_{-n\leq k\leq n}\nu_{k}x^{k}$ that would be uniformly null on the circle $|x|=1$, implying that the polynomial $x^{n}.f(x)$ has an infinite number of roots ; hence all its coefficients $\nu_{k}$ are null. Moreover since

[TABLE]

we then have that $Span\ \mathcal{F}_{n}=Span\ \mathcal{G}_{n}$, and in particular the space spanned by $\mathcal{F}_{n}$ is of dimension $2n+1$.

Consider any matrix $\mathsf{A}^{\mathcal{F}}\in\mathbb{C}^{(n+1)(n+2)/2\times N_{p}}$ defined for some $\{\theta_{l}\}_{1\leq l\leq N_{p}}\in\left(\mathbb{R}\right)^{N_{p}}$, with $N_{p}>2n+1$, by

[TABLE]

Its rank is at most $2n+1$. This is a simple consequence of the fact that the dimension of $Span\ \mathcal{F}_{n}$ is $2n+1<(n+1)(n+2)/2$: indeed, this implies that there exists a matrix $\mathsf{C}\in\mathbb{C}^{((n+1)(n+2)/2-2n-1)\times(n+1)(n+2)/2}$ of rank $(n+1)(n+2)/2-2n-1$ such that

[TABLE]

and therefore $\mathsf{C}\mathsf{A}^{\mathcal{F}}=0_{((n+1)(n+2)/2-2n-1)\times N_{p}}$ ; as a result the $N_{p}$ columns of $\mathsf{A}^{\mathcal{F}}$ belong to the kernel of $\mathsf{C}$, which is of dimension $2n+1$; so the rank of $\mathsf{A}^{\mathcal{F}}$ is at most $2n+1$. In particular the matrix $\mathsf{M}_{n}^{R}$ introduced in Definition 9 is such a matrix $\mathsf{A}^{\mathcal{F}}$, and is therefore of rank at most $2n+1$.

We know that $\mathsf{M}^{C}_{n}=\mathsf{D}_{n}^{RC}\mathsf{M}^{R}_{n}$ and $\mathsf{D}_{n}^{RC}$ is non-singular, so $rk(\mathsf{M}_{n}^{C})=rk(\mathsf{M}^{R}_{n})$. The rank of $\mathsf{M}^{C}_{n}$ is at most equal to $2n+1$ for any choice of angles $\{\theta_{l}\in\mathbb{R},1\leq l\leq p\}$. It was previously proved in Lemma 2 from [13] that for $p=2n+1$ and directions such that $\{\theta_{l}\in[0,2\pi),1\leq l\leq p,l_{1}\neq l_{2}\Rightarrow\theta_{l_{1}}\neq\theta_{l_{2}}\}$ the matrix $\mathsf{M}_{n}^{C}$ has rank $2n+1$. A trivial corollary of this proof is that, for any choice of $p$ distinct angles in $[0,2\pi)$,

[TABLE]

In [13] we also proved that the space $\mathbb{F}$ for the constant coefficient Helmholtz operator is equal to the range of $\mathsf{M}_{n}^{C}$ for the corresponding wave number $\kappa$. As a direct consequence, a space $\mathbb{V}_{h}^{\kappa}=Span\{\exp\mathrm{i}\kappa(\cos\theta_{l}(x-x_{0})+\sin\theta_{l}(y-y_{0})),1\leq l\leq p\}$ for any choice of distinct angles in $[0,2\pi)$ satisfies the interpolation property (30) for the Helmholtz equation if and only if $p\geq 2n+1$.

4.2 Generalized Plane Wave case

In order to prove that a GPW space $Span\ \mathbb{V}_{\alpha,p,q}^{0}$ (introduced in Definition 8) satisfies the interpolation property (30), we will rely on Proposition 6 to study the rank of the matrix (31) built from GPWs. As in the Helmholtz case, the proof relates the GPW matrix to the reference matrix, but here via an intermediate transition matrix.

Definition 10**.**

Consider a point $(x_{0},y_{0})\in\mathbb{R}^{2}$, a given $q\in\mathbb{N}^{*}$, a given $M\in\mathbb{N}$, $M\geq 2$, a given set of complex-valued functions $\alpha=\{\alpha_{k,l}\in\mathcal{C}^{q-1}\text{ at }(x_{0},y_{0}),0\leq k+l\leq M\}$, and the corresponding partial differential operator $\mathcal{L}_{M,\alpha}$. For the corresponding set of GPWs, $\mathbb{V}_{\alpha,p,q}^{0}=\{\varphi_{l}:=\exp P_{l},\forall l\in\mathbb{N}^{*},l\leq p,\theta_{l}\in[0,2\pi),\theta_{l_{1}}\neq\theta_{l_{2}}\ \forall l_{1}\neq l_{2},\kappa\in\mathbb{C}^{*}\}$, we define the corresponding matrix (31), denoted $\mathsf{M}_{n}$, as well as the transition matrix $\mathsf{M}_{n}^{Tr}$, by

[TABLE]

We first relate the transition matrix $\mathsf{M}_{n}^{Tr}$ to the reference matrix $\mathsf{M}_{n}^{R}$.

Lemma 7**.**

Consider an open set $\Omega\subset\mathbb{R}^{2}$, $(x_{0},y_{0})\in\Omega$, a given $(M,n,p,q)\in(\mathbb{N}^{*})^{4}$, $M\geq 2$, and a given set of complex-valued functions $\alpha=\{\alpha_{k_{1},k_{2}}\in\mathcal{C}^{q-1}(\Omega),0\leq k_{1}+k_{2}\leq M\}$, the corresponding partial differential operator $\mathcal{L}_{M,\alpha}$ and set of GPWs $\mathbb{V}_{\alpha,p,q}^{0}$. There exists a block diagonal non-singular matrix $\mathsf{D}_{n}^{RT}$ such that $\mathsf{M}_{n}^{Tr}=\mathsf{D}_{n}^{RT}\mathsf{M}^{R}_{n}$, independently of the number $p$ of GPWs in $\mathbb{V}_{\alpha,p,q}^{0}$.

Proof.

As long as there are four complex numbers $a,b,c,d$ such that

[TABLE]

then the diagonal blocks of $\mathsf{D}_{n}^{RT}=diag(d_{K}^{RT},K\text{ from }0\text{ to }n)$ of increasing size $d_{K}^{RT}\in\mathbb{C}^{(K+1)\times(K+1)}$ can be built thanks to the following binomial formula

[TABLE]

since the coefficient of this linear combination of trigonometric functions are independent on $l$. ∎

The following step is naturally to relate the GPW matrix $\mathsf{M}_{n}$ to the reference matrix $\mathsf{M}_{n}^{R}$.

Proposition 7**.**

Consider an open set $\Omega\subset\mathbb{R}^{2}$, $(x_{0},y_{0})\in\Omega$, a given $(M,n,p,q)\in(\mathbb{N}^{*})^{4}$, $M\geq 2$, $q\geq n-1$, and a given set of complex-valued functions $\alpha=\{\alpha_{k_{1},k_{2}}\in\mathcal{C}^{\max(n,q-1)}(\Omega),0\leq k_{1}+k_{2}\leq M\}$, the corresponding partial differential operator $\mathcal{L}_{M,\alpha}$ and set of GPWs $\mathbb{V}_{\alpha,p,q}^{0}$. There exists a lower triangular matrix $L_{n}^{R}$, whose diagonal coefficients are equal all non-zero and whose other non-zero coefficients depend only on derivatives of the PDE coefficients $\alpha$ evaluated at $(x_{0},y_{0})$, such that

[TABLE]

As a consequence $rk(\mathsf{M}_{n})=rk(\mathsf{M}_{n}^{R})$ independently of the number $p$ of GPWs in $\mathbb{V}_{\alpha,p,q}^{0}$, and both $\|\mathsf{L}_{n}^{R}\|$ and $\|(\mathsf{L}_{n}^{R})^{-1}\|$ are bounded by a constant depending only on the PDE coefficients $\alpha$.

Remark 4**.**

If $n=1$, then the various matrices $\mathsf{M}$ belong to $\mathbb{C}^{3\times 3}$, and we have $\mathsf{M}_{n}=\mathsf{M}_{n}^{Tr}$ independently of the value of $q$.

Proof.

Let’s first relate $\mathsf{M}_{n}$ to $\mathsf{M}_{n}^{Tr}$. The polynomials $R_{i,j}\in\mathbb{C}[X,Y]$ obtained in Proposition 6 have degree $\textrm{d}R_{i,j}\leq i+j-1$ and satisfy

[TABLE]

In order to apply this to all entries in the matrix $\mathsf{M}_{n}$, it is sufficient for $q$ to satisfy $n\leq q+1$, which explains the assumption on $q$. Therefore each entry $(i,j)$ of the matrix $\mathsf{M}_{n}$ can be written as the sum of the $(i,j)$ entry of $\mathsf{M}_{n}^{Tr}$ and a linear combination of entries $(k,j)$ of $\mathsf{M}_{n}^{Tr}$ for $k<i$. In other words, the existence of a lower triangular matrix $\mathsf{L}_{n}^{T}$, whose diagonal coefficients are 1 and whose other non-zero coefficients depend only on the derivatives of the coefficients $\alpha$ evaluated at $(x_{0},y_{0})$, such that $\mathsf{M}_{n}=\mathsf{L}_{n}^{T}\cdot\mathsf{M}_{n}^{Tr}$ is guaranteed since the coefficients of $R_{i,j}$ are independent of $l$ and any monomial in $R_{i,j}(\lambda_{1,0},\lambda_{1,0})$ has a degree lower than $i+j$.

As a consequence, the existence of $\mathsf{L}_{n}^{R}$ is guaranteed by Lemma 7 since $\mathsf{L}_{n}^{R}:=\mathsf{L}_{n}^{T}\cdot\mathsf{D}_{n}^{RT}$ satisfies the desired properties. ∎

Everything is now in place to state and finally prove the necessary and sufficient condition on the number $p$ of GPWs for the space $\mathbb{V}_{\alpha,p,q}^{0}$ to satisfy the interpolation property (30). We here turn to the specific case of second order operators.

Theorem 1**.**

Consider an open set $\Omega\subset\mathbb{R}^{2}$, $(x_{0},y_{0})\in\Omega$, $M=2$, a given $(n,p,q)\in\mathbb{(}N^{*})^{3}$, $n\geq M$, $q\geq n-1$ and a given set of complex-valued functions $\alpha=\{\alpha_{k_{1},k_{2}}\in\mathcal{C}^{n}(\Omega),0\leq k_{1}+k_{2}\leq M\}$, the corresponding partial differential operator $\mathcal{L}_{M,\alpha}$ and set of GPWs $\mathbb{V}_{\alpha,p,q}^{0}$. The space $\mathbb{V}_{h}^{G}:=span\mathbb{V}_{\alpha,p,q}^{0}$ satisfies the property

[TABLE]

if and only if $p\geq 2n+1$.

Proof.

Combining (32) with Proposition 7, if $q\geq n-1$, we see immediately that $rk(\mathsf{M}_{n})=2n+1$ if and only if $p\geq 2n+1$.

It is then sufficient to prove that the space

[TABLE]

satisfies $\mathbb{F}_{\alpha}\subset\mathcal{R}(\mathsf{M}_{n})$, the range of $\mathsf{M}_{n}$. To this end, we now define the space

[TABLE]

We can now see that

• •

$\mathcal{R}(\mathsf{M}_{n})\subset\mathfrak{K}$ independently of the value of $p$, since by construction of GPWs, as long as $q\geq n-1$, each column of $\mathsf{M}_{n}$ belongs to $\mathfrak{K}$;

• •

$\mathbb{F}_{\alpha}\subset\mathfrak{K}$, by definition of $\mathbb{F}_{\alpha}$;

• •

$\dim\mathfrak{K}=2n+1$, since - from the condition involving the Taylor expansion coefficients of $\mathcal{L}_{2,\alpha}f$ of order up to $n-2$ at $x_{0},y_{0})$ set to zero - $\mathfrak{K}\subset\mathbb{C}^{(n+1)(n+2)/2}$ is the kernel of a matrix $\mathsf{A}\in\mathbb{C}^{n(n-1)/2\times(n+1)(n+2)/2}$ with

[TABLE]

so that $\mathsf{A}$ is of maximal rank while its kernel has dimension $\frac{(n+1)(n+2}{2}-\frac{n(n-1)}{2}=2n+1$.

Therefore, if $p\geq 2n+1$, we obtain that $\mathcal{R}(\mathsf{M}_{n})=\mathfrak{K}$ and as a consequence $\mathbb{F}_{\alpha}\in\mathcal{R}(\mathsf{M}_{n})$ as expected. This concludes the proof. ∎

The necessary and sufficient condition on the number $p$ of GPWs for the space $\mathbb{V}_{\alpha,p,q}^{0}$ to satisfy the interpolation property (30) when $M>2$ are still unknown.

Remark 5**.**

As in [13], the theorem holds for the Helmholtz equation with sign changing.

5 Numerical experiments

In [13], GPWs where constructed and studied for the Helmholtz equation (1) with a variable and sign-changing coefficient $\beta$. The numerical experiments presented there were restricted to the Helmholtz equation at one point $(x_{0},y_{0})\in\mathbb{R}^{2}$, but considered a propagative case i.e. $\beta(x_{0},y_{0})<0$, an evanescent case i.e. $\beta(x_{0},y_{0})>0$, a cut-off case i.e. $\beta(x_{0},y_{0})=0$. They also considered a case not covered by the convergence theorem, but important for future applications: considering GPWs centered at points $(x_{0},y_{0})$ at a distance $h$ from the cut-off.

Here, we are interested in illustrating the results presented in Theorem 1. Since the well known case of classical PW for the constant-coefficient Helmholtz equation is included by the hypotheses of the theorem, we cannot expect any improvement on the required number of basis functions $p$. However, we are interested in exploring the impact of the order of approximation $q$ on the convergence of (34), in particular for anisotropic problems.

5.1 Test cases

We propose here four different test cases. Each test case consists of a partial differential operator of second order $\mathcal{L}$, an exact solution $u$ of the equation $\mathcal{L}u=0$, as well as a computational domain $\Omega\subset\mathbb{R}^{2}$, such that Hypotheses 1 and 2 hold at all $(x_{0},y_{0})\in\Omega$. The characteristics of the partial differential operators that we consider here are:

• •

polynomial coefficients $\alpha$,

• •

non-polynomial coefficients $\alpha$,

• •

anisotropy in the first order terms as $\overrightarrow{a}(x,y)\cdot\nabla$ for a vector-valued function $\overrightarrow{a}$;

• •

anisotropy in the second order terms as $\nabla\cdot(A(x,y)\nabla)$ for a matrix-valued function $A$.

We consider one partial differential operator is isotropic with polynomial coefficients:

[TABLE]

We have $\mathcal{L}_{Ad}u_{Ad}=0$ on $\mathbb{R}^{2}$, all the coefficients of $\mathcal{L}_{Ad}$ belong to $\mathcal{C}^{\infty}(\mathbb{R}^{2})$ and the coefficients $\{\alpha_{k,2-k}^{Ad};k=0,1,2\}$ satisfy

[TABLE]

so $\mathcal{L}_{Ad}$ satisfies Hypotheses 1 and 2 on $\mathbb{R}^{2}$. Note that the sign of the coefficient $\alpha_{0,0}^{Ad}(x,y)=2(x+y)$ changes in the computational domain along the curve $x+y=0$.

We consider a partial differential operator with non-polynomial coefficients of the terms of order $1$ and [math], and anisotropy in the first order term:

[TABLE]

We have $\mathcal{L}_{Jc}u_{Jc}=0$ on $(0,\infty)\times\mathbb{R}$, all the coefficients of $\mathcal{L}_{Jc}$ belong to $\mathcal{C}^{\infty}(\mathbb{R}^{+}\times\mathbb{R})$ and the coefficients $\{\alpha_{k,2-k}^{Jc};k=0,1,2\}$ satisfy

[TABLE]

so $\mathcal{L}_{Jc}$ satisfies Hypotheses 1 and 2 as long as $x>0$.

We consider a partial differential operator with polynomial coefficients and anisotropy in the first and second order terms:

[TABLE]

We have $\mathcal{L}_{JJ}u_{JJ}=0$ on $(\mathbb{R}^{+})^{2}$, all the coefficients of $\mathcal{L}_{JJ}$ belong to $\mathcal{C}^{\infty}((\mathbb{R}^{+})^{2})$ and the coefficients $\{\alpha_{k,2-k}^{JJ};k=0,1,2\}$ satisfy

[TABLE]

so $\mathcal{L}_{JJ}$ satisfies Hypotheses 1 and 2 as long as $xy\neq 0$.

Finally we consider a partial differential operator with non-polynomial coefficients and anisotropy in the second order term:

[TABLE]

We have $\mathcal{L}_{cs}u_{cs}=0$ on $\mathbb{R}^{2}$, all the coefficients of $\mathcal{L}_{cs}$ belong to $\mathcal{C}^{\infty}(\mathbb{R}^{2})$ and the coefficients $\{\alpha_{k,2-k}^{cs};k=0,1,2\}$ satisfy

[TABLE]

so $\mathcal{L}_{cs}$ satisfies Hypotheses 1 and 2 on $\mathbb{R}^{2}$.

5.2 Implementation of the construction algorithm

For a linear second order operator

[TABLE]

the associated operator $\mathfrak{L}_{2,\alpha}^{A}$ is defined by

[TABLE]

The implementation of Algorithm 1 simply requires, at each level $\mathfrak{L}$, the evaluation of $\{N_{I,\mathfrak{L}-I},0\leq I\leq\mathfrak{L}\}$ to apply formula (23). At each level $\mathfrak{L}$ the coefficeints $\{\mu_{ij},(i,j)\in\mathbb{N}^{2},i+j\leq q+1\}$ of $Q_{\mathfrak{L}}:=\sum_{0\leq i+j\leq M+\mathfrak{L}-1}\lambda_{i,j}(x-x_{0})^{i}(y-y_{0})^{j}$ are computed as

[TABLE]

and for $0\leq I\leq\mathfrak{L}$ the different contributions to $N_{I,\mathfrak{L}-I}$ can be described as follows:

• •

the linear contributions from first order terms $T_{3}$

[TABLE]

• •

the non-linear contributions from the terms $T_{2}$

[TABLE]

• •

the linear contributions from the second order terms $T_{1}$

[TABLE]

• •

the contribution from the zeroth order term $\alpha_{0,0}$

[TABLE]

Moreover, all experiments are conducted with the following choice of angles $\theta_{l}$ and $\kappa$ parameters to build the GPW space $\mathbb{V}_{\alpha,p,q}^{0}$:

[TABLE]

5.3 Numerical results

The $h$-convergence results presented in Theorem 1 are stated as local properties at a given point. In order to illustrate them, for each test case, we consider the following procedure.

• •

At each of 50 random points $(x_{0},y_{0})$ in the computational domain $\Omega$

Construct the set of GPWs from Algorithm 1 with the normalization proposed in section 3. 2. 2.

Compute $u_{a}$ the linear combination of GPWs studied in the theorem’s proof, matching its Taylor expansion to that of the exact solution.

• •

Estimate as a function of $h$ the maximum $L^{\infty}$ error on a disk of radius $h$ centered at the random point: $\max_{(x_{0},y_{0})\in\Omega}\|u-u_{a}\|_{L^{\infty}(\{(x,y)\in\mathbb{R}^{2},|(x,y)-(x_{0},y_{0})|<h\})}$.

We always consider a space $\mathbb{V}_{\alpha,p,q}^{0}$ of $p=2n+1$ GPWs. According to the theorem, we expect to observe convergence of order $n+1$ if the approximation parameter $q$ in the construction of the basis functions is at least equal to $n-1$. For each of the four test cases proposed, we present: on the one hand results for $n$ from $1$ to $5$ with $q=\max(1,n-1)$ (Left panel); on the other hand results for $q$ from $1$ to $4$ with $n=4$ (Right panel). Hence with the first choice of parameters the theorem predicts convergence of order $n+1$, while with the second choice the theorem does not cover these cases.

The results are presented in Figure 4 for the approximation of $u_{Ad}$, Figure 5 for the approximation of $u_{Jc}$, Figure 6 for the approximation of $u_{JJ}$, and Figure 7 for the approximation of $u_{cs}$. We observe on Figures 4 and 7 the effect of the large condition number of the matrix $\mathsf{M}_{n}$: even though the expected orders of convergence are observed for large values of $h$, the error stagnates at a threshold for smaller values of $h$. We also observe, on the left panels of Figures 4, 5, 6 and 7, that for all of our test cases the constant $C(\Omega,n)$ from (34) in Theorem 1 does not seem to depend on $n$, even though the Theorem predicts that it does.

We summarize in the following table the orders of convergence observed, always using $\mathbb{V}_{\alpha,p,q}^{0}$ with $p=2n+1$. The bold entries correspond to cases covered by Theorem 1 i.e. $n+1$ for $q\leq n-1$, and the red entries correspond to cases with order of convergence observed higher than the theorem predicts.

[TABLE]

We can see from this table that in all cases covered by the theorem, we observe a convergence rate equal or slightly better than predicted. But it would seem that the hypotheses of the theorem are sharp.

6 Conclusion

In this work we have considered local properties in the neighborhood of a point $(x_{0},y_{0})\in\mathbb{R}^{2}$, for an operator $\mathcal{L}_{M,\alpha}$. To summarize, we followed the steps announced in the introduction:

construction of GPWs $\varphi$ such that $\mathcal{L}_{M,\alpha}\varphi(x,y)=O\left(\|(x,y)-(x_{0},y_{0})\|^{q}\right)$

1. (a)

choose an ansatz for $\displaystyle\varphi(x,y)=\exp\sum_{0\leq i+j\leq dP}\lambda_{ij}(x-x_{0})^{i}(y-y_{0})^{j}$ 2. (b)

identify the corresponding $N_{dof}=\frac{(dP+1)(dP+2)}{2}$ degrees of freedom, and $N_{eqn}=\frac{q(q+1)}{2}$ constraints, namely respectively

[TABLE] 3. (c)

for $dP=q+M-1$, the number of degrees of freedom is $N_{dof}=\frac{(M+q)(M+q+1)}{2}>N_{eqn}$ and this ensures that there are linear terms in all the constraints 4. (d)

identify $N_{dof}-N_{eqn}=Mq+\frac{M(M+1)}{2}$ additional constraints, namely

[TABLE]

to obtain a global system that can be split into a hierarchy of linear triangular subsystems 5. (e)

compute the remaining $N_{eqn}$ degrees of freedom by forward substitution for each triangular subsystem, therefore at minimal computational cost 2. 2.

interpolation properties

1. (a)

thanks to the normalization, in particular $\{\lambda_{i,j}=0,(i,j)\in\mathbb{N}^{2},i+j<M+q\text{ and }i<M,i+j\neq 1\}$, study the properties of the remaining $N_{eqn}$ degrees of freedom, that is $\{\lambda_{i,j},(i,j)\in\mathbb{N}^{2},i+j<M+q\text{ and }i\geq M\}$, with respect to $(\lambda_{1,0},\lambda_{0,1})$ 2. (b)

identify a simple reference case depending only on two parameters, that is basis functions $\displaystyle\phi(x,y)=\exp\lambda_{1,0}(x-x_{0})+\lambda_{0,1}(y-y_{0})$ depending only on the choice of $(\lambda_{1,0},\lambda_{0,1})$, independently of $\phi$ being an exact solution to the constant coefficient equation 3. (c)

study the interpolation properties of this reference case with classical PW techniques 4. (d)

relate the general case to the reference case thanks to 2a 5. (e)

prove the interpolation properties of the GPWs from those of the reference case

This construction process guarantees that the GPW function $\varphi$ satisfies the approximate Trefftz property $\mathcal{L}_{M,\alpha}\varphi(x,y)=O\left(\|(x,y)-(x_{0},y_{0})\|^{q}\right)$ independently of the normalization, that is the values chosen for $\{\lambda_{i,j},(i,j)\in\mathbb{N}^{2},i+j<M\}$, while the proof of interpolation properties heavily rely on the normalization.

This work focuses on interpolation of solutions of a PDE, and is limited to local results, in the neighborhood of a given point. In order to address the convergence of a numerical method for a boundary value problem on a domain $\Omega$ with a GPW-discretized Trefftz method, on a mesh $\mathcal{T}_{h}$ of $\Omega$, we will consider a space $\mathbb{V}_{h}$ of GPWs built element-wise, at the centroid $(x_{0},y_{0})=(x_{K},y_{K})$ of each element $K\in\mathcal{T}_{h}$, to study interpolation properties on $\Omega$. In particular, meshing the domain $\Omega$ to respect any discontinuity in the coefficients, the interpolation error on $\Omega$, $\|(I-P_{\mathbb{V}_{h}})\|$, will converge at the same order as the local interpolation error on each element, and the crucial point will be to investigate the behavior of the constant $C(\Omega,n)$ from Theorem 1. Related computational aspects of the construction of GPWs proposed in this work are currently under study.

Appendix A Chain rule in dimension 1 and 2

For the sake of completeness, this section is dedicated to describing the formula to derive a composition of two functions, in dimensions one and two. A wide bibliography about this formula is to be found in [25]. It is linked to the notion of partition of an integer or the one of a set. The 1D version is not actually used in this work but is displayed here as a comparison with a 2D version, mainly concerning this notion of partition.

A.1 Faa Di Bruno Formula

Faa Di Bruno formula gives the $m$th derivative of a composite function with a single variable. It is named after Francesco Faa Di Bruno, but was stated in earlier work of Louis F.A. Arbogast around 1800, see [6].

If $f$ and $g$ are functions with sufficient derivatives, then

[TABLE]

where the sum is over all different solutions in nonnegative integers $(b_{k})_{k\in[\![1,m]\!]}$ of $\sum_{k}kb_{k}=m$. These solutions are actually the partitions of $m$.

A.2 Bivariate version

The multivariate formula has been widely studied, the version described here is the one from [5] applied to dimension $2$. A linear order on $\mathbb{N}^{2}$ is defined by: $\forall(\mu,\nu)\in\left(\mathbb{N}^{2}\right)^{2}$, the relation $\mu\prec\nu$ holds provided that

$\mu_{1}+\mu_{2}<\nu_{1}+\nu_{2}$; or 2. 2.

$\mu_{1}+\mu_{2}=\nu_{1}+\nu_{2}$ and $\mu_{1}<\nu_{1}$.

If $f$ and $g$ are functions with sufficient derivatives, then

[TABLE]

where the partitions of $(i,j)$ are defined by the following sets: $\forall\mu\in[\![1,i+j]\!]$, $\forall s\in[\![1,i+j]\!]$, $p_{s}((i,j),\mu)$ is equal to

[TABLE]

See [11] for a proof of the formula interpreted in terms of collapsing partitions.

Appendix B Faa di Bruno

The multivariate formula has been widely studied, the version described here is the one from [5] applied to dimension $2$. A linear order on $\mathbb{N}^{2}$ is defined by: $\forall(\mu,\nu)\in\left(\mathbb{N}^{2}\right)^{2}$, the relation $\mu\prec\nu$ holds provided that

$\mu_{1}+\mu_{2}<\nu_{1}+\nu_{2}$; or 2. 2.

$\mu_{1}+\mu_{2}=\nu_{1}+\nu_{2}$ and $\mu_{1}<\nu_{1}$.

If $f$ and $g$ are functions with sufficient derivatives, then

[TABLE]

[TABLE]

where the partitions of $(i,j)$ are defined by the following sets: $\forall\mu\in[\![1,i+j]\!]$, $\forall s\in[\![1,i+j]\!]$, $p_{s}((i,j),\mu)$ is equal to

[TABLE]

Note that $s$ is the number of different terms appearing in the product, while $\mu$ is the number of terms in the product counting multiplicity, $k_{m}$ is the multiplicity of the $m$th term in the product, while $p_{s}$ represents the possible partitions of $(i,j)$.

Note that since $k_{m}>0$, the condition $\sum_{m=1}^{s}k_{m}=\mu$ implies that $\mu=\sum_{m=1}^{s}k_{m}\geq\sum_{m=1}^{s}1=s$.

Appendix C Polynomial formulas

Here are two important comments. The first one concerns the product of polynomials. Assume that $\min(D_{1},D_{2})\geq q$. Then the product of two polynomials, respectively of degree $D_{1}$ and $D_{2}$, satisfies:

[TABLE]

Since in particular the summation indices are such that $0\leq\tilde{i}\leq i$, $0\leq i-\tilde{i}\leq i$, $0\leq\tilde{j}\leq j$, and $0\leq j-\tilde{j}\leq j$, the only coefficients $p_{i,j}$ and $q_{i,j}$ appearing in the $(I_{0},J_{0})$ coefficient of the product have a length of the multi-index $i+j\leq I_{0}+J_{0}$. As a consequence, the only coefficients of several polynomials appearing in the $(I_{0},J_{0})$ coefficient of the product these several polynomials have a length of the multi-index $i+j\leq I_{0}+J_{0}$. The second comment turns to the derivative of a polynomial:

[TABLE]

In particular the only coefficients $p_{i,j}$ appearing in the $(I_{0},J_{0})$ coefficient of the derivative has a length of the multi-index $i+j=I+J+I_{0}+J_{0}$.