Analysis and numerical simulation of the three-dimensional Cauchy problem for quasi-linear elliptic equations

Accepted Manuscript Analysis and numerical simulation of the three-dimensional Cauchy problem for quasi-linear elliptic equations

Quoc Viet Tran, Mokhtar Kirane, Nguyen Huy Tuan, Van Thinh Nguyen

PII: DOI: Reference:

S0022-247X(16)30477-2 http://dx.doi.org/10.1016/j.jmaa.2016.08.045 YJMAA 20681

To appear in:

Journal of Mathematical Analysis and Applications

Received date:

18 January 2016

Please cite this article in press as: Q.V. Tran et al., Analysis and numerical simulation of the three-dimensional Cauchy problem for quasi-linear elliptic equations, J. Math. Anal. Appl. (2016), http://dx.doi.org/10.1016/j.jmaa.2016.08.045

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Analysis and numerical simulation of the three-dimensional Cauchy problem for quasi-linear elliptic equations Quoc Viet Tran1 , Mokhtar Kirane

2,3

, Nguyen Huy Tuan

4 ˚

, Van Thinh Nguyen

1

1

Department of Civil and Environmental Engineering, Seoul National University, Republic of Korea. 2 LaSIE, Facult´e des Sciences et Technologies , Universi´e de La Rochelle, Avenue M. Cr´epeau, 17042 La Rochelle Cedex, France 3 Nonlinear Analysis and Applied Mathematics (NAAM) Research Group, Department of Mathematics, Faculty of Science, King Abdulaziz University, P.O. Box 80203, Jeddah 21589, Saudi Arabia 4 Applied Analysis Research Group, Faculty of Mathematics and Statistics, Ton Duc Thang University, Ho Chi Minh City, Viet Nam

Abstract This work is concerned with solving the Cauchy problem for quasilinear elliptic equations whose exponential instability is manifestly seen by the catastrophic growth in the representation of the exact solution. Our proposed regularization procedure consists in damping the unbounded terms in the representation. Moreover, we show that its solution converges to the exact solution uniformly and strongly in L2 under a priori assumptions on the exact solution. In order to verify our analysis and the accuracy of the numerical procedures, we exhibit two numerical examples. Our main tools for simulation are the trigonometric polynomial approximation, and the fast Fourier transform in combination with the cubic Hermite interpolation. Keywords and phrases: Elliptic equation, Ill-posed problem, Cauchy problem, Contraction principle, Regularization method, Fast Fourier Transform.

1

Introduction

It is well known that the Cauchy problem for second-order elliptic partial differential equations is ill-posed in the sense of Hadamard [4], i.e., the solution does not exist, or even it exists, it does not depend continuously on the Cauchy data. An answer to the ill-posedness issue of such Cauchy problem can be found in Faker Bin Belgacem et al. [1]. Even though the ill-posedness causes difficulty for numerical computation, such problems usually appear in many engineering applications related to propagating waves in different environments, such as acoustic, hydrodynamic and electromagnetic waves [5, 6]. Besides, most of those problems are usually solved in 3D domains with inhomogeneous sources, and particularly, source functions depend on the unknown function. Hence a necessity of investigating Cauchy problems for nonlinear elliptic equations motivates this study. Let a, b, c ą 0 and Ω “ p´a, aq ˆ p´b, bq be a rectangle in R2 with the boundary BΩ. We seek a ˚

Corresponding author: [email protected]

1

function u that satisfies: Δu “ f pu, x, y, zq, upx, y, zq “ 0, ε

¯ ˆ r0, cs, px, y, zq P Ω

(1)

px, y, zq P BΩ ˆ r0, cs,

(2)

