Method for solving hyperbolic systems with initial data on non-characteristic manifolds with applications to the shallow water wave equations

TL;DR

This paper introduces a new method for transforming initial data on non-characteristic manifolds for hyperbolic PDEs, with applications to modeling long wave run-up in fluid mechanics, particularly in inclined bays.

Contribution

The paper presents a novel method for converting initial conditions on non-characteristic manifolds into standard initial conditions for hyperbolic systems, enabling new solutions in fluid mechanics.

Findings

Effective transformation method for hyperbolic PDE initial data

Complete solution for long wave run-up in inclined bays

Application to fluid mechanics problems involving nonzero initial velocity

Abstract

We are concerned with hyperbolic systems of order-one linear PDEs originated on non-characteristic manifolds. We put forward a simple but effective method of transforming such initial conditions to standard initial conditions (i.e. when the solution is specified at an initial moment of time). We then show how our method applies in fluid mechanics. More specifically, we present a complete solution to the problem of long waves run-up in inclined bays of arbitrary shape with nonzero initial velocity.

Equations

$${\mathbf{u}}_t+\sum _{i=1}^n{A}_i{\mathbf{u}}_{x_i}+A_0\mathbf{u}$$

$${\mathbf{u}|}_{\Gamma }$$

$$\Gamma =\{(\mathbf{x},\tau (\mathbf{x}))|\mathbf{x}\in {\mathbb{R}}^n\}$$

$$A\left(\mathbf{x},t,\xi \right):=\sum _{i=1}^n{\xi }_i{A}_i\left(\mathbf{x},t\right),\xi =\left({\xi }_i\right)\in {\mathbb{R}}^n,$$

$$\det \left\{I-A\left(\mathbf{x},t,\nabla \tau \left(\mathbf{x}\right)\right)\right\}\ne 0\text{in }{\mathbb{R}}^{n+1}.$$

$$D\mathbf{u}:=\left[\sum _{i=1}^n{A}_i\left(\mathbf{x}\right)\frac{\partial }{\partial x_i}+A_0\left(\mathbf{x}\right)\right]u.$$