ε

}g ´ up¨, ¨, 0q} ` }h ´ Bz up¨, ¨, 0q} ď ε,

(3)

where Δ is the three dimensional Laplacian, Bz is the z-partial derivative, f is a given function depending on the unknown u, both g ε and hε are given functions in the Lebesgue space L2 pΩq that is endowed with the natural inner product ¨, ¨ and with the norm }¨}, the positive number ε is the data error of pg ε , hε q from the exact Cauchy data ¯ px, yq P Ω.

pg, hq “ pu, Bz uq|z“0 ,

(4)

Throughout this paper, we shall use the following notations. The space H m pΩq is the Sobolev space involving all functions that are in L2 pΩq as well as their s-th order derivatives for all s ď m, H01 pΩq contains all the functions of H 1 pΩq such that their traces vanish on BΩ. We shall adopt the notation pCpr0, cs; L2 pΩqq, |||¨|||q for the Banach space of continuous functions mapping r0, cs to L2 pΩq, wherein |||¨||| stands for the supremum norm. Let λ2mn and ψmn denote the eigenvalues and corresponding eigenfunctions of A :“ ´Δ defined on its domain DpAq Ă H01 pΩq, which are c ˙ ˆ ˙ ˆ nπpy ` bq mπpx ` aq π ´ m ¯2 ´ n ¯2 sin , (5) ` , ψmn px, yq “ sin λmn “ 2 a b 2a 2b for any pm, nq P N2 . Here after, we shall denote the Fourier coefficients of the functions v “ vpx, yq, ˆmn pzq “ κ wp¨, ¨, zq, ψmn  and fˆmn pw, zq “ w “ wpx, y, zq and f “ f pw, x, y, zq by vˆmn “ κ v, ψmn , w ´2 ˆmn pzq will stand κ f pwp¨, ¨, zq, ¨, ¨, zq, ψmn , respectively, where κ “ }ψmn } “ 1{pabq. Similarly, Bz w for κ Bz wp¨, ¨, zq, ψmn . Using the method of separation of variables, the exact solution of the problem (1)-(4) satisfies upx, y, zq “

8 ´ 8 ÿ ÿ

¯ p mn pg, h, zq ` Jpmn pu, zq ψmn px, yq, G

(6)

m“1 n“1

where zλmn p mn pg, h, zq “ e G 2

Jpmn pu, zq “

1 2λmn

˜

˜ ¸ ¸ p p hmn hmn e´zλmn gpmn ` gpmn ´ ` , λmn 2 λmn żz´ ¯ epz´sqλmn ´ eps´zqλmn fpmn pu, sq ds.

(7) (8)

0

ˆ mn pg, h, zq and Jˆmn pu, zq in Eqs. (7) and (8) blow up quickly as λmm tends to It can be seen that G infinity due to the fast increasing of exppzλmn q. Therefore, a numerical calculation of Eqs. (6)-(8) ˆ mn , fˆmn q may tend to zero in practice is impossible, even when the exact Fourier coefficients pˆ gmn , h rapidly. It is well-known that the Cauchy problem for elliptic equations is severely ill-posed in the sense of Hadamard, i.e, a small pertubation in the given Cauchy data may cause a very large error in the output solution upx, y, zq for z P r0, cs. The instability increases with the increase of the distance z from the boundary z “ 0. Therefore, it is very difficult to solve the problem by using classical numerical methods of inversion. In order to overcome this instability, regularization methods are naturally required. In the past, there were many studies on the Cauchy problem for specific forms of the elliptic equation (1). In the absence of the source function, i.e. f “ 0, the problem is known as the Cauchy problem for the Laplace equation, see [2, 3]. For f “ ´k 2 u (k is a constant), (1) reduces to the homogeneous Helmholtz equation; it has been extensively studied and many results concerning regularization methods have been investigated, see, e.g. the recent work by Reginska and her group [7, 8].

2

Recently, Nguyen et al. [12] considered Eq. (1) in 2D for a modified Helmholtz equation (or the Yukawa equation) with an inhomogeneous source, i.e. the source function is the form f “ k 2 u ` r (r is a function). Later, Tran et al. [13] extended the results of his group ([12]) to solve the problem in 3D for Helmholtz-type equation, i.e. Eq. (1) with a source function f “ ˘k 2 u ` r associated with both homogeneous Dirichlet and Neumann boundary condition. The aforementioned equations were also generalized to abstract settings, for example, [9, 10, 11] proposed many regularization strategies to the operator equation. In [9], Elden et al. applied the cutoff method to obtain the stabilized solution, and treated numerical procedures by the rational Krylov method. In [10], the authors proposed a modified method by constructing new bounded kernels in which they substitute the unbounded terms in the representation of the solution, and obtained various error estimates corresponding to several a priori conditions on the exact solution. To our knowledge, there are very few results on the Cauchy problem for nonlinear elliptic equations in three dimensions. In this paper, we consider problem (1)-(4) in the case of a function f satisfying the Lipschitz condition DK ą 0, @w, v P L2 pΩq, @z P R, }f pw, zq ´ f pv, zq} ď K }w ´ v} .

(9)

This paper is organized as follows. To regularize a solution of problem (1)-(4), a regularization problem based on damping the unbounded terms exppzλmn q in Eqs. (7) and (8) will be introduced at the beginning of Section 2. In this section, we shall prove that the regularized problem is well-posed (Theorem 1) and its solution converges to the exact solution uniformly under certain assumptions (Theorem 2). Section 3 is reserved to numerical experiments, which is divided into three subsections. In section 3.1, we model the Cauchy data from its discrete data set by using the trigonometric polynomial approximation [14]. This application allows us to use the fast Fourier transform (FFT). In section 3.2, we explain calculation procedures combining FFT with the cubic Hermite interpolation [15] to solve the regularized problem numerically. In section 3.2, two test cases with numerical results and discussions will be presented to verify the stability and convergence of the regularized solution as well as the accuracy of the numerical procedures. Finally, the conclusion is presented in Section 4.

2

Theoretical Results

Regularization methods. Let T ě c be a constant. For each α ą 0, we seek for a function uε,α that satisfies 8 ´ 8 ÿ ¯ ÿ α p αmn pg, h, zq ` Jpmn puα , zq ψmn px, yq, G uα px, y, zq “

(10)

m“1 n“1

where ˜ ¸ ¸ ˜ p p hmn hmn e´zλmn e´pT ´zqλmn gpmn ´ ` , gpmn ` λmn 2 λmn 2 pαλmn ` e´T λmn q ¸ ż z ˜ ´pT ´z`sqλmn e 1 α Jpmn puα , zq “ ´ eps´zqλmn fpmn puα , sq ds, 2λmn 0 αλmn ` e´T λmn

p α pg, h, zq “ G mn

(11) (12)

for px, y, zq P Ω ˆ r0, T s. Here, α ą 0 is a regularization parameter depending on . The main theoretical results rely on the following lemmas. Lemma 1. Let p ě 0 and q ą 0. We have ´ ´ q ¯¯ p ´1 e´pλ q 1´ pq ďq , α ln α αλ ` e´qλ for any λ ą 0 and α P p0, qq. Moreover, if q ě p, one has ˆ ˆ ˙˙ p ´1 q D e´pλ ď B α ln , ´qλ α αλ ` e for any λ ą 0 and α P p0, Dq, where B “ maxt1, qu and D “ mint1, qu.

3

(13)

(14)

Proof. The proof of Lemma 1 can be found in [16]. Lemma 2. Let φ, ϕ, σ, ς P L2 pΩq and w, v P Cpr0, T s; L2 pΩqq arbitrarily. Then the following estimates hold ˆ ˆ ˙˙´ 2z ˆˇ ˙ ˇ ˇ2 ˇ2 T D ˇ pα ˇ ˇ ˇ 2 α p p pmn ´ φmn ˇ ` |p (15) ςmn ´ σ pmn | , ˇGmn pϕ, ς, zq ´ Gmn pφ, σ, zqˇ ď C1 α ln ˇϕ α and ˇ2 ˇ ˇ ˇ pα α pv, zqˇ ď ˇJmn pw, zq ´ Jpmn

z λ2mn

żz˜ 0

¸

e´2pT ´z`sqλmn 2

pαλmn ` e´T λmn q

` e2ps´zqλmn

ˇ ˇ2 ˇp ˇ ˇfmn pw, sq ´ fpmn pυ, sqˇ ds, (16)

(  2 2 ´2 where C1 “ 2B 2 max 1, λ´2 min and C2 “ 2B K λmin , λmin is the smallest value of λmn defined in (5), B and D are defined in Lemma 1, and K is the Lipschitz constant in (9). Proof. Substituting pϕ, ςq and pφ, σq into (11), we obtain ¸ ¸˜ ˜ ϕ pmn ´ φˆmn e´pT ´zqλmn ´zλmn α α p p `e Gmn pϕ, ς, zq ´ Gmn pφ, σ, zq “ 2 αλmn ` e´T λmn ˜ ¸ˆ ˙ e´pT ´zqλmn ςpmn ´ σ pmn ´zλmn ` . ´e 2λmn αλmn ` e´T λmn Using Cauchy-Schwarz inequality, we get ¸˜ ¸ ˜ ˇ ˇ2 |p ˇ ˇ2 ςmn ´ σ pmn |2 e´2pT ´zqλmn ˇ pα ˇ ˇ ˇ α ´2zλmn p p pmn ´ φmn ˇ ` , ˇGmn pϕ, ς, zq ´ Gmn pφ, σ, zqˇ ď ˇϕ 2 `e λ2mn pαλmn ` e´T λmn q in which we use (14), with p :“ T ´ z and q “ T gives ˆ ˆ ˙˙´ 2z T D e´2pT ´zqλmn ´2zλmn 2 α ln ` e ď 2B , @α P p0, Dq, 2 α pαλmn ` e´T λmn q ` ` ˘˘´ 2z T . This implies (15). By substituting w and v into (12) and using since e´2zλmn ď 1 ď α ln D α Cauchy-Schwarz and H¨older’s inequalities, we get (16). 2 Theorem 1 (Well-posedness of problem (10)-(12)). For any α P p0, ` Dq and 2φ, σ ˘P L pΩq, the problem α (10)-(12) associated with data pφ, σq has a unique solution w P C r0, T s; L pΩq . Furthermore, if v α is the solution corresponding to the data pϕ, ςq P L2 pΩq ˆ L2 pΩq, then one has ˆ ˆ ˙˙´ z T D α α p}ϕ ´ φ} ` }ς ´ σ}q , (17) }w pzq ´ v pzq} ď Q1 α ln α

where Q1 is a positive constant; the problem (10)-(12) is well-posed. Proof. The proof of the existence of a unique solution is based on references [16, 17]. Let us introduce a function F as follows 8 ´ 8 ÿ ¯ ÿ α p αmn pφ, σ, zq ` Jpmn F rwspzq “ pw, zq ψmn , (18) G m“1 n“1

˘ p α and Jpα are defined by (11) and (12). Using the for w P C r0, T s; L2 pΩq and z P r0, T s, where G mn mn induction principle, we shall prove that d › › pQzq2l › l › @l P N, ›F rwspzq ´ F l rvspzq› ď |||w ´ v||| , (19) 1.3.5...p2l ´ 1q `

4

for any w, v P Cpr0, T s; L2 pΩqq, where Q is independent of l. Using the estimate (16), we get żz 2 2 α α 2 }F rwspzq ´ F rvspzq} “ }J pw, zq ´ J pv, zq} ď Q z }wpsq ´ vpsq}2 ds ď pQzq2 |||w ´ v|||2 , ?

0

` D ˘˘

and C2 is the constant in Lemma 2. Thus (19) holds for l “ 1. Now where Q “ C2 { α ln α assume that (19) holds for l “ k, we shall prove that (19) holds for l “ k ` 1. Indeed, we have › ›2 › ” ı ” ı ›2 › k`1 › › › rwspzq ´ F k`1 rvspzq› “ ›F F k rws pzq ´ F F k rvs pzq› ›F żz› ›2 › k › 2 ďQ z ›F rwspsq ´ F k rvspsq› ds `

0

ď Q2 z

żz 0

ď

pQsq2k |||w ´ v|||2 ds 1.3.5...p2k ´ 1q

pQzq2pk`1q |||w ´ v|||2 . 1.3.5...p2k ´ 1qp2k ` 1q

Thus (19) holds for l “ k ` 1. Hence using the induction principle, we can state that the estimate (19) hold for any l P N. Whereupon d ˇˇˇ ˇˇˇ pQT q2l ˇˇˇ ˇˇˇ l l |||w ´ v||| , @w, v P Cpr0, T s; L2 pΩqq. (20) ˇˇˇF rws ´ F rvsˇˇˇ ď 1.3.5...p2l ´ 1q c c pQT q2l 1.3.5...p2l´1q

Since limlÑ8

“ 0, there exists an integer l0 ą 0 such that

pQT q2l0 1.3.5...p2l0 ´1q

ă 1, and thus,

F l0

the operator is a contraction. Applying the Banach fixed-point theorem, we conclude that equation F l0 rws “ w has a unique solution in Cpr0, T s; L2 pΩqq, say wα , which is the fixed point of the operator F l0 . Then, we have F pF l0 pwα qq “ F pwα q, or, F l0 pF pwα qq “ F pwα q. According to the uniqueness of the fixed point, we obtain F pwα q “ wα . In other words, the problem (10)-(12) has a unique solution wα P Cpr0, T s; L2 pΩqq. In order to prove the stability condition (17), we substitute pϕ, ςq and pφ, σq into (10), then take the inner product with ψmn ; we obtain α p αmn pφ, σ, zq ` Jpmn wα pzq, ψmn  “ G pwα , zq, p α pϕ, ς, zq ` Jpα pv α , zq. v α pzq, ψmn  “ G mn

mn

We clearly have }wα pzq ´ v α pzq}2 ď 2

8 ˇ 8 ˇ 8 ÿ 8 ÿ ˇ2 ˇ2 ÿ ÿ ˇ ˇ pα ˇ pα α p αmn pφ, σ, zqˇˇ ` 2 pv α , zqˇ . ˇGmn pϕ, ς, zq ´ G ˇJmn pwα , zq ´ Jpmn

m“1 n“1

m“1 n“1

First, using Lemma 2 and Parseval’s identity for the latter estimate, we conclude that ˆ ˆ ˙˙´ 2z 8 ˇ 8 ÿ ˇ2 ÿ T D ˇ ˇ pα α p p}ϕ ´ φ} ` }ς ´ σ}q2 . ˇGmn pϕ, ς, zq ´ Gmn pφ, σ, zqˇ ď C1 α ln α m“1 n“1

(21)

Using Lemma 1, we have e´2pT ´z`sqλmn pαλmn `

2 e´T λmn q

`e

2ps´zqλmn

ď 2B

2

ˆ

α ln

ˆ

D α

˙˙´ 2pz´sq T

, @α P p0, Dq.

(22)

It follows from (22) that ˆ ˙˙´ 2pz´sq 8 ˇ 8 ÿ ˇ2 2B 2 z ż z ˆ ÿ T D ˇ ˇ pα α α α p }f pwpsq, sq ´ f pvpsq, sq}2 ds. α ln ˇJmn pw , zq ´ Jmn pv , zqˇ ď 2 α λ 0 min m“1 n“1 (23)

5

Combining (21) and (23) leads to α

α

ˆ

2

}w pzq ´ v pzq} ď C3 α ln ` C4

D α

ˆ

˙˙´ 2z

żzˆ

T

α ln

ˆ

0

D α

p}ϕ ´ φ} ` }ς ´ σ}q2

˙˙´ 2pz´sq T

}wα psq ´ v α psq}2 ds,

(24)

where C3 “ 2C1 and C4 “ 2C2 c, here C1 and C2 are the constants in Lemma 2. The inequality (24) can be written in the following form ˆ

α ln

ˆ

D α

˙˙ 2z T

α

α

2

2

}w pzq ´ v pzq} ď C3 p}ϕ ´ φ} ` }ς ´ σ}q ` C4

żzˆ

α ln

0

ˆ

D α

˙˙ 2s T

}wα psq ´ v α psq}2 ds.

Gronwall’s inequality leads to ˆ ˆ ˙˙ 2z T D α ln }wα pzq ´ v α pzq}2 ď C3 eC4 z p}ϕ ´ φ} ` }ς ´ σ}q2 . α

(25)

ˆ ˆ ˙˙´ z T D }wα pzq ´ v α pzq} ď Q1 α ln p}ϕ ´ φ} ` }ς ´ σ}q , α

(26)

Therefore

where Q1 “

a

C3 exppC4 T q. This completes the proof.

Note that the well-posed problem (10)-(12) was defined for z P r0, T s, where T ě c. Particularly, in case T ą c, the solution wα in Theorem 1 is quite smooth on p0, cs. Indeed, taking derivative of (10) with respect to z, one can check that 8 8 ÿ ÿ m“1 m“1

α |Bz wmn pzq|2 ă 8 and

8 ˇ 8 ÿ ÿ ˇ2 α ˇλ2mn wmn pzqˇ ă 8,

m“1 m“1

for any z P p0, T q, and thus, Bz wα pzq P L2 pΩq and wα pzq P H 2 pΩq X H01 pΩq. The following theorem shows that restriction of wα on r0, cs converges uniformly to the analytical solution of problem (1)-(4) under some a priori assumptions on the solution u. Theorem 2 (Regularization strategy for the problem (1)-(4)). Assume that the problem (1)-(4) admits a unique exact solution u P Cpr0, cs; H01 pΩqq X C 1 pr0, cs; L2 pΩqq. For each ε P p0, Dq and pg ε , hε q satisfying (3), a regularized function uε,α P Cpr0, cs; L2 pΩqq can be constructed from the problem (10)(12) associated with the Cauchy data pg ε , hε q. Choosing αpεq such that ε bounded, εÑ0 αpεq

lim αpεq “ 0, and lim

εÑ0

(27)

we morever have (1) If u satisfies 8 żc 8 ÿ ÿ m“1 n“1 0

e2pc´sqλmn |fˆmn pu, sq|2 ds ă 8,

(28)

then one has ε,α

|||u

ˆ

´ u||| ď Q1 ε α ln

ˆ

D α

˙˙´ z

T

6

` P1 α

c´z T

ˆ ˆ ˙˙´ T `z´c T D ln . α

(29)

(2) If u satisfies Dγ ą 0,

sup

8 8 ÿ ÿ

zPr0,cs m“1 n“1

e2pγ`c´zqλmn |λmn u ˆmn pzq ` Bz u ˆmn pzq|2 ă 8,

(30)

then setting up the problem (10)-(12) with T ě c ` γ, we obtain |||u

ε,α

ˆ

´ u||| ď Q1 ε α ln

ˆ

D α

˙˙´ z

T

` P2 α

γ`c´z T

ˆ

ˆ ln

D α