$$\left\{\begin{array}[c]{c}\mathbf{u}_{t}=-D\mathbf{u}\\ \mathbf{u}|_{\Gamma}=\mathbf{g}\end{array}\right..$$

$$\left\{\begin{array}[c]{c}\mathbf{u}_{t}=-D\mathbf{u}\\ \mathbf{u}|_{t=0}=\widetilde{\mathbf{g}}\end{array}\right.$$

$$\mathbf{u}(\mathbf{x},0)=\sum _{k=0}^j\frac{1}{k!}\frac{{\partial }^k\mathbf{u}(\mathbf{x},t)}{\partial t^k}(-t{)}^k+E_j\left(\mathbf{x},t\right)$$

$${\mathbf{u}\left(\mathbf{x},t\right)|}_{t=0}$$

$$=:{\mathbf{g}}_j(\mathbf{x})+{E_j|}_{\Gamma }.$$

$$\mathbf{u}{|}_{t=0}={\mathbf{g}}_j+O\left(\epsilon \right).$$

$${\mathbf{g}}_0(\mathbf{x})=\mathbf{u}(\mathbf{x},t){|}_{\Gamma }=\mathbf{g}\left(\mathbf{x}\right)\text{ (IC in (refIC)).}$$

$${{\mathbf{u}}_t|}_{\Gamma }=-{(D\mathbf{u})|}_{\Gamma }.$$

$$D\mathbf{g}$$

$$=\sum _{i=1}^n{A}_i\left({{\mathbf{u}}_{x_i}|}_{\Gamma }+{\tau }_{x_i}{{\mathbf{u}}_t|}_{\Gamma }\right)+A_0\mathbf{g}\text{(by the chain rule)}$$

$$={\left(\sum _{i=1}^n{A}_i{\mathbf{u}}_{x_i}+A_0\mathbf{u}\right)|}_{\Gamma }+\sum _{i=1}^n{\tau }_{x_i}{A}_i{{\mathbf{u}}_t|}_{\Gamma }$$

$$={\left(D\mathbf{u}\right)|}_{\Gamma }-\left(\sum _{i=1}^n{\tau }_{x_i}{A}_i\right){\left(D\mathbf{u}\right)|}_{\Gamma }\text{ (by (ref1st coeff))}$$

$$=\left(I-\sum _{i=1}^n{\tau }_{x_i}{A}_i\right){\left(D\mathbf{u}\right)|}_{\Gamma }.$$

$${\left(D\mathbf{u}\right)|}_{\Gamma }={\left(I-\sum _{i=1}^n{\tau }_{x_i}{A}_i\right)}^{-1}D\mathbf{g},$$

$${{\mathbf{u}}_t|}_{\Gamma }=-{\left(I-A\right)}^{-1}D\left(\mathbf{u}{|}_{\Gamma }\right)=-{\left(I-A\right)}^{-1}D\mathbf{g},$$

$$A:=\sum _{i=1}^n{\tau }_{x_i}\left(\mathbf{x}\right)A_i\left(\mathbf{x}\right)=A\left(\mathbf{x},\nabla \tau \left(\mathbf{x}\right)\right).$$

$${\left({\mathbf{u}}_t\right)}_t=-D\left({\mathbf{u}}_t\right)$$

$${{\mathbf{u}}_t|}_{\Gamma }=-{\left(I-A\right)}^{-1}D\mathbf{g}.$$

$$\left({\frac{{\partial }^k\mathbf{u}}{\partial t^k}|}_{\Gamma }\right)={\left(-{\left(I-A\right)}^{-1}D\right)}^k\mathbf{g},k=0,1,2,...$$

$${\mathbf{g}}_j=\mathbf{g}+\sum _{k=0}^j\frac{{\tau }^k}{k!}{\left({\left(I-A\right)}^{-1}D\right)}^k\mathbf{g}$$

$${\mathbf{u}}_t+\sum _{i=1}^n{A}_i\left(\mathbf{x}\right){\mathbf{u}}_{x_i}+A_0\left(\mathbf{x}\right)\mathbf{u}=0\text{ in }{\mathbb{R}}^{n+1},$$

$${\mathbf{u}|}_{\Gamma }=\mathbf{g}\left(\mathbf{x}\right).$$

$$\underset{\mathbf{x}\in {\mathbb{R}}^n}{\max }\Vert \frac{{\tau }^{j+1}\left(\mathbf{x}\right)}{\left(j+1\right)!}{\left\{{\left[I-A\left(\mathbf{x},\nabla \tau \left(\mathbf{x}\right)\right)\right]}^{-1}D\right\}}^{j+1}g\left(\mathbf{x}\right)\Vert <\epsilon ,$$

$$A\left(\mathbf{x},\nabla \tau \left(\mathbf{x}\right)\right):=\sum _{i=1}^n{\tau }_{x_i}\left(\mathbf{x}\right)A_i\left(\mathbf{x}\right),D:=\sum _{i=1}^n{A}_i\left(\mathbf{x}\right)\frac{\partial }{\partial x_i}+A_0\left(\mathbf{x}\right).$$

$${\tilde{\mathbf{u}}|}_{t=0}=\mathbf{g}+\sum _{k=0}^j\frac{{\tau }^k}{k!}{\left({\left(I-A\right)}^{-1}D\right)}^k\mathbf{g}$$

$$\left\{\begin{array}[c]{ccc}\eta_{t}+\left(1+\eta_{x}\right)u+c\left(x+\eta\right)u_{x}=0&&\text{(continuity equation)}\\ u_{t}+uu_{x}+\eta_{x}=0&&\text{(momentum equation)}\\ \eta\left(x,0\right)=\eta_{0}\left(x\right),\ \ u\left(x,0\right)=u_{0}\left(x\right)&&\text{(initial conditions)}\end{array}\right.,$$

$$\begin{array}[c]{cc}\varphi\left(\sigma,\tau\right)=u\left(x,t\right),&\psi\left(\sigma,\tau\right)=\eta\left(x,t\right)+u^{2}\left(x,t\right)/2,\\ \sigma=x+\eta\left(x,t\right),&\tau=t-u\left(x,t\right),\end{array}$$

Full text

The author is supported in part by the NSF grant DMS-1716975.

1. Introduction

Let $A_{i}=A_{i}\left(\mathbf{x},t\right),i=0,1,...,n$ be known $m\times m$ matrix-valued functions of $\mathbf{x}=\left(x_{i}\right)\in\mathbb{R}^{n}$ and $t\in\mathbb{R}$; $\mathbf{f}\left(\mathbf{x},t\right),\mathbf{g}\left(\mathbf{x}\right)$ are known $\mathbb{R}^{m}$-valued functions; and $\mathbf{u}\left(\mathbf{x},t\right)$ is an unknown $\mathbb{R}^{m}$-valued function. We are concerned with the Cauchy problem on manifolds for hyperbolic systems of linear first-order PDEs

[TABLE]

where

[TABLE]

is a n-dimensional manifold in $\mathbb{R}^{n+1}$ defined by a known scalar function $\tau(\mathbf{x})$. In particular, if $n=1$ then $\Gamma$ is a curve.

Recall that the system (1.1) is called hyperbolic if the matrix

[TABLE]

is diagonalizable for each $\mathbf{x},\mathbf{\xi}\in\mathbb{R}^{n}$ and $t\in\mathbb{R}$ with real eigenvalues (playing roles of characteristic velocities) (see, e.g. [7]).

If $\tau\left(\mathbf{x}\right)=0$ then (1.2) becomes the standard initial condition (IC) $\left.\mathbf{u}\right|_{t=0}=\mathbf{g}$ for which the initial value problem (IVP) is well-posed for a large class of $A_{i},\mathbf{f},\mathbf{g}$ and is very well studied. If $\tau\left(\mathbf{x}\right)\neq 0$ then (1.1)-(1.2) is well-posed [7] only if

[TABLE]

Such manifolds are called non-characteristic, i.e. tangent planes to $\Gamma$ and to characteristic manifolds are nowhere parallel. The problem however is that solving (1.1)-(1.2) becomes a notoriously hard problem even for constant $A_{i}$. On the other hand, such problems routinely appear while linearizing certain quasi-linear PDEs by the hodograph transform, a technique involving switching independent and dependent variables (see e.g. [5] for an extensive list of PDEs linearizable by the hodograph and [14] for the current literature). As it is simply put in [8] “Interchanging the dependent and independent variables simplifies the governing equations, but complicates the boundary/initial conditions.” This is a serious issue deterring from using the hodograph transform and when they do use it the problem of dealing with entangled initial/boundary conditions is solved on an ad hoc basis relying on various intuitive assumptions with no error analysis.

In the present note we put forward an elementary but very effective method of converting conditions on manifolds to standard conditions. Once it is done, one can proceed with any applicable method of solving the standard IBVP. Note our approach easily extends to higher order linear hyperbolic PDEs since they can be written as systems. We demonstrate effectiveness of our approach on the run-up problem of tsunami waves, a well known problem of fluid mechanics, by improving on (in chronological order) [4], [10], [11], [2], and [13] (see also the literature cited therein).

We owe the following comment to one of the referees. The book [17] provides more examples of the applicability of our work to other nonlinear hyperbolic equations. In particular, to Section 2.2 and Chapter 5 of [17] where our method could be used to eliminate some of the simplifications used to form initial conditions post applying a hodograph transform.

2. The method of data projection

In this section we introduce what we call the* *method of data projection. For simplicity we assume that $A_{i}$ depend only on $\mathbf{x}$ and the nonhomogeneous term $\mathbf{f}$ is absent. It will be clear however that the general case merely results in more complicated formulas. Introduce the following differential operator

[TABLE]

Our IVP then reads

[TABLE]

Our idea is, given any accuracy $\epsilon$, we find a standard initial condition $\mathbf{u}|_{t=0}=\widetilde{\mathbf{g}}$ such that the solution $\widetilde{\mathbf{u}}$ to

[TABLE]

is within $O\left(\varepsilon\right)$ from the actual solution to (2.1) uniformly in $\left(\mathbf{x},t\right)$ in the domain of interest. I.e., the problems (2.1) and (2.2) are equivalent up to $O\left(\varepsilon\right)$.

We call the map $\mathbf{g}\longrightarrow$ $\widetilde{\mathbf{g}}$ the* data projection * of the manifold $\Gamma$ onto the hyperplane $R^{n}\times\left\{t=0\right\}$. Finding $\widetilde{\mathbf{g}}$ is based upon the Taylor formula in ”reverse”.

Apply the Taylor formula in one variable $t$ to the (still unknown) solution $\mathbf{u}\left(\mathbf{x},t\right)$ of (2.1). For each fixed point $(\mathbf{x},t)$ we then have

[TABLE]

with some error $E_{j}$. I.e., (2.3) is the Taylor formula about $\left(\mathbf{x},t\right)$ evaluated at $\left(\mathbf{x},0\right)$ (not the other way around). Taking in (2.3) $t=\tau\left(\mathbf{x}\right)$ yields

[TABLE]

We call $\mathbf{g}_{j}$ the $j^{\text{th}}$ order projection of initial data $\mathbf{g}$ onto $R^{n}\times\left\{t=0\right\}$. One can see now that if we are able to properly compute all $\left.\dfrac{\partial^{k}\mathbf{u}(\mathbf{x},t)}{\partial t^{k}}\right|_{\Gamma}$ then $\mathbf{g}_{j}$ produces a desirable standard IC $\widetilde{\mathbf{g}}$ . Indeed, given error $\varepsilon$ (no matter how small), we take $j$ so large as $\left|E_{j}\left(\mathbf{x},t\right)\right|<\varepsilon$ uniformly in the domain of interest and hence

[TABLE]

Thus the solution $\widetilde{\mathbf{u}}$ to (2.2) with $\widetilde{\mathbf{g}}=\mathbf{g}_{j}$ will coincide with the solution $\mathbf{u}$ of the original problem (2.1) up to $O\left(\varepsilon\right)$.

We now show how to compute the Taylor coefficients in (2.4). The zeroth one is obvious

[TABLE]

All of the other Taylor coefficients in (2.4) can also be explicitly computed. Start with the first-order one:

[TABLE]

We now compute $\left.(D\mathbf{u})\right|_{\Gamma}$ by the chain rule. One has

[TABLE]

It follows that

[TABLE]

and thus for the first order Taylor coefficient we finally have

[TABLE]

where, recalling (1.4),

[TABLE]

Our computation of higher order Taylor coefficients will be based on the following observation. Since $\partial/\partial t$ and $D$ commute, $\mathbf{u}_{t}$ is also a solution to

[TABLE]

with the initial condition given by (2.6)

[TABLE]

Thus, if $\mathbf{u}$ is the solution originated from $\mathbf{g}$ then $\mathbf{u}_{t}$ is the solution originated from the initial data $-\left(I-A\right)^{-1}D\mathbf{g}$. By induction one concludes that $\dfrac{\partial^{k}\mathbf{u}}{\partial t^{k}}$ is the solution originated from $\left[-\left(I-A\right)^{-1}D\right]^{k}\mathbf{g}$, $k=2,3,4,...$

Therefore, we get the following nice formula

[TABLE]

Substituting this into (2.4) one has

[TABLE]

and we finally arrive at

Theorem 2.1**.**

Let $\Gamma=\{(\mathbf{x},\tau(\mathbf{x}))|\mathbf{x}\in$ $\mathbb{R}^{n}\}$ be a non-characteristic manifold for the hyperbolic system

[TABLE]

and

[TABLE]

Fix accuracy $\varepsilon$ and let $j$ be chosen to satisfy ($\left\|\cdot\right\|$ is the Euclidean norm in $\mathbb{R}^{m}$)

[TABLE]

[TABLE]

If $\widetilde{\mathbf{u}}$ is the solution to (2.8) subject to the initial condition

[TABLE]

then the solution $\mathbf{u}$ to (2.8)-(2.9) is subject to $\mathbf{u}\left(\mathbf{x},t\right)=\widetilde{\mathbf{u}}\left(\mathbf{x},t\right)+O\left(\varepsilon\right).$

The map $\mathbf{g}\longrightarrow\mathbf{g}_{j}$ is linear and well-defined as long as the matrix $I-A$ is non-singular, which holds if (1.5) does. An important feature of our data projection method is that, by construction, both $\mathbf{u}$ and $\widetilde{\mathbf{u}}$ solve (exactly) the same system (1.1) but satisfy different (equivalent) IC conditions $\mathbf{g}$ and $\widetilde{\mathbf{g}}=\mathbf{g}_{j}$. Of course $\mathbf{u}$ and $\widetilde{\mathbf{u}}$ can be made as close as one wishes (while $\mathbf{g}$ and $\widetilde{\mathbf{g}}$ need not be close at all).

3. Applications to the non-linear shallow water wave system

In this section we apply our formalism to the hyperbolic 1+1 quasi-linear shallow water system (see, e.g. [12]) modeling e.g. the tsunami wave run-up and run-down. For the so-called inclined bathymetries this system reads (in dimensionless variables)

[TABLE]

where $\eta$ is the unknown water elevation over unperturbed level, $u$ is the unknown cross-section averaged flow velocity, and $c\geq 0$ is a known function solely encoding the information about the shape of our inclined bathymetry [16], [15]. The point $x=0$ in (3.1) corresponds to the unperturbed coast line and the $x$-axis is directed off shore. The main feature of this problem is a moving (wet/dry) boundary also know as run-up/run-down. The substitution

[TABLE]

referred to as the Carrier-Greenspan (CG) transform or CG hodograph, turns (3.1) into the linear (strictly) hyperbolic system

[TABLE]

where $\Gamma$ and $\mathbf{\phi}_{0}$ will be given later. The CG transform was introduced in [3] and has become a standard tool in the study of the run-up/run-down process (see, e.g. [6], [16], [18] and the extensive literature cited therein). The form (3.2) is taken from our [15]. Besides linearizing (3.1), the CG turns the moving boundary into the fixed point $\sigma=0$. It has however a serious drawback: the IC in (3.5) is no longer standard. Indeed, under (3.2), the (horizontal) line $t=0$ in the plane $\left(x,t\right)$ becomes the parametric curve $\Gamma=\left(x+\eta_{0}(x),-u_{0}(x)\right)$ in the $(\sigma,\tau)$ plane and, denoting the inverse of $x+\eta_{0}(x)$ by $\gamma\left(\sigma\right)$, one has

[TABLE]

which defines $\mathbf{\phi}_{0}$ in (3.5). We see from (3.13) that the IC (3.20) is standard iff the initial velocity $u_{0}=0.$ The latter has been a typical assumption in much of the previous literature (which however does not give a valid inundation picture caused by a tsunami wave). That is why it has been a good open problem since [3] how to make the CG transform run for general IC. We refer to [4], [10], [11], [2] where this problem was addressed under certain assumptions of relative smallness of the IC (and only in the context of the plane beach). We discovered our method in [13] in a particular case while trying to put previous works on a solid footing.

We now apply the method of data projection to (3.5). Given accuracy $\varepsilon$, choose $j$ so that

[TABLE]

and construct the $j$th order projection

[TABLE]

We can then solve (3.5) with the standard IC $\left.\mathbf{\phi}\right|_{\tau=0}=\mathbf{\phi}_{j}\left(\sigma\right)$ by any suitable method. In fact, for the so-called power shaped bays, $c\left(\sigma\right)\sim\sigma$ and it can be solved by the Hankel transform in terms of Bessel functions [13]. Performing the inverse CG transform (3.2) solves the original problem (3.1) in the physical space. The latter is, in general, not explicit but can easily be done numerically without affecting the total accuracy, which remains $O\left(\varepsilon\right).$ In fact, we can call our method exact as the error it introduces can be made negligible comparing with the one inherited by the shallow water approximation leading to the very system (3.1).

Note that while (3.21) looks unwieldy (the reader is invited to amuse him/herself with trying to unzip it even for $j=1$), its numerical implementation is not a problem. It was the matrix form (3.5) that made our derivation quite transparent.

We emphasize that an important feature of our method is that the procedure works smoothly as long as the wave non-breaking condition [16] is satisfied and no extra assumptions of smallness (made in the previous literature) are needed. Extensive numerical verification and simulations (which will be published elsewhere) show that our method is robust and can be effectively used for rapid forecasting of characteristics of inundation zone due to large-amplitude sea waves.

In the conclusion, we mention that our method can be adapted to treat boundary conditions and improve on the relevant results of [1], [2], [9], [18]. We plan to address this elsewhere.

4. Acknowledgments

We are grateful to Efim Pelinovsky for stimulating email discussions and an annonemous referee for useful literature hints.