˙˙ γ`c´z ´1 T

,

(31)

here Q1 is defined in Theorem 1, P1 is defined in (41), P2 is defined in (44), D “ mint1, cu and |||¨||| is the supremum norm of Cpr0, cs; L2 pΩqq. Remark 1. 1. Using (27) and 0 ď z ď c ď T , it is easy to check that the right hand sides of (29) and (31) tend to zero as ε Ñ 0. Hence, our method is a regularization method. Moreover, we give a rule for choosing αpεq such that (27) holds. We can choose α “ εa , 0 ă a ď 1 or α “ εa

´ a ¯b 0 , 0 ă a ă 1, b ě 0, a0 ě 0, a1 ą 0. lnp aε1 q

2. If we choose α “ ε, then (27) holds. In this case, the estimate (29) becomes ˜ˆ ¸ˆ ˆ ˙˙1´ c ˆ ˙˙ c p1´ z q T T c D D maxtQ1 , P1 u ε,α `D˘ . ε ln `1 ε ln }u pzq ´ upzq} ď ε ε ln ε

(32)

For T ě c, we have ε lnpD{εq ă 1 for any ε P p0, Dq due to (14). This implies that }uε,α pzq ´ upzq} is of order ln 1D . pεq 3. If we also choose α “ ε then the estimate (31) becomes ¨ ˛ ˆ ˆ ˙˙ c´z γ T T T ´pc`γq D ε P 2 ε,α }u pzq ´ upzq} ď ` ` ˘˘ c ˝Q1 ε T ` ` ` ˘˘ T ´pc`γq ‚ ε ln . ε T ln Dε T ln Dε Setting M “ maxtQ1 , P2 u and using the inequality (14), i.e. ε lnpD{εq ă 1 for any ε P p0, Dq, the latter estimate yields ¨ ˛ ˆ ˆ ˙˙ pc`γq´T γ T T T ´pc`γq D M ε ‚, @z P r0, cs. }uε pzq ´ upzq} ď ` ` ˘˘ c ˝ε T ` ln (33) ε ln Dε T Remark 2. The assumptions (28) and (30) are meaningful in the spirit of abstract Gevrey spaces (but obviously rather strict than the Gevrey regularity). Indeed, if the exact solution is supposed to be a real analytic function in r0, cs with values in Gevrey spaces, particularly ´ ¯¯ ´ ´ ¯¯ ´ 1{2 1{2 X C 1 r0, cs ; D eT p´Aq C r0, cs ; H01 pΩq X D p´Aq1{2 eT p´Aq ¯ ´ ¯ ´ 1{2 1{2 and D p´Aq1{2 eT p´Aq with T ě c`γ in which we recall the norms of the Gevrey classes D eT p´Aq by the following: }U }D´eT p´Aq1{2 ¯

› › 1{2 › › “ ›eT p´Aq U › “

˜

8 8 ÿ ÿ m“1 n“1

ˇ ˇ ˇ ˆ ˇ2 e2T λmn ˇU mn ˇ

7

¸1{2 ,

for

´ ¯ 1{2 U P D eT p´Aq ,

˜ }U }D

1{2 p´Aq1{2 eT p´Aq

´

¯

“

8 8 ÿ ÿ m“1 n“1

λ2mn e2T λmn

ˇ ˇ ˇ ˆ ˇ2 ˇUmn ˇ

¸1{2 ,

for

´ ¯ 1{2 U P D p´Aq1{2 eT p´Aq ,

then the assumption (30) is satisfied. Notice that ´ ¯¯ ´ ´ ¯¯ ´ 1{2 1{2 X C 1 r0, cs ; D eT p´Aq C r0, cs ; H01 pΩq X D p´Aq1{2 eT p´Aq ˘ ˘ ` ` is a closed subspace of C r0, cs ; H01 pΩq X C 1 r0, cs ; L2 pΩq . In´ the same´ vein, (28) ¯¯holds when the exact solution is such that the source function f belongs 1{2 T p´Aq to C r0, cs ; D e for T ě c. On the other hand, if f p0, x, y, zq “ 0 for all px, y, zq P ¯ Ω´ ˆ r0, cs, we need ´ only 1{2 ¯¯ to make a reasonable Gevrey assumption on the exact solution, that is u P T p´Aq C r0, cs ; D e , because ˜

8 8 ÿ ÿ m“1 n“1

e

2T λmn

ˇ2 ˇ ˇ ˇˆ puq f ˇ ˇ mn

¸1{2 “ }f puq}D´eT p´Aq1{2 ¯ ď K }u}D´eT p´Aq1{2 ¯

by the global Lipschitz continuity of f . Remark 3. In practice, it is impossible to compute the regularized solution uα by the infinite series as the formula (10)-(12). To approximate uε,α , we need a truncation solution given by uε,α,M px, y, zq “

λmn ďM ÿ m,n“1

¯ α p αmn pg ε , hε , zq ` Jpmn puε,α,M , zq ψmn px, yq. G

´

(34)

The level of error }uε,α,M pzq ´ upzq} should be proved. This is a task for a future work. Proof. The outline of the proof is described as follows. Let uα and uε,α be the solutions of the problem (10)-(12) associated with the Cauchy data pg, hq and pg ε , hε q, respectively. Firstly, under the given assumptions, we show that uα converges to u uniformly as α tends to zero. Secondly, uα approximates uε,α via a choice of αpεq based on the stability of the well-posed problem (10)-(12), which property was featured by Theorem 1. Finally, we obtain the convergence estimate of uε in Cpr0, cs; L2 pΩqq; as a consequence, the solution u is unique. Proof of (1).

We shall estimate }uα pzq ´ upzq} as follows. Using (6)-(8) and (10)-(12), we have upzq ´ uα pzq “

8 ´ 8 ÿ ÿ m“1 n“1

¯ α α pu, zq ´ Jpmn puα , zq ψmn , Hmn ` Jpmn

(35)

where α p mn pg, h, zq ´ G p αmn pg, h, zq ` Jpmn pu, zq ´ Jpmn pu, zq, z P r0, cs. Hmn “ G

(36)

Subtracting (11) from (7) and (12) from (8), respectively, then grouping the result with a care, (36) can be rewritten as follows ¸ ˆ ˙˜ ż z ´sλmn p hmn e e´T λmn 1 1´ gpmn ` ` (37) fpmn pu, sq ds ezλmn . Hmn “ 2 λmn λmn αλmn ` e´T λmn 0 On the other hand, differentiating (7) and (8) with respect to z (see also Lemma 2 [12]), we deduce that ¸ ˜ ż z ´sλmn p hmn e ˆmn pzq Bz u ` pmn pzq ` . (38) fpmn pu, sq ds ezλmn “ u gpmn ` λmn λ λmn mn 0

8

Likewise, substituting z :“ c into (38), a simple computation gives ˆ ˙ żc α e´pc´zqλmn pc´sqλmn ˆ Hmn “ pmn pcq ` Bz u pmn pcq ´ λmn u e fmn pu, sq ds . 2 αλmn ` e´T λmn z Now, applying Lemma 1 with p :“ c ´ z and q :“ T , we get ˜ ¸2 ˇ ˇ2 żc ˇ ˇ e´pc´zqλmn α2 2 pc´sqλmn ˆ ˇλmn u ˇ |Hmn | ď ˆ pcq ` B u ˆ pcq ´ e pu, sq ds f mn z mn mn ˇ ˇ ´T λ mn 4 αλmn ` e z ˆ ˙ żc α2 p3B 2 {4q 2 2 2pc´sqλmn ˆ 2 |λ ď` u ˆ pcq| ` |B u ˆ pcq| ` c e | f pu, sq| ds . mn mn z mn mn ` ˘˘ 2c p T ´1` z q T c c z α ln D α Therefore, we have the following estimate 8 8 ÿ ÿ m“1 n“1

|Hmn |2 ď `

α ln

M1 α 2 ` D ˘˘ 2c p T ´1` z q , T

α

c

c

where M1 “ 3B 2 E1 maxtκ, cu{4, E1 “ }Bx upcq}2 ` }By upcq}2 ` }Bz upcq}2 `

8 żc 8 ÿ ÿ m“1 n“1 0

e2pc´sqλmn |fˆmn pu, sq|2 ds,

(39)

due to the assumption (28). It should be mentioned that we used the following estimate 8 8 ÿ ÿ

´ ¯ |λmn u ˆmn pcq|2 ď κ }Bx upcq}2 ` }By upcq}2 ,

(40)

m“1 n“1

which can be obtained from Parseval’s theorem and Green’s first identity, as upzq P H01 pΩq for any z P r0, cs. Now, applying Parseval’s identity to (35) and using the estimate (16), we deduce that 8 8 2 ÿ ÿ |Hmn |2 ` 2 }J α pu, zq ´ J α puα , zq}2 κ m“1 n“1 ˆ ˙˙´ 2pz´sq żzˆ T D 2M1 κ´1 α2 ` 2C α ln ď` c }upsq ´ uα psq}2 ds. 2 2c T z ` D ˘˘ p ´1` q α T c c 0 α ln α

}upzq ´ uα pzq}2 ď

Thus, ˆ ˆ ˙˙ 2z ˆ ˙˙ 2s żzˆ T T D D 2M1 κ´1 α2 2 α α ln ` 2C2 c }upzq ´ u pzq} ď ` }upsq ´ uα psq}2 ds. α ln 2c T ` ˘˘ ´1 p q α α T c 0 α ln D α Applying Gronwall’s inequality to the latter estimate, we get ˆ

α ln

ˆ ˙˙ 2z T T P12 α2 }upzq ´ uα pzq}2 ď ` ` ˘˘ 2c p T ´1q , α T c α ln D α

where P1 “

a

2M1 κ´1 exp p2C2 cq .

(41)

Therefore, we obtain α

}upzq ´ u pzq} ď P1 α

c´z T

9

ˆ

ˆ ln

D α

˙˙´ T `z´c T

.

(42)

By virtue of the well-posedness studied in Theorem 1, we obtain }u

ε,α

˙˙´ z ´ ¯ T D pzq ´ u pzq} ď Q1 α ln }g ε ´ g} ` }hε ´ h} α ˆ ˆ ˙˙´ z T D , ď Q1 ε α ln α ˆ

α

ˆ

(43)

due to (17) and (3). Combining (42) and (43), we obtain }uε,α pzq ´ upzq} ď }uε,α pzq ´ uα pzq} ` }uα pzq ´ upzq} ˆ ˆ ˙˙´ T `z´c ˆ ˆ ˙˙´ z T T c´z D D ď Q1 ε α ln ` P1 α T ln . α α The proof of (1) is completed. Proof of (2). This part is similar to the proof above. To estimate }uα pzq ´ upzq}, we substitute (38) into (37); it yields Hmn “

α e´pγ`c´zqλmn pλmn u ˆmn pzq ` Bz u ˆmn pzqq epγ`c´zqλmn . 2 αλmn ` e´T λmn

In this case, we set up the problem (10)-(12) with T ě c ` γ. Applying Lemma 1 with p :“ γ ` c ´ z and q :“ T gives ¸2 e´pγ`c´zqλmn |λmn u ˆmn pzq ` Bz u ˆmn pzq|2 e2pγ`c´zqλmn αλmn ` e´T λmn ˆ ˆ ˙˙ 2pγ`c´zq ´2 T 2pγ`c´zq D |λmn u ln ď B2α T ˆmn pzq ` Bz u ˆmn pzq|2 e2pγ`c´zqλmn . α

α2 |Hmn | ď 4 2

˜

Hence, we get the following estimate 8 8 ÿ ÿ

2

|Hmn | ď M2 α

2pγ`c´zq T

m“1 n“1

ˆ

ˆ ln

D α

˙˙ 2pγ`c´zq ´2 T

,

ř ř8 2pγ`c´zqλmn |λmn u ˆmn pzq ` Bz u ˆmn pzq|2 ă 8 due to the assumption where M2 “ supzPr0,cs 8 m“1 n“1 e (30). Next, applying Parseval’s identity to (35) and using the estimate (16), we deduce that 8 8 2 ÿ ÿ |Hmn |2 ` 2 }J α pu, zq ´ J α puα , zq}2 κ m“1 n“1 ˆ ˙˙´ 2pz´sq żzˆ T D 2M2 κ´1 α2 }upsq ´ uα psq}2 ds. ď` ` D ˘˘2´ 2pγ`c´zq ` 2C2 c 0 α ln α T α ln α

}upzq ´ uα pzq}2 ď

Thus, we get ˆ

α ln

ˆ

D α

˙˙ 2z T

2pγ`cq

2M2 κ´1 α T }upzq ´ u pzq} ď ` ` ˘˘ 2pγ`cq ` 2C2 c 2´ T ln D α α

2

żzˆ

α ln

0

ˆ

D α

˙˙ 2s T

}upsq ´ uα psq}2 ds.

Applying Gronwall’s inequality to the latter estimate, we obtain ˆ ˆ ˙˙ 2z ˆ ˆ ˙˙ 2pγ`cq ´2 T T T D 2 α 2 2pγ`cq α ln }upzq ´ u pzq} ď P2 α T ln , α α

10

where P2 “

a

2M2 κ´1 exp p2C2 cq .

(44)

Therefore α

}upzq ´ u pzq} ď P2 α

γ`c´z T

ˆ

ˆ ln

D α

˙˙ γ`c´z ´1 T

.

(45)

Moreover, combining (45) with (43), we get }uε,α pzq ´ upzq} ď }uε,α pzq ´ uα pzq} ` }uα pzq ´ upzq} ˆ ˆ ˙˙ γ`c´z ´1 ˆ ˆ ˙˙´ z T T γ`c´z D D T ln ` P2 α . ď Q1 ε α ln α α The proof of (2) is completed. In addition, we also conclude that the exact solution u satisfying the given assumptions is unique due to the completeness of Cpr0, cs; L2 pΩqq. The proof of Theorem 2 is complete.

3

Numerical experiments

Throughout this section, the Cauchy data pg ε , hε q in (3) is defined on a mesh such that their pointwise values are differed randomly with values of the exact data pg, hq in (4). This amounts to modeling discrete data so that the accuracy of numerical calculation can be controllable and guaranteed by a proper mesh resolution. This part is described as follows.

3.1

Modeling of data

Assume that the exact data pg, hq belong to H 2 pΩq ˆ H 2 pΩq. The main purpose of this subsection is to construct the noisy data for pg, hq . Using a uniform rectangular grid with a resolution of M ˆ N in the xy-plane, which is defined by nodal interior points pxi , yj q as xi “ iδx ´ a, yj “ jδy ´ b,

2a , i “ 1, M , M P N, M `1 2b , j “ 1, N , N P N. δy “ N `1 δx “

(46) (47)

we define the data input ε gij “ gij ` ε1,i,j ,

hεij “ hij ` ε2,i,j ,

(48)

where gij “ gpxi , yj q, hij “ hpxi , yj q, ε1,i,j and ε2,i,j play the role of random noise of measurement. ε , hε q, we construct the noisy data pg ε , hε q as follows Now using pgij ij g ε px, yq “

N M ÿ ÿ m“1 n“1

ε

h px, yq “

N M ÿ ÿ m“1 n“1

ε gpmn sin

ˆ

p hεmn sin

ˆ

mπpx ` aq 2a

˙

mπpx ` aq 2a

˙

ˆ sin

ˆ sin

nπpy ` bq 2b

˙

nπpy ` bq 2b

˙

,

(49)

,

(50)

where ε gpmn

p hεmn

˙ ˆ ˙ ˆ M N nπj mπi 2 ÿÿ ε 2 sin , “ g sin M ` 1 N ` 1 i“1 j“1 ij M `1 N `1 ˙ ˆ ˙ ˆ M N nπj mπi 2 ÿÿ ε 2 sin . “ h sin M ` 1 N ` 1 i“1 j“1 ij M `1 N `1

We obtain that pg ε , hε q as in (49), (50) is a pertubed data of pg, hq by the following Theorem.

11

(51)

(52)

Theorem 3. We have the following estimate }g ε ´ g} ` }hε ´ h} ď ε, ? a a where ε “ 4ε0 ab ` 2Cpδx2 ` δy2 q maxt }gxx }2 ` }gyy }2 , }hxx }2 ` }hyy }2 u and C is the positive constant in Lemma 3. Proof. First, according to the results of Hesthaven et al. [14], Chapter 2, we have the following Lemma. Lemma 3. Let v P H 2 pΩq X H01 pΩq. Denote vij “ vpxi , yj q. A function v˜ is called the interpolating function of data sets tvij u, if it satisfies v˜px, yq “

N M ÿ ÿ

ˆ

mπpx ` aq 2a

˙

ˆ

nπpy ` bq 2b

˙

,

(53)

˙ ˆ ˙ ˆ M N nπj mπi 2 ÿÿ 2 sin . “ vij sin M ` 1 N ` 1 i“1 j“1 M `1 N `1

(54)

vˆmn sin

m“1 n“1

sin

where vˆmn Moreover, we have vij “

M ÿ N ÿ

vˆmn sin

ˆ

i“1 j“1

mπi M `1

˙

ˆ sin

nπj N `1

˙

.

(55)

Then the error between v and v˜ is bounded by }˜ v ´ v} ď C

δx2

`

`

δy2

˘b }vxx }2 ` }vyy }2 ,

(56)

where C is a positive constant independently of v, δx and δy . The relationship of vij and vˆmn via Eqs. (54) and (55) can be referred to Press et al. [19], Chapter ˜ be the interpolating functions of data sets tgij u and thij u, respectively. Using discrete 12. Let g˜ and h form of Parseval’s identity, we obtain N M ÿ ÿ m“1 n“1

ε |ˆ gmn ´ gˆmn |2 “

M N ˇ2 2 ÿ ÿ ˇˇ ε 2 gij ´ gij ˇ . M ` 1 N ` 1 i“1 j“1

Thus, combining the latter equation with the definition (53), we obtain }g ε ´ g˜}2 “ ab

N M ÿ ÿ m“1 n“1

ε |ˆ gmn ´ gˆmn |2 “

N M ÿ ÿ 4ab |ε1,i,j |2 pM ` 1qpN ` 1q m“1 n“1

˜ so far the due to Parseval’s theorem and Eq. (48). Thus, applying the same procedure to hε ´ h, latter equation yields ? ? ˜ ď 2ε0 ab, }hε ´ h} }g ε ´ g˜} ď 2ε0 ab, where ε0 “ maxi,j t ε1,i,j , ε2,i,j u. Now using Eq. (56) and the triangle inequality, we get ˜ ` }h ˜ ´ h} g ´ g} ` }hε ´ h} }g ε ´ g} ` }hε ´ h} ď }g ε ´ g˜} ` }˜ b ? ? ` ˘ ď 2ε0 ab ` C δx2 ` δy2 }gxx }2 ` }gyy }2 ` 2ε0 ab ` ˘b ` C δx2 ` δy2 }hxx }2 ` }hyy }2 , and this completes the proof of Theorem 3.

12

3.2

Calculation procedures

This section is devoted to numerical procedures of the regularized problem (10)-(12), which is a nonlinear Volterra integral equation of the second kind. We shall adopt the Picard iteration to solve ε , hε q we need to approximate both the Fourier coefficients and it numerically. For a given datum pgij ij the double summation in Eq. (10). Since we used the trigonometric polynomial approximation to ε , hε q, the Fourier coefficients pp ε ,p gmn hεmn q in Eq. (11) can model the Cauchy data pg ε , hε q from pgij ij be calculated exactly by the formula (54). This fact suggests us to adopt the main idea of Fouriercollocation method [14], in which we equate Eq. (10) at grid points for each iteration step. In other words, we seek for a numerical solution in form of a finite Fourier expansion, and thus, the double summation in Eq. (10) can be replaced by its truncation formed by Eq. (53). Efficiency of the method is relied on the FFT technique. The main duty of this section is splited into two tasks: the first task is to explain the calculation of the fast 2D sine transforms in details, the second task is to approximate the numerical integration (12). In our practice, the calculation of Eqs. (54) and (55) is performed in a natural way. For instance, the sine transform pvij q ÞÑ pˆ vmn q (Eq. (54)) can be computed into two steps: Step 1: Loop i “ 1, M , wni

˙ ˆ N nπj 2 ÿ , @n “ 1, N . :“ vij sin N ` 1 j“1 N `1

(57)

Step 2: Loop n “ 1, N , vˆmn :“

˙ ˆ M mπi 2 ÿ , @m “ 1, M . wni sin M ` 1 i“1 M `1

(58)

Here subroutines sint1f and sint1b of FFTPACK5 [18] can be applied for this calculation. For instance, given a two-dimensional array Vpi, jq size of M ˆ N storing values of tvij u, the following pseudoFORTRAN code is to compute (57) Do i = 1 to M call sint1f ( N, M, V(i,1), M*N, List of auxiliaries variables ) end do and similarly, to compute (58), we invoke Do n = 1 to N call sint1f ( M, 1, V(1,n), M, List of auxiliaries variables ) end do Herein the first i-loop computes a sequence of N elements with M -increment between locations in the array Vp:, :q with size of M ˆ N , in which beginning at Vpi, 1q for each i index, this procedure calculates wmi (Eq. (57)) and returns its values to the V’s storage. Besides, the second n-loop computes a sequence of M elements stored continuously at Vp1, nq for each n index, this procedure computes vˆmn (Eq. (58)) and returns its values to the V’s storage. This performance is based on the array model of FORTRAN programming language, in which an element Vpi, jq of the two-dimensional array Vp:, :q is stored at location index i ` pj ´ 1qM of one-dimensional storage Vp:q with length of M N , this fact holds for any i “ 1, M and j “ 1, N . Thus, the total computational burden in both i- and n-loops (Eqs. (57) and (58)) is appreciated by ¨ OpM log M q „ OpM N logpM N qq, M ¨ OpN log N q ` N looooooooomooooooooon loooooooomoooooooon i-loop

n-loop

which is equal to the number of operations on a vector with M N components. Finally, the calculation of the inverse transform pˆ vmn q ÞÑ pvij q (Eq. (55)) can be performed with subroutine sint1b by the same manner.

13

As shown in the theoretical section, the numerical solution of the problem (10)-(12) can be found by the fixed-point iteration, which is proceeded until convergence met. For instance, to calculate a α pu, zq (Eq. (12)) from the previous u-profile at any u-profile, we need to compute the integral Jˆmn z P r0, T s. Therefore, the computation should be performed on a fixed mesh in z-direction. We define a uniform mesh as follows T zk “ pk ´ 1qδz , δz “ , k “ 1, K, K P N. (59) K ´1 α pu, z q can be simplified as the following For each pm, nq P N2 , the integral Jˆmn k ż zk α pu, zk q “ wpsq ds, Jˆmn

(60)

z1

where wplq is defined by 1 wl :“ 2λmn

˜

e´pT ´zk `zl qλmn ´ epzl ´zk qλmn αλmn ` e´T λmn

¸ fˆmn pupzl q, zl q, @l “ 1, k.

(61)

Now we need a procedure to estimate the integral (60) from twl u above. According to traditional α pu, z q by the methods for integral equations (e.g. see [19]), for k “ 2 and k “ 3, we estimate Jˆmn k trapezoidal and Simpson rules, respectively. For k ě 4, first we fit the integrand w from its discrete values pwl q by using Hermite interpolation, i.e. wpsq “ h1 prqwj ` h2 prqδwj ` h3 prqwj`1 ` h4 prqδwj`1 , s P rzj , zj`1 s, for j “ 1, k ´ 1 and r “

s´zj zj`1 ´zj ,

where $ if l “ 1, & w 2 ´ w1 , 1 pw ´ w q, if 2 ď l ď k ´ 1, δwl “ l´1 % 2 l`1 if l “ k. wk ´ wk´1

The Hermite basis functions h1 , h2 , h3 , h4 are given by h1 prq “ 2r3 ´ 3r2 ` 1,

h2 prq “ r3 ´ 2r2 ` r,

h3 prq “ ´2r3 ` 3r2 ,

h4 prq “ r3 ´ r2 .

It should be mentioned that we have wpzl q “ wl and w1 pzl q “ δwl {δz for l “ 1, k. Second, substituting α pu, z q can be approximated wpsq into Eq. (60) and calculating polynomial integrations manually, Jˆmn k by $ 11w1 ` 14w2 ´ w3 , if j “ 2, k ż zj k ÿ ÿ δz & ´wj´2 ` 13wj´1 ` 13wj ´ wj`1 , if 3 ď j ď k ´ 1, ¨ wpsq ds “ 24 % j“2 zj´1 j“2 if j “ k. ´wk´2 ` 14wk´1 ` 11wk , To close this section, Eq. (10) is identified on the collocation grid pxi , yj , zk q by the following uij pzk q “ Gαij pg, h, zk q ` Jijα pu, zk q,

(62)

for every i “ 1, M , j “ 1, N and k “ 1, K. Here again the subscript ”ij” denotes the values of u, Gα , J α at the grid point pxi , yj q. Alternatively, we can replace the latter equation by its sine transformation to frequency domain, i.e. α ˆ αmn pg, h, zk q ` Jˆmn u ˆmn pzk q “ G pu, zk q,

(63)

for every m “ 1, M , n “ 1, N and k “ 1, K. In that way, Eqs. (54) and (55) were adopted to transfer data between the collocation grid and the frequency domain. Finally, to calculate (61), we used Eq. (54) as follows ˙ ˆ ˙ ˆ M ÿ N ÿ nπj mπi 4 ˆ sin . f pupxi , yj , zl q, xi , yj , zl q sin fmn pupzl q, zl q “ pM ` 1qpN ` 1q i“1 j“1 M `1 N `1

14

3.3

Test cases

In this section, we propose two test cases in differed scenarii to verify the stability and convergence of the regularized solution corresponding to both of exact and disturbed data. Firstly, to set up a test case, we choose a test function U that plays the role of the exact solution, i.e. gpx, yq “ U px, y, 0q,

(64)

hpx, yq “ Bz U px, y, 0q.

(65)

The source function f is determined such that U satisfies Eq. (1). Secondly, to illustrate the sensitivity of the computational accuracy to noise of the data, we repeated calculations with a variety of perturbed data. The perturbation was defined as ε randp¨, ¨q where each random term randp¨, ¨q was uniformly determined on r´1, 1s, and ε plays the role of an amplitude, i.e. ε gij “ gpxi , yj q ` ε randpxi , yj q,

(66)

hεij “ hpxi , yj q ` ε randpxi , yj q,

(67)

for any i “ 1, M , j “ 1, N , where xi , yj and M, N were described by (46) and (47). Let zk be defined by (59). The aim of the numerical experiments is to observe a relative errors δ ε,α , i.e. g f řN řM f i“1 j“1 |uα pxi , yj , cq ´ U pxi , yj , cq|2 ε,α δ pcq :“ e , (68) řN řM 2 i“1 j“1 |U pxi , yj , cq| as α tends to zero in the following cases of ε: ε , hε q is the exact data. 1. ε “ 0. The couple of pgij ij ε , hε q plays the role of the measured data with a random noise. 2. ε ą 0. The couple of pgij ij

Next we introduce two examples. Example 1:

Choose ˙ ˆ ˙„ ˆ ˙ ˆ ˙j ˆ 2πpy ` bq 5πpx ` aq 6πpy ` bq πpx ` aq sin 1 ` z cos cos . U px, y, zq “ sin 2a 2b 2a 2b

(69)

The source function f is taken in the sine-Gordon form as follows f pu, x, y, zq “ sinpuq ` Rpx, y, zq, where R “ ΔU ´ sinpU q. The domain of problem Ω ˆ r0, cs is defined by a “ b “ 5 and c “ 1.5. Mesh resolution is M ˆ N ˆ K “ 511 ˆ 511 ˆ 511. Here we can see that the test function U p¨, ¨, zq shown in Eq. (69) satisfies not only Eqs. (1) and (2) but also the assumption (28), since it can be presented by a finite expansion of the sine functions. Setting up the problem (10)-(12) with T “ c, we expect ε , hε q can converge stably to the exact that the regularized solution uα obtained from given datum pgij ij solution U as ε tends to zero. Example 2:

We setup a test case using the function U px, y, zq “

4 ÿ

e´ri ppx´xi zq

2 `py´y

i zq

2

q cospx2 ` y 2 q,

(70)

i“1

where r1 “ 1, ri`1 “ ri ` 0.1 and x1 “ x4 “ y1 “ y2 “ ´1.1, x2 “ x3 “ y3 “ y4 “ 1.1. We define the Cauchy data by Eqs. (66) and (67), and the source function f pu, x, y, zq “

u ` Rpx, y, zq, u2 ` 1

15

where R “ ΔU ´ U {pU 2 ` 1q. In this example, we consider the problem (10)-(12) in domain the Ω ˆ r0, cs (i.e. T “ c), where a “ b “ 5 and c “ 1; the mesh sizes are given by M “ N “ K “ 511. It is easy to see that the test function U shown in Eq. (70) does not satisfy the homogeneous Dirichlet boundary condition (2) for every z P r0, cs. Moreover, U is not the analytical solution of problem ¯ ˆ r0, cs. Therefore, assumption (28) is not true in this case. We have to find the solution (1)-(4) in Ω ε , hε q, we U of well-posed problem (10)-(12) by an approximating method. Using given data set pgij ij construct approximate solution uε,α for approximation of U . Figure 1 shows graphs of test functions, i.e. Figs. 1(a) and 1(b) show U p¨, ¨, cq in Eqs. (69) and (70), respectively. These figures will be used to compare with numerical solutions. Figure 2 shows log10 scaled relative error estimate (log10 pδuε,α q) of two contiguous solutions produced by the fixed-point iteration for solving the problem (10)-(12), wherein g ˇ f řN řM ˇˇ α f i“1 j“1 uold pxi , yj , cq ´ uαnew pxi , yj , cqˇ2 ε,α e , (71) δu “ ˇ2 řN řM ˇˇ α ˇ i“1 j“1 uold pxi , yj , cq for both ε “ 0 and ε ą 0 (ε “ 5 ¨ 10´2 , 1 ¨ 10´2 , 5 ¨ 10´3 , 1 ¨ 10´3 , 5 ¨ 10´4 , 1 ¨ 10´4 ), α “ 10´s for s “ 1, 7; here we initialized uα “ 0 for any cases. In the figure, the first and second columns show graphs of δuε,α from Examples 1 and 2. For each sub-figure corresponding to a value of ε, the solid lines show graphs of δuε,α corresponding to α “ 10´s , for all s “ 1, 7; here the calculating processes were terminated as δuε,α ă 10´9 . As shown in the figure, the reduction rate of δuε,α does not depend so much on either the quality of data (ε) or the magnitude of regularization parameter (α), since the curves are shaped almost linearly and parallel with each other among sub-figures. Quantitatively, the value of δuε,α is upper estimated by β ¨ 10´κl , where β is a constant, κ « 9{8, and l is number of iteration. In these examples, we can say that contraction rate of the operator F in Eq. (18) was proportional to 10´9{8 . Tables 1 and 2 show the relative error estimate δ ε,α of Example 1 and 2, respectively. Here δ ε,α was defined by Eq. (68), the 3D mesh resolution was determined by M “ N “ K “ 511. The error estimate in these tables were calculated with both cases of Cauchy data: ε “ 0 and ε ą 0 (ε “ 5 ¨ 10´2 , 1¨10´2 , 5¨10´3 , 1¨10´3 , 5¨10´4 , 1¨10´4 ); which are arranged from the left to right of Tables. For each such a noise amplitude ε, from the top to bottom of Tables, the regularization parameter α decreased from 10´1 to 10´7 . This arrangement is to show the fact that was predicted in the theoretical section: the error estimate δ ε,α converges to zero while α decreased until α “ ε; for smaller values of α, the convergence of regularized solution might not be guaranteed, or even, the divergence is attained. This phenomena is also simulated by Figures 3 and 4 as the following. Figure 3 shows the numerical solution of Example 1. Here the first and second columns show the graph and its contours of the regularized solution corresponding to the exact and disturbed data, i.e. ε “ 0 and ε ą 0. As shown in the first column (ε “ 0), the regularized solution obtained from the exact data converges to the exact solution stably as α tends to zero. However, in the second column (ε ą 0), the regularized solution calculated from the disturbed data converges to the exact solution only when α ě ε. For α ă ε, it start to diverge due to the collapse of regularization strategy. Figure 4 shows the numerical solution of Example 2. Similarly to Example 1, the first and second columns show the graph and its contours of the regularized solution corresponding to the exact and disturbed data. In the first column (ε “ 0), the regularized solution calculated from the exact data is getting closer to the test function U quite well as α decreases to zero. Meanwhile, not similarly to Example 1, the contours within sub-figures in the second column (ε ą 0) show a sensitive dependence of data disturbance of the regularized solution, e.g. their contours were strongly fluctuated near the boundaries. Although this phenomena was not seen in Fig. (3) of Example 1, the convergence of uα in this example shares the same process with the previous example: the regularized solution uα converges stably to U as α tends to zero and α ě ε, it will diverge once α ă ε.

16

4

Conclusion

A regularization method for a quasi-linear elliptic Cauchy problem has been successfully proposed in this paper. We stress that the regularized problem is well-posed, and its solution converges to the exact solution strongly in L2 where some a priori assumptions are considered. The remarkable thing is that the sets of functions in L2 that satisfy such conditions have a direct relation to abstract Gevrey spaces. Thus, our analysis is more concrete. In order to see how our method works in numerical procedures, two examples are provided under many powerful tools and the implementation is programmed by FORTRAN.

(a) Example 1, Eq. (69)

(b) Example 2, Eq. (70)

Figure 1: Graphs of the test function U px, y, cq in two examples shown with mesh size of M ˆ N “ 511 ˆ 511.

aa

ε α aa a 1.0E-01 1.0E-02 1.0E-03 1.0E-04 1.0E-05 1.0E-06 1.0E-07

0.0E+00

5.0E-02

1.0E-02

5.0E-03

1.0E-03

5.0E-04

1.0E-04

3.68E-01 2.79E-01 1.21E-01 1.83E-02 1.97E-03 2.02E-04 2.03E-05

3.68E-01 2.80E-01 3.65E-01 3.32E+00 3.22E+01 3.13E+02 3.05E+03

3.68E-01 2.79E-01 1.42E-01 6.75E-01 6.51E+00 6.32E+01 6.14E+02

3.68E-01 2.79E-01 1.27E-01 3.36E-01 3.25E+00 3.15E+01 3.06E+02

3.68E-01 2.79E-01 1.22E-01 6.97E-02 6.52E-01 6.33E+00 6.16E+01

3.68E-01 2.79E-01 1.21E-01 3.82E-02 3.23E-01 3.13E+00 3.05E+01

3.68E-01 2.79E-01 1.21E-01 1.95E-02 6.49E-02 6.30E-01 6.14E+00

Table 1: Example 1, relative error δ ε,α pcq defined by Eq. (68). Mesh resolution is M ˆ N ˆ K “ 511 ˆ 511 ˆ 511

aa

ε α aa a 1.0E-01 1.0E-02 1.0E-03 1.0E-04 1.0E-05 1.0E-06 1.0E-07

0.0E+00

5.0E-02

1.0E-02

5.0E-03

1.0E-03

5.0E-04

1.0E-04

4.85E-01 3.62E-01 1.69E-01 4.37E-02 8.96E-03 1.90E-03 7.14E-04

4.85E-01 3.82E-01 1.23E+00 1.16E+01 1.11E+02 1.07E+03 1.04E+04

4.86E-01 3.63E-01 2.99E-01 2.32E+00 2.23E+01 2.16E+02 2.09E+03

4.86E-01 3.62E-01 2.07E-01 1.16E+00 1.11E+01 1.08E+02 1.04E+03

4.86E-01 3.62E-01 1.71E-01 2.37E-01 2.24E+00 2.16E+01 2.09E+02

4.86E-01 3.62E-01 1.70E-01 1.24E-01 1.11E+00 1.07E+01 1.04E+02

4.85E-01 3.62E-01 1.69E-01 4.97E-02 2.22E-01 2.15E+00 2.09E+01

Table 2: Example 2, relative error δ ε,α pcq defined by Eq. (68). Mesh resolution is M ˆ N ˆ K “ 511 ˆ 511 ˆ 511.

17

1.0E+04

Log base 10 scale the relative error between new and old solution

Log base 10 scale the relative error between new and old solution

1.0E+02

1.0E+00

1.0E-02

1.0E-04

1.0E-06

1.0E-08

1.0E-10

1.0E-12

1.0E+02

1.0E+00

1.0E-02

1.0E-04

1.0E-06

1.0E-08

1.0E-10

1.0E-12 0

2

4

6

8

10

0

2

Numeber of iteration

(a) Example 1, ε “ 0.

8

10

8

10

8

10

1.0E+02

Log base 10 scale the relative error between new and old solution

Log base 10 scale the relative error between new and old solution

6

(b) Example 2, ε “ 0.

1.0E+04

1.0E+02

1.0E+00

1.0E-02

1.0E-04

1.0E-06

1.0E-08

1.0E-10

1.0E-12

1.0E+00

1.0E-02

1.0E-04

1.0E-06

1.0E-08

1.0E-10

1.0E-12 0

2

4

6

8

10

0

2

Numeber of iteration

4

6

Numeber of iteration

(c) Example 1, ε “ 10´2 .

(d) Example 2, ε “ 10´2 . 1.0E+02

Log base 10 scale the relative error between new and old solution

1.0E+04

Log base 10 scale the relative error between new and old solution

4

Numeber of iteration

1.0E+02

1.0E+00

1.0E-02

1.0E-04

1.0E-06

1.0E-08

1.0E-10

1.0E-12

1.0E+00

1.0E-02

1.0E-04

1.0E-06

1.0E-08

1.0E-10

1.0E-12 0

2

4

6

8

10

0

Numeber of iteration

2

4

6

Numeber of iteration

(e) Example 1, ε “ 10´3 .

(f) Example 2, ε “ 10´3 .

Figure 2: Log10 scale the relative error estimate δuε,α defined in Eq. (71). The calculation process is terminated as δuε,α ă 10´9 .

18

(a) ε “ 0, α “ 10´1

(b) ε “ 10´2 , α “ 10´1

(c) ε “ 0, α “ 10´2

(d) ε “ 10´2 , α “ 10´2

(e) ε “ 0, α “ 10´3

(f) ε “ 10´2 , α “ 10´3

Figure 3: Example 1. Graphs of uε,α px, y, cq with mesh size of M ˆ N “ 511 ˆ 511.

19

(a) ε “ 0, α “ 10´1

(b) ε “ 10´2 , α “ 10´1

(c) ε “ 0, α “ 10´2

(d) ε “ 10´2 , α “ 10´2

(e) ε “ 0, α “ 10´3

(f) ε “ 10´2 , α “ 10´3

Figure 4: Example 2. Graphs of uε,α px, y, cq with mesh size of M ˆ N “ 511 ˆ 511.

20

Acknowledgment This research is funded by Foundation for Science and Technology Development of Ton Duc Thang University (FOSTECT), website: ”http://fostect.tdt.edu.vn”, under Grant FOSTECT.2016.BR.20. The authors also desire to thank the handling editor and anonymous referees for their helpful comments on this paper. The corresponding author gratefully acknowledge stimulating discussions with Vo Anh Khoa in Gran Sasso Science Institute, Italy.

References [1] Faker Ben Belgacem, Why is the Cauchy problem severely ill-posed?, Inverse Problems, volume 23, number 2, p. 823, 2007, IOP Publishing. [2] Bin-Mohsin, B.; Lesnic, D, Numerical reconstruction of an inhomogeneity in an elliptic equation Inverse Probl. Sci. Eng. 22 (2014), no. 1, 184–198. [3] Bin-Mohsin, B.; Lesnic, D, The method of fundamental solutions for Helmholtz-type equations in composite materials. Comput. Math. Appl. 62 (2011), no. 12, 4377–4390. [4] J. Hadamard, Lectures on the Cauchy Problem in Linear Differential Equations, Yale University Press, New Haven, CT, 1923. [5] H. Jol and D. Smith, Ground Penetrating Radar of Northern Lacustrine Delta: Canadian Journal of Earth Sciences, 28, 1939-1947, 1991. [6] G. Hogan, Migration of Ground Penetrating Radar Data: a Technique for Locating Subsurface Targets. Proceedings of the Symposium on the Application of Geophysics to Engineering and Environmental Problems, U.S. Goeligical Survey, Golden, CO., U.S.A., March 28-31, Society of Engineering and Mineral Geophysics, pp. 809-822, 1988. [7] T. Reginska, A. Wakulicz, Wavelet moment method for the Cauchy problem for the Helmholtz equation, Journal of Computational and Applied Mathematics, Volume 223, 2009, Pages 218–229. [8] T. Reginska, T. Reginski, Approximate solution of a Cauchy problem for the Helmholtz equation, Inverse Problems 22 (2006), 975–989. [9] L. Elden, V. Simoncini, A numerical solution of a Cauchy problem for an elliptic equation by Krylov subspaces, Inverse Problems, volume 25, number 6, pages 065002, 2009, IOP Publishing. [10] Nguyen Huy Tuan, Le Duc Thang, Dang Duc Trong, and Vo Anh Khoa, Approximation of mild solutions of the linear and nonlinear elliptic equations, Inverse Probl. Sci. Eng. 23 (2015), no. 7, 1237–1266. [11] Nguyen Huy Tuan, Le Duc Thang, and Vo Anh Khoa, A modified integral equation method of the nonlinear elliptic equation with globally and locally Lipschitz source, Applied Mathematics and Computation, volume 265, pages 245–265, 2015. [12] Nguyen Huy Tuan , Quoc Viet Tran, Van Thinh Nguyen, Some remarks on a modified Helmholtz equation with inhomogeneous source, Appl. Math. Model. 37 (2013), no. 3, 793–814. [13] Quoc Viet Tran, Nguyen Huy Tuan , Van Thinh Nguyen, and Duc Trong Dang, A general filter regularization method to solve the three dimensional Cauchy problem for inhomogeneous Helmholtz-type equations: theory and numerical simulation, Appl. Math. Model. 38 (2014), no. 17-18, 4460–4479. [14] J. Hesthaven, S. Gottlieb and D. Gottlieb, Spectral methods for time-dependent problems, Cambridge UP, Cambridge, UK, 2007.

21

[15] D.R. Kincaid and E.W. Cheney, Numerical analysis: mathematics of scientific computing, volume 2, 2002, American Mathematical Society. [16] D.D. Trong, P.H. Quan, and N.H. Tuan, A quasi-boundary value method for regularizing nonlinear ill-posed problems, Electronic J. Differential Equations 2009(109) (2009), pp. 1-16. [17] N.H. Tuan and D.D. Trong, A nonlinear parabolic equation backward in time: Regularization with new error estimates, Nonlinear Anal. Theor. Meth. Appl. 73(6) (2010), pp. 1842-1852. [18] P. N. Swarztrauber, FFTPACK5, Computational Information Systems Laboratory, University Corporation for Atmospheric Research, 2011. [19] W.H. Press et al., Numerical Recipes in Fortran 90, 2nd ed., Cambridge University Press, New York, 1996.

22