- Email: [email protected]
A numerical method for solving a Goursat–Cauchy boundary value problem
Applied Mathematics and Computation 220 (2013) 123–141
Contents lists available at SciVerse ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
A numerical method for solving a Goursat–Cauchy boundary value problem Mihaela-Cristina Drignei Division of Physical and Computational Sciences, University of Pittsburgh at Bradford, Bradford, PA 16701, USA
a r t i c l e
i n f o
a b s t r a c t
Keywords: Second order hyperbolic partial differential equation Inverse Sturm–Liouville problem Goursat boundary value problem Cauchy boundary value problem
In this paper we propose a method to construct numerically the solution-pair of a second order hyperbolic partial differential equation when two types of boundary conditions on a triangular domain are imposed. Specifically, the PDE is W xx ðx; tÞ W tt ðx; tÞ qðxÞWðx; tÞ ¼ dqðxÞLðx; tÞ, where the unknown functions are Wðx; tÞ and dqðxÞ. One set of boundary conditions refer to one set of the triangle’s boundaries (we call them Goursat type), and the other set of boundary conditions both refer to the third boundary of the triangle (we call them Cauchy type). To solve numerically this boundary value problem we make a change of dependent variables and place a mesh on the triangular domain with lines that are the characteristic curves of the hyperbolic PDE. Then we shall integrate along the characteristic curves of the hyperbolic PDE in such a way that numerical values of Wðx; tÞ and dqðxÞ will be produced at mesh nodes located on vertical lines. The calculations will proceed from the vertical boundary of the domain going backwards to the origin O of the Oxt coordinate system, from one vertical line to another. The trapezoidal formula will be used for numerical integration. Numerical examples will be provided. Ó 2013 Elsevier Inc. All rights reserved.
1. Introduction Our aim in this paper is to construct numerically the solution-pair fdqðxÞ; Wðx; tÞg of the boundary value problem in the triangle D0 :¼ fðx; tÞ 2 R2 : 0 6 t 6 x 6 ag (see Fig. 1):
W xx ðx; tÞ W tt ðx; tÞ qðxÞWðx; tÞ ¼ dqðxÞLðx; tÞ; ðx; tÞ 2 D0 ; Wðx; xÞ ¼
1 2
Z
ð1:1Þ
x
dqðsÞ ds;
0 6 x 6 a;
ð1:2Þ
0
Wðx; 0Þ ¼ 0;
0 6 x 6 a;
Wða; tÞ ¼ f ðtÞ; W x ða; tÞ ¼ gðtÞ;
ð1:3Þ
0 6 t 6 a;
ð1:4Þ
0 6 t 6 a;
ð1:5Þ 2
1
2
where qðxÞ; Lðx; tÞ; f ðtÞ; gðtÞ are known functions in the spaces L ð0; aÞ; CðD0 Þ; H ð0; aÞ; L ð0; aÞ, respectively, and f ð0Þ ¼ 0 ¼ gð0Þ. They are assumed real-valued. Here, E-mail address: [email protected] 0096-3003/$ - see front matter Ó 2013 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.amc.2013.05.041
124
M.-C. Drignei / Applied Mathematics and Computation 220 (2013) 123–141
Fig. 1. The domain of the problem (1.1)–(1.5).
CðD0 Þ is the space of continuous functions on D0 ; L2 ð0; aÞ :¼
u:
Z
a
u2 ðsÞ ds < 1 ;
0
Z H1 ð0; aÞ :¼ u 2 C½0; a : uðxÞ ¼ a þ
x
v ðsÞ ds;
for some a 2 R and some
v 2 L2 ð0; aÞ real-valued
:
0
Note that the requirement f ð0Þ ¼ 0 ¼ gð0Þ is needed in order to make (1.4) and (1.5) compatible with (1.3). The existence and uniqueness of the solution-pair fdqðxÞ; Wðx; tÞg to the problem (1.1)–(1.5) were proved in [[20] Theorem 4.17, pp. 152–154]. Our search for the numerical solution to the above boundary value problem was motivated by the fact that such a problem is to be solved at each stage of the Newton iteration process in order to obtain the numerical solution qðxÞ to the nonlinear equation:
FðqÞ ¼ f~f ðtÞ; g~ðtÞg;
ð1:6Þ
for a given pair of functions f~f ðtÞ; g~ðtÞg in H1 ð0; aÞ L2 ð0; aÞ. Here the map F is defined by
FðqÞ :¼ fKða; t; qÞ; Kx ða; t; qÞg;
ð1:7Þ
where Kðx; t; qÞ is the week solution to the Goursat problem in the triangle D0 :
8 ðx; tÞ 2 D0 ; > tt ðx; tÞ qðxÞKðx; tÞ ¼ 0; < Kxx ðx; tÞ K Rx Kðx; xÞ ¼ 12 0 qðsÞ ds; 0 6 x 6 a; > : Kðx; 0Þ ¼ 0; 0 6 x 6 a:
ð1:8Þ
In the theory of inverse spectral problems, solving (1.6) is equivalent to finding the potential function qðxÞ of the canonical ~ðtÞg obtained from two known spectral data of this operSturm–Liouville operator from knowledge of the Cauchy data f~f ðtÞ; g ator (see Section Appendix A). Then the Newton iterates are:
q0 ¼ initial guess; qnþ1 ¼ qn þ dqn ;
for n ¼ 0; 1; 2; . . . ;
ð1:9Þ
where fdqn ðxÞ; W n ðx; tÞg is the solution-pair to (1.1)–(1.5) with
8 qðxÞ :¼ qn ðxÞ; > > > < Lðx; tÞ :¼ Kðx; t; q Þ; n ~f ðtÞ Kða; t; q Þ; > f ðtÞ :¼ > n > : gðtÞ :¼ g~ðtÞ Kx ða; t; qn Þ:
ð1:10Þ
Complete details about the application of the Newton method to (1.6) will be provided in a different paper. We associated the name Goursat–Cauchy to the boundary value problem (1.1)–(1.5) because of its resemblance to the problem
8 ðx; tÞ 2 D0 ; > tt ðx; tÞ qðxÞWðx; tÞ ¼ 0; < W xx ðx; tÞ W R 1 x Wðx; xÞ ¼ 2 0 qðsÞ ds; 0 6 x 6 a; > : Wðx; 0Þ ¼ 0; 0 6 x 6 a; known in literature as the Goursat problem in the triangle D0 (see the statement right before Theorem 4.15 in [20]), and to the problem
M.-C. Drignei / Applied Mathematics and Computation 220 (2013) 123–141
125
8 > < W xx ðx; tÞ W tt ðx; tÞ qðxÞWðx; tÞ ¼ 0; ðx; tÞ 2 D0 ; Wða; tÞ ¼ f ðtÞ; 0 6 t 6 a; > : W x ða; tÞ ¼ gðtÞ; 0 6 t 6 a; known in literature as the Cauchy problem in the triangle D0 (see the statement right before Theorem 4.16 in [20]). For the past several decades the following boundary value problem, known as the Goursat problem on a rectangular domain has been abundantly studied (e.g. [1–4,9,12,16,15,17,22,23,25,26,28–31]):
8 > < uxy ¼ f ðx; y; u; ux ; uy Þ; 0 6 x 6 a and 0 6 y 6 b; uðx; 0Þ ¼ /ðxÞ; 0 6 x 6 a; > : uð0; yÞ ¼ wðyÞ; 0 6 y 6 b;
ð1:11Þ
where /ð0Þ ¼ wð0Þ, and various kinds of regularity on the known functions f ; /; w were assumed. In [14], the authors study a Goursat problem on a triangular domain. Specifically, they investigate the numerical resolvability of the problem
8 0 6 x 6 1 and x 6 y 6 1; > < uxy ¼ f ðx; y; u; ux ; uy Þ; uð0; yÞ ¼ sðyÞ; 0 6 y 6 1; > : H1 ux ðx; xÞ þ H2 uy ðx; xÞ þ uðx; xÞ ¼ rðxÞ; 0 6 x 6 1;
ð1:12Þ
where the functions f ; s; r are known and the real constants H1 ; H2 are also known, and point out that their method ’’can be used to approximately solve some inverse problems, which are further applied in areas such as geoscience and biomedical engineering’’. In the appendix of [14] the authors describe the connection between their problem (1.12) and an inverse spectral problem: recovery of the potential function of the canonical Sturm–Liouville operator from knowledge of the potential function only on half of its interval (i.e. ’partially known potential’). One difference between (1.12) and our current problem consists in the fact that they are applicable to different inverse spectral problems: ours is used for the recovery of the potential function from two spectral data via the Cauchy data, and (1.12) is applicable to the ’partially known potential’ case. Another difference is that, even with a change of variables, our problem (1.1)–(1.5) cannot be put into the form of (1.12), nor of (1.11). Indeed, by the change of the independent variables
ðx; tÞ $ ðn; gÞ; given by
(
n ¼ xþt ; 2 g ¼ xt ; 2
our domain D0 is transformed into a new domain D0 :¼ f0 6 g 6 n 6 a; and n þ g 6 ag, which is also a triangle (see Fig. 2), the PDE we study, namely (1.1) is transformed into a new PDE:
wng ðn; gÞ qðn þ gÞwðn; gÞ ¼ dqðn þ gÞlðn; gÞ; ðn; gÞ 2 D0 ;
ð1:13Þ
where wðn; gÞ :¼ Wðx; tÞ and lðn; gÞ :¼ Lðx; tÞ, and the boundary conditions (1.2)–(1.5) are transformed respectively into the following boundary conditions:
wðn; 0Þ ¼
1 2
Z
n
dqðsÞ ds;
0 6 n 6 a;
0
Fig. 2. The domain of the problem (1.13)–(1.17).
ð1:14Þ
126
M.-C. Drignei / Applied Mathematics and Computation 220 (2013) 123–141
wðn; nÞ ¼ 0;
06n6
a ; 2
wðn; a nÞ ¼ f ð2n aÞ;
ð1:15Þ a 6 n 6 a; 2
wn ðn; a nÞ þ wg ðn; a nÞ ¼ 2gð2n aÞ;
ð1:16Þ a 6 n 6 a; 2
ð1:17Þ
Thus, our original problem (1.1)–(1.5) is equivalent to the problem (1.13)–(1.17) for the unknown pair of functions fdqðn þ gÞ; wðn; gÞg. Note that (1.13) is not a particular case of the PDE in (1.11), nor in (1.12). Indeed, the PDE (1.13) is to be solved for two unknowns wðn; gÞ and dqðn þ gÞ, while those in (1.11) and in (1.12) are to be solved for only one, namely uðx; yÞ. The boundary conditions (1.14)–(1.17) refer to all boundaries of the domain in our case, whereas in all of the above mentioned papers they refer only to two boundaries of the rectangular or triangular domains. Also, the methods of numerical resolvability of the proposed problems are different in our paper and in Ref. [14]. In the present paper we make a change of dependent variables
Wðx; tÞ $ fZ 1 ðx; tÞ; Z 2 ðx; tÞg and place a mesh on the triangular domain D0 that consists of lines that are the characteristic curves of the hyperbolic PDE (1.1). We integrate along these characteristic curves to produce numerical values of Z 1 ðx; tÞ; Z 2 ðx; tÞ (from which Wðx; tÞ will Rx be obtained), dqðxÞ and rðxÞ :¼ 0 dqðsÞ ds at the grid points. The integration is performed in such a way that we progress from one vertical line to another, starting from the vertical boundary of the domain D0 and going backwards to the origin O of the Oxt coordinate system. Specifically, our each working cell is a triangle with two vertexes on one vertical line, and the third vertex on the neighboring vertical to the left (details and graphical illustrations are provided in Section 2). We utilize the previously calculated values of Z 1 ; Z 2 ; dq; r at the two vertical nodes of the triangular cell to produce their values at the third vertex of the cell. In Ref. [14], the mesh consists of lines parallel to the coordinate axes (x-axis and y-axis), the working cell is a square, the computations proceed from one horizontal line to another, starting from the origin O and going upwards, and from left to right on each horizontal line. Gou and Sun in [14] first find approximate values of u; ux ; uy at all nodes on the vertical boundary x ¼ 0, then at the mesh node ðx1 ; y1 Þ on the slanted boundary y ¼ x. Then for each j P 2 the computations move from the horizontal line y ¼ yj1 to the horizontal line y ¼ yj : for each square cell spanning between the lines y ¼ yj1 and y ¼ yj , approximate values of u; ux ; uy at the upper-right corner are produced from the previously calculated values at the other three corners. The values of u; ux ; uy at the grid point ðxj ; yj Þ on the slanted boundary are obtained either from the values at the grid points ðxi ; yj Þ, with i ¼ 0; 1; . . . ; j 1, or else from the values at the grid points ðxj1 ; yj Þ and ðxj1 ; yj1 Þ, depending on how the impedance boundary parameters H1 ; H2 are chosen. Theoretical approaches to the existence of global solutions to some non-linear Cauchy–Goursat problems for the wave equation have been developed for example in [18,19]. The equations considered are respectively utt uxx juja ¼ f ðx; tÞ, and utt uxx juja1 u ¼ f ðx; tÞ, in fðx; tÞ : 0 < x < t; 0 < t < Tg, with the conditions ux ð0; tÞ ¼ 0 and uðt; tÞ ¼ 0 for 0 6 t 6 T. The function f ðx; tÞ and the parameter a > 0 are known. 2. The method First, we note that the characteristic curves of the hyperbolic PDE (1.1) are the curves
x þ t ¼ constant; x t ¼ constant:
They are obtained by making a change of independent variables
ðx; tÞ $ ðn; gÞ that puts (1.1) into an equivalent, but easier to integrate form - at least for theoretical purposes. Setting wðn; gÞ :¼ Wðx; tÞ, straightforward calculations based on the change of variables n ¼ nðx; tÞ; g ¼ gðx; tÞ lead to the equivalent PDE:
~ðn; gÞwðn; gÞ ¼ dq ~ðn; gÞlðn; gÞ; ðn2x n2t Þwnn þ 2ðnx gx nt gt Þwng þ ðg2x g2t Þwgg þ ðnxx ntt Þwn þ ðgxx gtt Þwg q
ð2:1Þ
~ðn; gÞ :¼ qðxÞ; dq ~ðn; gÞ :¼ dqðxÞ; lðn; gÞ :¼ Lðx; tÞ. Choosing where q
(
n ¼ xþt ; 2
g ¼ xt ; 2 Eq. (2.1) takes the simpler form (1.13). Hence, we have justified the equations of the characteristic curves. However, we shall not pursue Eq. (1.13), but rather (1.1). We make the change of dependent variables
Wðx; tÞ $ fZ 1 ðx; tÞ; Z 2 ðx; tÞg given by:
ð2:2Þ
M.-C. Drignei / Applied Mathematics and Computation 220 (2013) 123–141
Z 1 ðx; tÞ :¼ W x ðx; tÞ þ aðxÞWðx; tÞ þ W t ðx; tÞ;
127
ð2:3Þ
Z 2 ðx; tÞ :¼ W x ðx; tÞ þ aðxÞWðx; tÞ W t ðx; tÞ; where aðxÞ is to be chosen conveniently later, and set
rðxÞ :¼
Z
x
dqðsÞ ds:
ð2:4Þ
0
From (2.3) we can write based on (1.1) that:
(
ðZ 1;x Z 1;t Þðx; tÞ ¼ ða0 ðxÞ a2 ðxÞ þ qðxÞÞWðx; tÞ þ aðxÞZ 2 ðx; tÞ þ dqðxÞLðx; tÞ; ðZ 2;x þ Z 2;t Þðx; tÞ ¼ ða0 ðxÞ a2 ðxÞ þ qðxÞÞWðx; tÞ þ aðxÞZ 1 ðx; tÞ þ dqðxÞLðx; tÞ:
ð2:5Þ
We shall choose aðxÞ such that the Riccati equation is satisfied:
a0 ðxÞ a2 ðxÞ þ qðxÞ ¼ 0; for 0 6 x 6 a:
ð2:6Þ
Inserting (2.6) into (2.5) and taking t ¼ x þ c in the first equation of (2.5), and t ¼ x c in the second equation of (2.5), we obtain:
d Z 1 ðx; x þ cÞ ¼ aðxÞZ 2 ðx; x þ cÞ þ dqðxÞLðx; x þ cÞ; dx
for 0 6 x 6 a;
ð2:7Þ
and respectively,
d Z 2 ðx; x cÞ ¼ aðxÞZ 1 ðx; x cÞ þ dqðxÞLðx; x cÞ; dx
for 0 6 x 6 a:
ð2:8Þ
Utilizing (1.4) and (1.5) in (2.3) we obtain two new boundary conditions for Z 1 and Z 2 on the vertical boundary of the triangular domain D0 :
Z 1 ða; tÞ ¼ gðtÞ þ aðaÞf ðtÞ þ f 0 ðtÞ;
for 0 6 t 6 a;
ð2:9Þ
for 0 6 t 6 a:
ð2:10Þ
and respectively
Z 2 ða; tÞ ¼ gðtÞ þ aðaÞf ðtÞ f 0 ðtÞ; Utilizing (1.3) in (2.3) we write:
Z 1 ðx; 0Þ þ Z 2 ðx; 0Þ ¼ 0;
for 0 6 x 6 a:
ð2:11Þ
Finally, using (1.2) and (2.4) in the first equation of (2.3) we obtain:
Z 1 ðx; xÞ ¼
1 ðdqðxÞ þ aðxÞrðxÞÞ; 2
for 0 6 x 6 a:
ð2:12Þ
Note that (2.9)–(2.11) do not lead to a value for aðaÞ. They rather show a compatibility relationship, for if we make t ¼ 0 in (2.9) and (2.10), and add them together we obtain
Z 1 ða; 0Þ þ Z 2 ða; 0Þ ¼ 2ðgð0Þ þ aðaÞf ð0ÞÞ; which on grounds of (2.11) with x ¼ a leads to
gð0Þ þ aðaÞf ð0Þ ¼ 0:
ð2:13Þ
The hypothesis f ð0Þ ¼ 0 ¼ gð0Þ (see (1.1)–(1.5)) makes (2.13) true for any value of aðaÞ. We do need a value for aðxÞ either at x ¼ 0, or else at x ¼ a in order to solve uniquely (2.6). We shall prescribe a final condition, aðaÞ. For a discussion on this matter see Subsection 6.1.2 in [10]. Also note that Eqs. (2.7) and (2.8) clearly show the necessity of a mesh that consists of characteristic curves for the hyperbolic PDE (1.1). Specifically, the mesh lines are:
x t ¼ c; where c ¼ x1 ; x3 ; x5 ; . . . ; xN ; x þ t ¼ c; where c ¼ x1 ; x3 ; x5 ; . . . ; xN ; xNþ2 ; . . . ; x2N1 :
ð2:14Þ
Here 0 ¼ x1 < x2 < x3 < < xN ¼ a < xNþ1 < < x2N1 ¼ 2a are equally spaced points, and N is chosen to be an odd number. See Fig. 3. For each fixed c > 0, Eq. (2.7) integrates along the characteristic curve t ¼ x þ c that is parallel to the secondary diagonal (i.e. the calculation of Z 1 progresses parallel to the secondary diagonal of the coordinate system Oxt), and Eq. (2.8) integrates along the characteristic curve t ¼ x c that is parallel to the main diagonal (i.e. the calculation of Z 2 progresses parallel to the main diagonal). See Fig. 4. In the following subsections we shall use (2.4) and (2.7)–(2.12) to find numerical values of Z 1 ; Z 2 ; r; dq at the mesh nodes. But first we shall solve (2.6).
128
M.-C. Drignei / Applied Mathematics and Computation 220 (2013) 123–141
Fig. 3. The mesh on the domain D0 .
2.1. The Riccati equation Eq. (2.6) along with a prescribed value of aðaÞ is a final value problem for the unknown function aðxÞ. We can call a MATLAB ode solver (e.g. ODE45) to solve it. Another way is as follows. Let
uðxÞ :¼
Z
x
qðsÞ ds:
ð2:15Þ
0
Then, using the notations uj :¼ uðxj Þ, for j ¼ 1; 2; . . . ; N we have:
(
u1 ¼ 0; because x1 ¼ 0; ðqj þ qjþ1 Þ; for j ¼ 1; 2; . . . N 1ðtrapezoid formulaÞ: ujþ1 uj ¼ dx 2
ð2:16Þ
The elementary trapezoid formula can be found in [[13], pp 153], or in any textbook on Calculus. Next, using (2.15) into (2.6) we write:
ða þ uÞ0 ðxÞ ¼ a2 ðxÞ;
for 0 6 x 6 a:
Letting bðxÞ :¼ ða þ uÞðxÞ, the above ODE becomes: 0
2
b ðxÞ ¼ ðb uÞ ðxÞ;
for 0 6 x 6 a:
ð2:17Þ
The final value bðaÞ ¼ aðaÞ þ uðaÞ ¼ aðaÞ þ uN is known since aðaÞ is prescribed and uN is found from (2.16). Hence, numerical values of bðxÞ at x ¼ x1 ; x2 ; . . . ; xN can be obtained via (2.17):
Fig. 4. How the calculations of Z 1 and Z 2 at the mesh nodes are to be performed.
M.-C. Drignei / Applied Mathematics and Computation 220 (2013) 123–141
8 b ¼ aðaÞ þ uN ; > < N 0 2 bN1 ¼ bN ðdxÞbN ¼ bN ðdxÞðbN uN Þ ðbackward difference schemeÞ; > : 0 2 bj2 ¼ bj ð2dxÞbj1 ¼ bj ð2dxÞðbj1 uj1 Þ ; for j ¼ N; N 1; . . . ; 3 ðcentered difference schemeÞ:
129
ð2:18Þ
Finally, from (2.16) and (2.18), and the definition of bðxÞ the numerical values of aðxÞ at x ¼ x1 ; x2 ; . . . ; xN will be obtained:
aN ¼ aðaÞ ¼ prescribed; aj ¼ bj uj ; for j ¼ N 1; N 2; . . . ; 2; 1:
ð2:19Þ
2.2. The data on the vertical boundary x ¼ xN ¼ a Numerical values for Z 1 and Z 2 at the grid points on the vertical boundary are obtained via (2.9) and (2.10), because aðaÞ is prescribed by us and f ðtÞ; gðtÞ are given functions. For the calculations of f 0 ðt j Þ; j ¼ 1; 2; . . . ; N t we shall use a forward, centered or backward difference scheme derivable from the Taylor’s series expansion:
8 0 f2 f1 f ¼ dt ; > > < 1 f f j ¼ 2; 3; . . . ; Nt 1; fj0 ¼ jþ12dtj1 ; > > : 0 fNt fNt 1 fNt ¼ dt ; where we used the abbreviations fj :¼ f ðt j Þ and fj0 ¼ f 0 ðtj Þ, and 0 ¼ t1 < t2 < . . . < t Nt ¼ a is the partition in t-direction with equally spaced points. Note that dt ¼ 2dx, and so N t ¼ Nþ1 (hence the necessity of N being an odd number). The forward 2 2 and backward difference scheme are of order OðdtÞ, whereas the centered difference scheme is of order Oðdt Þ. The numerical values of r and dq at the grid points on this boundary x ¼ xN ¼ a (i.e. r N :¼ rðxN Þ; dqN :¼ dqðxN Þ) will be calculated in Subsection 2.3. 2.3. The data on the vertical line x ¼ xN1 We integrate (2.7) with c ¼ x2N3 and (2.8) with c ¼ x1 ¼ 0 over ½xN1 ; xN and use the simple trapezoid rule (see Fig. 5). We obtain:
Z 1 ðGÞ Z 1 ðHÞ ¼
dx ðaN1 Z 2 ðHÞ þ aN Z 2 ðGÞ þ dqN1 LðHÞ þ dqN LðGÞÞ; 2
ð2:20Þ
Z 2 ðFÞ Z 2 ðHÞ ¼
dx ðaN1 Z 1 ðHÞ þ aN Z 1 ðFÞ þ dqN1 LðHÞ þ dqN LðFÞÞ: 2
ð2:21Þ
From (2.12) with x ¼ xN1 and respectively x ¼ xN we obtain:
Z 1 ðHÞ ¼
1 ðdqN1 þ aN1 r N1 Þ; 2
ð2:22Þ
Z 1 ðFÞ ¼
1 ðdqN þ aN rN Þ: 2
ð2:23Þ
Fig. 5. Z 1 ; Z 2 ; dq; r for the highest triangular cell between the lines x ¼ xN and x ¼ xN1 .
130
M.-C. Drignei / Applied Mathematics and Computation 220 (2013) 123–141
From (2.4) we have
rðxN Þ rðxN1 Þ ¼
Z
xN
dqðsÞ ds; xN1
which after the application of the simple trapezoid rule we obtain:
rN r N1 ¼
dx ðdqN1 þ dqN Þ: 2
ð2:24Þ
Taking x ¼ xN ¼ a in (2.4) and using (1.2) with x ¼ a and (1.4) with t ¼ a, we obtain:
rN ¼ 2f N :
ð2:25Þ
We used the abbreviations dqi :¼ dqðxi Þ; r i :¼ rðxi Þ; ai :¼ aðxi Þ. We note that (2.20)–(2.25) is a linear system with six equations and six unknowns: Z 1 ðHÞ; Z 2 ðHÞ; dqN1 ; rN1 ; dqN ; r N , because Z 1 ðGÞ; Z 2 ðGÞ; Z 1 ðFÞ; Z 2 ðFÞ are known from the previous step (see Subsection 2.2). We can solve this system easily. Next, we look at an arbitrary triangular cell between the vertical lines x ¼ xN and x ¼ xN1 . See Fig. 6. For each 2 6 i 6 N1 , 2 we integrate (2.7) with c ¼ x2N2i1 and (2.8) with c ¼ x2i1 over ½xN1 ; xN and use the simple trapezoid rule. We obtain respectively:
Z 1 ðBÞ Z 1 ðCÞ ¼
dx ðaN1 Z 2 ðCÞ þ aN Z 2 ðBÞ þ dqN1 LðCÞ þ dqN LðBÞÞ; 2
ð2:26Þ
Z 2 ðAÞ Z 2 ðCÞ ¼
dx ðaN1 Z 1 ðCÞ þ aN Z 1 ðAÞ þ dqN1 LðCÞ þ dqN LðAÞÞ: 2
ð2:27Þ
We solve the linear system (2.26) - (2.27) for the unknowns Z 1 ðCÞ; Z 2 ðCÞ, because Z 1 ðAÞ; Z 2 ðAÞ; Z 1 ðBÞ; Z 2 ðBÞ were obtained in Subsection 2.2 and dqN1 ; dqN were obtained from (2.20)–(2.25). 2.4. The data on the vertical line x ¼ x2j1 We shall obtain numerical values of Z 1 ; Z 2 ; r; dq at the grid points located on the vertical line x ¼ x2j1 from the values of Z 1 ; Z 2 ; r; dq at the grid points on the vertical line x ¼ x2j we calculated previously. We start with the top triangular cell spanned between the vertical lines x ¼ x2j and x ¼ x2j1 . See Fig. 7. We integrate (2.7) with c ¼ x4j3 , and (2.8) with c ¼ x1 ¼ 0 over ½x2j1 ; x2j and use the simple trapezoid rule. We obtain respectively:
Z 1 ðPÞ Z 1 ðMÞ ¼
dx ða2j1 Z 2 ðMÞ þ a2j Z 2 ðPÞ þ dq2j1 LðMÞ þ dq2j LðPÞÞ; 2
ð2:28Þ
Z 2 ðKÞ Z 2 ðMÞ ¼
dx ða2j1 Z 1 ðMÞ þ a2j Z 1 ðKÞ þ dq2j1 LðMÞ þ dq2j LðKÞÞ: 2
ð2:29Þ
From (2.12) with x ¼ x2j1 we write:
Z 1 ðMÞ ¼
1 ðdq2j1 þ a2j1 r2j1 Þ: 2
From (2.4) we have
Fig. 6. Z 1 ; Z 2 for an arbitrary triangular cell between the lines x ¼ xN and x ¼ xN1 .
ð2:30Þ
M.-C. Drignei / Applied Mathematics and Computation 220 (2013) 123–141
131
Fig. 7. Z 1 ; Z 2 ; dq; r for the highest triangular cell between the lines x ¼ x2j and x ¼ x2j1 .
rðx2j Þ rðx2j1 Þ ¼
Z
x2j
dqðsÞ ds;
x2j1
which after the application of the simple trapezoid rule becomes:
r2j r2j1 ¼
dx ðdq2j1 þ dq2j Þ: 2
ð2:31Þ
We solve the linear system (2.28)–(2.31) for the unknowns Z 1 ðMÞ; Z 2 ðMÞ; r 2j1 ; dq2j1 . Next, we look at an arbitrary cell that spans between the vertical lines x ¼ x2j and x ¼ x2j1 . See Fig. 8. For each 2 6 i 6 j 1, we integrate (2.7) with c ¼ x4j2i1 , and (2.8) with c ¼ x2i1 over ½x2j1 ; x2j and use the simple trapezoid rule. We obtain respectively:
Z 1 ðSÞ Z 1 ðTÞ ¼
dx ða2j1 Z 2 ðTÞ þ a2j Z 2 ðSÞ þ dq2j1 LðTÞ þ dq2j LðSÞÞ; 2
ð2:32Þ
Z 2 ðRÞ Z 2 ðTÞ ¼
dx ða2j1 Z 1 ðTÞ þ a2j Z 1 ðRÞ þ dq2j1 LðTÞ þ dq2j LðRÞÞ: 2
ð2:33Þ
We solve the linear system (2.32) - (2.33) for the unknowns Z 1 ðTÞ; Z 2 ðTÞ. Note that dq2j1 was obtained from (2.28)–(2.31), and dq2j from the calculations done for the previous vertical line x ¼ x2j . Finally, we calculate Z 1 and Z 2 at the lowest grid point on the vertical line x ¼ x2j1 . See Fig. 9. From (2.11) with x ¼ x2j1 we write:
Z 1 ðUÞ þ Z 2 ðUÞ ¼ 0:
ð2:34Þ
Taking c ¼ x2j1 in (2.8) and integrating over ½x2j1 ; x2j we obtain:
Fig. 8. Z 1 ; Z 2 for an arbitrary triangular cell between the lines x ¼ x2j and x ¼ x2j1 .
132
M.-C. Drignei / Applied Mathematics and Computation 220 (2013) 123–141
Fig. 9. Z 1 ; Z 2 at the lowest grid point on the vertical line x ¼ x2j1 .
Z 2 ðVÞ Z 2 ðUÞ ¼
dx ða2j1 Z 1 ðUÞ þ a2j Z 1 ðVÞ þ dq2j1 LðUÞ þ dq2j LðVÞÞ: 2
ð2:35Þ
We solve the linear system (2.34) - (2.35) for Z2(U) for Z 1 ðUÞ; Z 2 ðUÞ. 2.5. The data on the vertical line x ¼ x2j2 We shall perform our calculations from the highest triangular cell to the lowest triangular cell spanned between the vertical lines x ¼ x2j1 and x ¼ x2j2 . The highest cell is illustrated in Fig. 10. Integrating (2.7) with c ¼ x4j5 , and (2.8) with c ¼ x1 ¼ 0 over ½x2j2 ; x2j1 and using the simple trapezoid rule, we obtain respectively:
Z 1 ðZÞ Z 1 ðXÞ ¼
dx ða2j2 Z 2 ðXÞ þ a2j1 Z 2 ðZÞ þ dq2j2 LðXÞ þ dq2j1 LðZÞÞ; 2
ð2:36Þ
Z 2 ðYÞ Z 2 ðXÞ ¼
dx ða2j2 Z 1 ðXÞ þ a2j1 Z 1 ðYÞ þ dq2j2 LðXÞ þ dq2j1 LðYÞÞ: 2
ð2:37Þ
From (2.12) with x ¼ x2j2 we obtain:
Z 1 ðXÞ ¼
1 ðdq2j2 þ a2j2 r 2j2 Þ: 2
From (2.4) we have:
rðx2j1 Þ rðx2j2 Þ ¼
Z
x2j1
dqðsÞ ds;
x2j2
which after applying the simple trapezoid rule gives:
Fig. 10. Z 1 ; Z 2 ; dq; r for the highest triangular cell between the lines x ¼ x2j1 and x ¼ x2j2 .
ð2:38Þ
M.-C. Drignei / Applied Mathematics and Computation 220 (2013) 123–141
r2j1 r 2j2 ¼
dx ðdq2j2 þ dq2j1 Þ: 2
133
ð2:39Þ
We solve the linear system (2.36)–(2.39) for the unknowns Z 1 ðXÞ; Z 2 ðXÞ; r 2j2 and dq2j2 . Note that Z 1 ðZÞ; Z 2 ðZÞ; Z 1 ðYÞ; Z 2 ðYÞ; r 2j1 ; dq2j1 were obtained in Subsection 2.4. Next we look at an arbitrary triangular cell spanned between the vertical lines x ¼ x2j1 and x ¼ x2j2 . See Fig. 11. For each 2 6 i 6 j 1, we integrate (2.7) with c ¼ x4j2i3 , and (2.8) with c ¼ x2i1 over ½x2j2 ; x2j1 . We obtain:
Z 1 ðEÞ Z 1 ðIÞ ¼
dx ða2j2 Z 2 ðIÞ þ a2j1 Z 2 ðEÞ þ dq2j2 LðIÞ þ dq2j1 LðEÞÞ; 2
ð2:40Þ
Z 2 ðDÞ Z 2 ðIÞ ¼
dx ða2j2 Z 1 ðIÞ þ a2j1 Z 1 ðDÞ þ dq2j2 LðIÞ þ dq2j1 LðDÞÞ: 2
ð2:41Þ
We solve the linear system (2.40) - (2.41) for the unknowns Z 1 ðIÞ; Z 2 ðIÞ. Note that dq2j2 was obtained from (2.36)–(2.39)), and Z 1 ðDÞ; Z 2 ðDÞ; Z 1 ðEÞ; Z 2 ðEÞ; dq2j1 were obtained in Subsection 2.4. 2.6. The data on the vertical line x ¼ x1 There is only one grid point on the vertical line x ¼ x1 , namely O – the origin of the coordinate system. See Fig. 12. Taking x ¼ x1 ¼ 0 in (2.11) we have:
Z 1 ðOÞ þ Z 2 ðOÞ ¼ 0:
ð2:42Þ
Integrating (2.8) with c ¼ x1 ¼ 0 over ½x1 ; x2 and using the simple trapezoid formula, we obtain:
Z 2 ðQ Þ Z 2 ðOÞ ¼
dx ða1 Z 1 ðOÞ þ a2 Z 1 ðQ Þ þ dq1 LðOÞ þ dq2 LðQ ÞÞ: 2
ð2:43Þ
From (2.12) with x ¼ x1 ¼ 0 we have:
Z 1 ðOÞ ¼
1 ðdq1 þ a1 r 1 Þ: 2
ð2:44Þ
Due to (2.4) we write:
rðx2 Þ rðx1 Þ ¼
Z
x2
dqðsÞ ds; x1
which after applying the simple trapezoid rule gives:
r2 r1 ¼
dx ðdq1 þ dq2 Þ: 2
We solve the linear system (2.42)–(2.45) for the unknowns Z 1 ðOÞ; Z 2 ðOÞ; r 1 and dq1 . 2.7. The calculation of W; W x ; W t at the mesh nodes From (2.3) we can obtain W t ðx; tÞ and W x ðx; tÞ þ aðxÞWðx; tÞ at all mesh nodes ðx; tÞ:
Fig. 11. Z 1 ; Z 2 for an arbitrary triangular cell between the lines x ¼ x2j1 and x ¼ x2j2 .
ð2:45Þ
134
M.-C. Drignei / Applied Mathematics and Computation 220 (2013) 123–141
Fig. 12. Z 1 ; Z 2 ; dq; r at the lowest grid point O.
W t ðx; tÞ ¼
Z 1 ðx; tÞ Z 2 ðx; tÞ ; 2
W x ðx; tÞ þ aðxÞWðx; tÞ ¼
Z 1 ðx; tÞ þ Z 2 ðx; tÞ : 2
ð2:46Þ
ð2:47Þ
Since ri :¼ rðxi Þ, for i ¼ 1; 2; . . . ; N were calculated in Subsections 2.3 - 2.6 we will also be able to calculate Wðx; xÞ via (1.2) and (2.4) at all mesh nodes on the slanted boundary. Letting Wðk; lÞ :¼ Wðxk ; xl Þ; W x ðk; lÞ :¼ W x ðxk ; xl Þ; W t ðk; lÞ :¼ W t ðxk ; xl Þ and using (2.46), (2.47), one way to obtain W t ; W; W x at all mesh nodes is as follows. For each k ¼ N; N 1; . . . ; 6; 5, perform calculations on the vertical line x ¼ xk going from top to bottom:
8 Z 1 ðk;kÞZ 2 ðk;kÞ ; > 2 < W t ðk; kÞ ¼ 1 Wðk; kÞ ¼ 2 r k ; > : 2 ðk;kÞ W x ðk; kÞ ¼ Z1 ðk;kÞþZ ak Wðk; kÞ; 2 8 Z 1 ðk;k2ÞZ 2 ðk;k2Þ ; > 2 < W t ðk; k 2Þ ¼ Wðk; k 2Þ ¼ Wðk; kÞ ð2dxÞW t ðk; kÞ; ðbackward difference schemeÞ; > : 2 ðk;k2Þ ak Wðk; k 2Þ; W x ðk; k 2Þ ¼ Z1 ðk;k2ÞþZ 2 and for each i ¼ k; k 2; k 4; . . . ; 7; 5 (if k is odd), or i ¼ k; k 2; k 4; . . . ; 8; 6 (if k is even)
8 Z 1 ðk;i4ÞZ 2 ðk;i4Þ ; > 2 < W t ðk; i 4Þ ¼ Wðk; i 4Þ ¼ Wðk; iÞ ð4dxÞW t ðk; i 2Þ; ðcentered difference schemeÞ; > : 2 ðk;i4Þ ak Wðk; i 4Þ: W x ðk; i 4Þ ¼ Z1 ðk;i4ÞþZ 2 For k ¼ 4 and k ¼ 3, perform calculations on the vertical line x ¼ xk going from top to bottom:
8 Z 1 ðk;kÞZ 2 ðk;kÞ ; > 2 < W t ðk; kÞ ¼ 1 Wðk; kÞ ¼ 2 r k ; > : 2 ðk;kÞ W x ðk; kÞ ¼ Z1 ðk;kÞþZ ak Wðk; kÞ; 2 8 Z 1 ðk;k2ÞZ 2 ðk;k2Þ ; > 2 < W t ðk; k 2Þ ¼ Wðk; k 2Þ ¼ Wðk; kÞ ð2dxÞW t ðk; kÞ; ðbackward difference schemeÞ; > : 2 ðk;k2Þ ak Wðk; k 2Þ; W x ðk; k 2Þ ¼ Z1 ðk;k2ÞþZ 2 For k ¼ 2, set
8 Z 1 ðk;kÞZ 2 ðk;kÞ ; > 2 < W t ðk; kÞ ¼ 1 Wðk; kÞ ¼ 2 r k ; > : 2 ðk;kÞ W x ðk; kÞ ¼ Z1 ðk;kÞþZ ak Wðk; kÞ; 2 For k ¼ 1, set
M.-C. Drignei / Applied Mathematics and Computation 220 (2013) 123–141
135
8 Z 1 ðk;kÞZ 2 ðk;kÞ > ; < W t ðk; kÞ ¼ 2 Wðk; kÞ ¼ 0; ðdue to ð1:3Þ with x ¼ x1 ¼ 0Þ; > : W x ðk; kÞ ¼ 0; ðdue to ð1:3Þ with x ¼ x1 ¼ 0Þ: 3. The algorithm Let N be an odd number. Let 0 ¼ x1 < x2 < < xN ¼ a < xNþ1 < xNþ2 < < x2N1 ¼ 2a be equally spaced points. Place the mesh (2.14) on the domain D0 . Find Z 1 and Z 2 at the grid points on the vertical boundary x ¼ xN ¼ a. See Subsection 2.2. Find Z 1 ; Z 2 at the grid points on the vertical line x ¼ xN1 and rN ; dqN ; rN1 ; dqN1 . See Subsection 2.3. For each j ¼ N1 ; N3 ; . . . ; 4; 3; 2 2 2 – find Z 1 ; Z 2 ; dq; r at the grid points on the vertical line x ¼ x2j1 (see Subsection 2.4), – find Z 1 ; Z 2 ; dq; r at the grid points on the vertical line x ¼ x2j2 (see Subsection 2.5). Find Z 1 ; Z 2 ; dq1 ; r1 at the grid points on the vertical line x ¼ x1 ¼ 0. (See Subsection 2.6), Find W; W x ; W t at the mesh nodes. (See Subsection 2.7). 4. Remarks
Remark 1. The ODE (2.8) with c ¼ x1 ¼ 0 can be further manipulated due to (2.12), (2.4) and (2.6):
d 1 Z 2 ðx; xÞ ¼ aðxÞZ 1 ðx; xÞ þ dqðxÞLðx; xÞ ¼ aðxÞ ðdqðxÞ þ aðxÞrðxÞÞ þ dqðxÞLðx; xÞ dx 2 1 1 0 2 ¼ ðaðxÞr ðxÞ þ a ðxÞrðxÞÞ þ dqðxÞLðx; xÞ ¼ ðaðxÞr 0 ðxÞ þ a0 ðxÞrðxÞ þ qðxÞrðxÞÞ þ dqðxÞLðx; xÞ; 2 2
for 0 6 x 6 a;
which is equivalent to
d 1 1 Z 2 ðx; xÞ aðxÞrðxÞ ¼ qðxÞrðxÞ þ dqðxÞLðx; xÞ; dx 2 2
for 0 6 x 6 a:
ð4:1Þ
Integrating (4.1) over appropriate intervals and using the simple trapezoid formula we produce alternative equations to (2.21), (2.29), (2.37), and (2.43), respectively:
1 1 dx dx Z 2 ðFÞ Z 2 ðHÞ ðaN r N aN1 r N1 Þ ¼ ðqN1 r N1 þ qN rN Þ þ ðdqN1 LðHÞ þ dqN LðFÞÞ; 2 2 2 2
ð4:2Þ
1 1 dx dx Z 2 ðKÞ Z 2 ðMÞ ða2j r 2j a2j1 r 2j1 Þ ¼ ðq2j1 r 2j1 þ q2j r 2j Þ þ ðdq2j1 LðMÞ þ dq2j LðKÞÞ; 2 2 2 2
ð4:3Þ
1 1 dx dx Z 2 ðYÞ Z 2 ðXÞ ða2j1 r2j1 a2j2 r2j2 Þ ¼ ðq2j2 r 2j2 þ q2j1 r2j1 Þ þ ðdq2j2 LðXÞ þ dq2j1 LðYÞÞ; 2 2 2 2
ð4:4Þ
1 1 dx dx Z 2 ðQ Þ Z 2 ðOÞ ða2 r 2 a1 r1 Þ ¼ ðq1 r1 þ q2 r2 Þ þ ðdq1 LðOÞ þ dq2 LðQ ÞÞ: 2 2 2 2
ð4:5Þ
There is however not a significant difference in the algorithm’s accuracy when using formulas (4.2)–(4.5) in place of (2.21), (2.29), (2.37), (2.43), respectively.
5. Numerical examples We illustrate the effectiveness of our numerical method/algorithm on the following two examples where analytical solutions are available, thus making the comparison with the numerical solutions possible. To obtain these examples, we tackled the boundary value problem (1.1)–(1.5) backwards: we rather started with a given formula for Wðx; tÞ and produced dqðxÞ; qðxÞ; Lðx; tÞ; f ðtÞ; gðtÞ such that (1.1)–(1.5) are met. qðxÞ t ð1þ2x2 2 Þ 2 2 2 2 Example 1. For Wðx; tÞ ¼ ex t, we obtain dqðxÞ ¼ 2ex ð1 þ 2x2 Þ; Lðx; tÞ ¼ ; f ðtÞ ¼ ea t; gðtÞ ¼ 2aea t. Thus, the 1þ2x2 2 hypothesis condition f ð0Þ ¼ 0 ¼ gð0Þ is satisfied. We chose qðxÞ ¼ 1þx3 and a ¼ 1, for the numerical illustration.
136
M.-C. Drignei / Applied Mathematics and Computation 220 (2013) 123–141
Fig. 13. Example 1 – the comparison for dqðxÞ, when N ¼ 15 x-partition points have been chosen.
Fig. 14. Example 1 – illustration of
Fig. 15. Example 1 – illustration of
W exact ðx;tÞW numeric ðx;tÞ , kW exact k2
when N ¼ 15 x-partition points have been chosen.
W exact x ðx;tÞW numeric ðx;tÞ x Þ, k2 kW exact x
when N ¼ 15 x-partition points have been chosen.
½1þqðxÞðsin xþcos xÞt Example 2. For Wðx; tÞ ¼ t ðsin x þ cos xÞ, we obtain dqðxÞ ¼ 2½ð1 xÞ sin x þ ð1 þ xÞ cos x; Lðx; tÞ ¼ 2½ð1xÞ ; f ðtÞ ¼ sin xþð1þxÞ cos x ðsin a þ cos aÞ t; gðtÞ ¼ ðcos a sin aÞ t, so that f ð0Þ ¼ 0 ¼ gð0Þ, as required. For the numerical illustration we took a ¼ 1 and 1 qðxÞ ¼ sin xþcos . x MATLAB figure-outputs are provided for our two examples with N ¼ 15 partition points in x-direction. See Figs. 13–20, where the exact solution dqðxÞ (in blue) and the numerical solution dqðxÞ (in red) are displayed for comparison, as well as the divided differences:
W exact ðx; tÞ W numeric ðx; tÞ ; kW exact k2
M.-C. Drignei / Applied Mathematics and Computation 220 (2013) 123–141
Fig. 16. Example 1 – illustration of
W exact ðx;tÞW numeric ðx;tÞ t t , k2 kW exact t
when N ¼ 15 x-partition points have been chosen.
Fig. 17. Example 2 – the comparison for dqðxÞ, when N ¼ 15 x-partition points have been chosen.
Fig. 18. Example 2 – illustration of
W exact ðx; tÞ W numeric ðx; tÞ x x ; k2 kW exact x W exact ðx; tÞ W numeric ðx; tÞ t t : k kW exact 2 t
W exact ðx;tÞW numeric ðx;tÞ , kW exact k2
when N ¼ 15 x-partition points have been chosen.
137
138
M.-C. Drignei / Applied Mathematics and Computation 220 (2013) 123–141
Fig. 19. Example 2 – illustration of
Fig. 20. Example 2 – illustration of
W exact ðx;tÞW numeric ðx;tÞ x x Þ, kW exact k2 x
W exact ðx;tÞW numeric ðx;tÞ t t , k2 kW exact t
when N ¼ 15 x-partition points have been chosen.
when N ¼ 15 x-partition points have been chosen.
The reason we chose to display the derivatives W x ðx; tÞ and W t ðx; tÞ in addition to the solution fdqðxÞ; Wðx; tÞg to (1.1)–(1.5) is because they are directly involved in the computation of Wðx; tÞ (see Subsection 2.7).
6. Conclusions We solved a boundary value problem on a triangular domain for two unknown functions. The numerical method consisted of placing a mesh on the triangular domain with lines that are the characteristic curves of the hyperbolic PDE in the boundary value problem. Then a change of dependent variables was made and the numerical values for the needed unknown functions were obtained at the mesh nodes by going from one vertical line to another, starting from the vertical boundary. Integration was performed along the characteristic curves of the original PDE, and the simple trapezoid formula was utilized to carry out these calculations. Numerical examples were provided to illustrate the effectiveness of the method/ algorithm. This algorithm can be further embedded in a Newton type algorithm through which the numerical solution to an inverse Sturm–Liouville problem can be obtained, provided two spectral data are known (i.e. an inverse two-spectra problem can be solved numerically). Appendix A In this section we present the connection with the inverse two-spectra problems for the canonical Sturm–Liouville operator. The canonical Sturm–Liouville equation on a finite interval ½0; a is
u00 ðxÞ þ qðxÞuðxÞ ¼ kuðxÞ;
0 6 x 6 a;
ð7:1Þ
where qðxÞ is the coefficient function (also called the potential function of the canonical Sturm–Liouville operator). Two boundary conditions are associated to (7.1), one per end-point. The type of boundary conditions gives the name of the boundary value problem so formed. To fix the ideas, let us consider the following two situations:
139
M.-C. Drignei / Applied Mathematics and Computation 220 (2013) 123–141
uð0Þ ¼ 0 ¼ uðaÞ
ð7:2Þ
uð0Þ ¼ 0 ¼ u0 ðaÞ:
ð7:3Þ
and
Thus, (7.1) & (7.2) is called the Dirichlet eigenvalue problem for the canonical Sturm–Liouville operator, because both boundary conditions are Dirichlet type, and (7.1) & (7.3) is refereed to as the Dirichlet–Neumann eigenvalue problem for the canonical Sturm–Liouville operator, because the boundary condition at the left end point x ¼ 0 is Dirichlet type, and the boundary condition at the right end point x ¼ a is Neumann type. The direct problem is formulated as follows: given qðxÞ in a certain space of functions (e.g. L2 ð0; aÞ; H1 ð0; aÞ), find the eigenpairs fk; uðxÞg of either one of the BVP mentioned above, or any other BVP that consists of (7.1) and two boundary conditions (one per end-point). The inverse problem formulates as follows: given information about the spectral data, find qðxÞ. Note that the spectral data do not mean the eigenpairs, because if so, then the problem would become trivial (i.e. even with the knowledge of only one eigenpair you can get qðxÞ from (7.1) as if you would solve a linear equation). Spectral data mean for instance, one sequence of eigenvalues and one sequence of norming constants. Or they can be two sequences of eigenvalues. We discuss the two-spectra case. This means, we know two sequences of eigenvalues. And since we decided to focus on the boundary conditions (7.2) and (7.3), we shall go ahead with the sequence of Dirichlet eigenvalues fkn gnP1 and the sequence of Dirichlet–Neumann eigenvalues fmn gnP1 . That means, fkn ; un ðxÞgnP1 satisfy (7.1) & (7.2), and fmn ; v n ðxÞgnP1 satisfy (7.1) & (7.3), for the same potential function qðxÞ. The facts that the eigenvalues are real (provided qðxÞ is a real-valued function), countably many, and satisfy specific asymptotic formulas were rigorously established in the literature devoted to the direct Sturm–Liouville problems. See for example [[20] Section 4.3]:
kn ¼
np2 þ ½q þ cn ; a
for n ¼ 1; 2; 3; . . . ;
where fcn gnP1 is an l2 sequence of real numbers (i.e.
mn
2 ðn 12Þp ¼ þ ½q þ dn ; a
ð7:4Þ P1
2
n¼1 ðcn Þ
< 1), and
for n ¼ 1; 2; 3; . . . ;
ð7:5Þ
Ra with fdn gnP1 2 l2 . Here ½q :¼ 1a 0 qðxÞ dx is the mean of the function q 2 L2R ð0; aÞ. Also, one can prove that k is a Dirichlet eigenvalue of the canonical Sturm–Liouville operator with potential function qðxÞ, if and only if k is a zero of the function
k 2 C ! Sða; q; kÞ;
ð7:6Þ
called the characteristic function of the BVP (7.1) & (7.2), and m is a Dirichlet–Neumann eigenvalue of the canonical Sturm– Liouville operator with potential function qðxÞ, if and only if m is a zero of the function
k 2 C ! S0 ða; q; kÞ;
ð7:7Þ 2
called the characteristic function of the BVP (7.1) & (7.3). Here, Sðx; q; kÞ stands for the weak solution (i.e. in H ð0; aÞ) of the initial value problem
8 00 > < u ðxÞ þ qðxÞuðxÞ ¼ kuðxÞ; uð0Þ ¼ 0; > : 0 u ð0Þ ¼ 1:
0 6 x 6 a; ð7:8Þ
pffiffi The notation SðÞ was chosen to resemble the function sine, specifically sinðpffiffikkxÞ which is the weak solution to (7.8) with qðxÞ ¼ 0. We choose to detail for the second case, the equivalence between being an eigenvalue and being a zero of the characteristic function. For the other case, the reasons are similar. So, let us assume that m is a Dirichlet–Neumann eigenvalue. Let v ðxÞ be an associated eigenfunction. That is, fm ; v ðxÞg satisfies (7.1) & (7.3). It follows that ðv Þ0 ð0Þ – 0, as otherwise (7.1) and the first boundary condition of (7.3) would imply that v ðxÞ ¼ 0, for all x 2 ½0; a, because the initial value problem (7.1) & uð0Þ ¼ 0 ¼ u0 ð0Þ has only one solution, and the identically zero function is a solution. We cannot have v ðxÞ ¼ 0, for all ~ ðxÞ :¼ ðvv ÞðxÞ ~ ðxÞ satisfies (7.8) x 2 ½0; a, because v ðxÞ is an eigenfunction. Thus, v is well-defined. It follows immediately that v 0 ð0Þ 0 ðaÞ 0 ~ ðxÞ ¼ Sðx; q; m Þ. Hence, S ða; q; m Þ ¼ v ~ 0 ðaÞ ¼ ððvv ÞÞ0 ð0Þ with k ¼ m , and therefore v ¼ 0. We showed that a Dirichlet–Neumann eigenvalue m is a zero of (7.7). Conversely, let m be a zero of (7.7). Then, due to the definition of Sðx; q; kÞ with k ¼ m , we have that fm ; Sðx; q; m Þg satisfies (7.1) & (7.3). This means nothing but m is a Dirichlet–Neumann eigenvalue of the canonical Sturm–Liouville operator with potential function qðxÞ. Hence, we have
Sða; q; kn Þ ¼ 0;
for n ¼ 1; 2; 3; . . . ;
ð7:9Þ
S0 ða; q; mn Þ ¼ 0;
for n ¼ 1; 2; 3; . . . :
ð7:10Þ
and
140
M.-C. Drignei / Applied Mathematics and Computation 220 (2013) 123–141
It is also known from the theory of direct Sturm–Liouville problems that Sðx; q; kÞ has an integral representation (see for instance [[20] Example 4.19 and Theorem 4.18]):
Sðx; q; kÞ ¼
pffiffiffi pffiffiffi Z x sinð kxÞ sinð ktÞ pffiffiffi þ dt; Kðx; t; qÞ pffiffiffi k k 0
for 0 6 x 6 a:
ð7:11Þ
Here Kðx; t; qÞ goes under the name of Gelfand-Levitan-Marchenko kernel. As seen in (7.11), the solution to one initial value problem, namely (7.8) with qðxÞ ¼ 0, is connected to the solution to another initial value problem (7.8) by means of a second kind Volterra integral operator with kernel Kðx; t; qÞ. The kernel satisfies the Goursat problem (1.8) in the triangle D0 :¼ fðx; tÞ 2 R2 : 0 6 t 6 x 6 ag. From (7.11) we obtain:
pffiffiffi pffiffiffi Z x pffiffiffi sinð kxÞ sinð ktÞ þ dt; S0 ðx; q; kÞ ¼ cosð kxÞ þ Kðx; x; qÞ pffiffiffi Kx ðx; t; qÞ pffiffiffi k k 0
for 0 6 x 6 a:
ð7:12Þ
From (7.9), (7.10), (7.11), (7.12), and the first boundary condition of (1.8) we obtain that Kða; t; qÞ is the solution ~f ðtÞ to
Z pffiffiffiffiffi sinð kn aÞ þ
a
pffiffiffiffiffi ~f ðtÞ sinð kn tÞ dt ¼ 0;
for n ¼ 1; 2; 3; . . . ;
ð7:13Þ
0
and Kx ða; t; qÞ is the solution g~ðtÞ to
pffiffiffiffiffi
pffiffiffiffiffi
mn cosð mn aÞ þ
pffiffiffiffiffi a½q sinð mn aÞ þ 2
Z
a
g~ðtÞ sinð
pffiffiffiffiffi
mn tÞ dt ¼ 0; for n ¼ 1; 2; 3; . . . :
ð7:14Þ
0
Note that in (7.14), ½q - the mean of the function qðxÞ can be estimated from the eigenvalues. Due to (7.4) and (7.5), we can write
kN1
2 2 ðN2 12Þp N1 p ½q mN2 ; a a
for N 1 ; N 2 large enough. Hence, everything in (7.13) and (7.14) depends only on the known quantities fkn gnP1 and fmn gnP1 . After solving the systems of integral equations (7.13) and (7.14), we can proceed to the determination of the function qðxÞ – the unknown in the inverse two-spectra problem, by applying the Newton method to the non-linear equation (1.6), where the map FðqÞ is defined by (1.7). ~ðtÞ with Note that, numerically we can solve (7.13) and (7.14) by writing a Fourier series expansion for each of ~f ðtÞ and g 1 2 appropriate bases of functions. Since in practice only the first several eigenvalues are available, say fkn gNn¼1 and fmn gNn¼1 , we truncate these Fourier series after the first N 1 , and respectively N 2 terms. Thus, we will need to keep the equations that correspond to n ¼ 1; 2; . . . N 1 from (7.13), and the equations that correspond to n ¼ 1; 2; . . . N 2 from (7.14). The two linear systems so obtained will be square N 1 N 1 and respectively N 2 N 2 linear systems for the first N 1 and N 2 Fourier coefficients of ~f ðtÞ and g~ðtÞ, respectively. Once these coefficients are calculated, we have estimates of ~f ðtÞ and g~ðtÞ and we can continue as outlined above. We close this section by saying that a canonical Sturm–Liouville equation models the infinitesimal, vertical vibration of an elastic string of negligible mass, and the boundary conditions associated to the equation are derivable from the constraints the ends of the string are subjected to (see for example [[10] Chapter 2 and Section 5.2]). The reader may find of interest the papers [5–8,11,14,21,24,27], which are either on inverse problems or on numerical partial differential equations. Appendix B. Supplementary data Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/ j.amc.2013.05.041. References [1] R.T. Dames, Stability and convergence for a numerical solution of the Goursat problem, Journal of Mathematical Physics 38 (1959–1960) 42–67. [2] J.T. Day, A Gaussian quadrature method for the numerical solution of the characteristic initial value problem uxy ¼ f ðx; y; uÞ, Mathematics of Computation 17 (1963) 438–441. [3] J.T. Day, A Runge-Kutta method for the numerical solution of the Goursat problem in hyperbolic partial differential equations, The Computer Journal 9 (1966) 81–83. [4] J.T. Day, On the numerical solution of the Goursat problem, Nordisk Tidskrift for Informationsbehandling (BIT) 11 (1971) 450–454. [5] M. Dehghan, On the solution of an initial-boundary value problem that combines Neumann and integral condition for the wave equation, Numerical Methods for Partial Differential Eqs. 21 (2005) 24–40. [6] M. Dehghan, A. Ghesmati, Combination of meshless local weak and strong (MLWS) forms to solve the two dimensional hyperbolic telegraph equation, Engineering Analysis with Boundary Elements 34 (2010) 324–336. [7] M. Dehghan, A. Mohebbi, High order implicit collocation method for the solution of two-dimensional linear hyperbolic equation, Numerical Methods for Partial Differential E qs. 25 (2009) 232–243. [8] M. Dehghan, R. Salehi, A method based on meshless approach for the numerical solution of the two-space dimensional hyperbolic telegraph equation, Mathematical Methods in the Applied Science 35 (2012) 1220–1233.
M.-C. Drignei / Applied Mathematics and Computation 220 (2013) 123–141
141
[9] J. Diaz, On an analogue of the Euler–Cauchy polygon method for the numerical solution of uxy ¼ f ðx; y; u; ux ; uy Þ, Archive for Rational Mechanics and Analysis 1 (1957) 154–180. [10] M.C. Drignei, Inverse Sturm–Liouville problems using multiple spectra, Ph.D Thesis, Iowa State University, 2008. [11] M.C. Drignei, Constructibility of an L2R ð0; aÞ solution to an inverse Sturm–Liouville problem using three Dirichlet spectra, Inverse Problems 26 (2010) 025003 (29pp). [12] D.J. Evans, B.B. Sanugi, Numerical solution of the Goursat problem by a nonlinear trapezoidal formula, Applied Mathematics Letters 1 (3) (1988) 221– 223. [13] W. Gautschi, Numerical Analysis an Introduction, Birkhauser, Boston-Besel-Berlin, 1997. [14] K. Gou, B. Sun, Numerical solution of the Goursat problem on a triangular domain with mixed boundary conditions, Applied Mathematics and Computation 217 (2011) 8765–8777. [15] N. Gradinaru, G. Teodoru, On certain numerical methods for the Goursat problem, (Romanian) Buletinul Institutului Politehnic din Iasi, Serie Noua 20 (24) (1974) 43–46. [16] E. Goursat, A Course in Mathematical Analysis, vol. 3, Dover, New York, 2006. [17] M.K. Jain, K.D. Sharma, Cubature method for the numerical solution of the characteristic initial value problem, Journal of the Australian Mathematical Society 8 (1968) 355–368. [18] O. Jokhadze, Cauchy–Goursat problem for one-dimensional semilinear wave equations, Communications in Partial Differential Equations 34 (4–6) (2009) 367–382. [19] O. Jokhadze, On existence and nonexistence of global solutions of Cauchy–Goursat problem for nonlinear wave equations, Journal of Mathematical Analysis and Applications 2 (2008) 1033–1045. [20] A. Kirsch, An introduction to the mathematical theory of inverse problems, Applied Mathematical Sciences 120, Springer-Verlag, New York, 1996. [21] M. Lakestani, M. Dehghan, The use of Chebyshev cardinal functions for the solution of a partial differential equation with an unknown time-dependent coefficient subject to an extra measurement, Journal of Computational and Applied Mathematics 235 (2010) 669–678. [22] R.H. Moore, A Runge-Kutta procedure for the Goursat problem in hyperbolic partial differential equations, Archive for Rational Mechanics and Analysis 7 (1961) 37–63. [23] M.A.S. Nasir, A.I.M. Ismail, Numerical solution of a linear Goursat problem: stability, consistency and convergence, WSEAS Transactions on Mathematics 3 (2004) 616–619. [24] A. Saadatmandi, M. Dehghan, Computation of two time-dependent coefficients in a parabolic partial differential equation subject to additional specifications, International Journal of Computer Mathematics 87 (2010) 997–1008. [25] N.P. Salihov, Solution of the Goursat problem by multipoint difference methods. I, (Russian) Vestnik Moskovskogo Universiteta. Serija I. Matematika, Mehanika 4 (1961) 25–34. [26] N.P. Salihov, Solution of the Goursat problem by multipoint difference methods. II, (Russian) Vestnik Moskovskogo Universiteta. Serija I. Matematika, Mehanika 6 (1961) 3–16. [27] M. Shamsi, M. Dehghan, Determination of a control function in three-dimensional parabolic equations by Legendre pseudospectral method, Numerical Methods for Partial Differential E qs. 28 (2012) 74–93. [28] H.J. Stetter, W. Torning, General Multistep Finite Difference Methods for the Solution of uxy ¼ f ðx; y; u; ux ; uy Þ, Rendiconti del Circolo Matematico di Palermo 12 (1963) 281–298. [29] A.M. Wazwaz, On the numerical solution of the Goursat problem, Applied Mathematics and Computation 59 (1993) 89–95. [30] A.M. Wazwaz, The decomposition method for approximate solution of the Goursat problem, Applied Mathematics and Computation 69 (1995) 299– 311. [31] A.M. Wazwaz, The variational iteration method for a reliable treatment of the linear and the nonlinear Goursat problem, Applied Mathematics and Computation 193 (2007) 455–462.
Contents lists available at SciVerse ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
A numerical method for solving a Goursat–Cauchy boundary value problem Mihaela-Cristina Drignei Division of Physical and Computational Sciences, University of Pittsburgh at Bradford, Bradford, PA 16701, USA
a r t i c l e
i n f o
a b s t r a c t
Keywords: Second order hyperbolic partial differential equation Inverse Sturm–Liouville problem Goursat boundary value problem Cauchy boundary value problem
In this paper we propose a method to construct numerically the solution-pair of a second order hyperbolic partial differential equation when two types of boundary conditions on a triangular domain are imposed. Specifically, the PDE is W xx ðx; tÞ W tt ðx; tÞ qðxÞWðx; tÞ ¼ dqðxÞLðx; tÞ, where the unknown functions are Wðx; tÞ and dqðxÞ. One set of boundary conditions refer to one set of the triangle’s boundaries (we call them Goursat type), and the other set of boundary conditions both refer to the third boundary of the triangle (we call them Cauchy type). To solve numerically this boundary value problem we make a change of dependent variables and place a mesh on the triangular domain with lines that are the characteristic curves of the hyperbolic PDE. Then we shall integrate along the characteristic curves of the hyperbolic PDE in such a way that numerical values of Wðx; tÞ and dqðxÞ will be produced at mesh nodes located on vertical lines. The calculations will proceed from the vertical boundary of the domain going backwards to the origin O of the Oxt coordinate system, from one vertical line to another. The trapezoidal formula will be used for numerical integration. Numerical examples will be provided. Ó 2013 Elsevier Inc. All rights reserved.
1. Introduction Our aim in this paper is to construct numerically the solution-pair fdqðxÞ; Wðx; tÞg of the boundary value problem in the triangle D0 :¼ fðx; tÞ 2 R2 : 0 6 t 6 x 6 ag (see Fig. 1):
W xx ðx; tÞ W tt ðx; tÞ qðxÞWðx; tÞ ¼ dqðxÞLðx; tÞ; ðx; tÞ 2 D0 ; Wðx; xÞ ¼
1 2
Z
ð1:1Þ
x
dqðsÞ ds;
0 6 x 6 a;
ð1:2Þ
0
Wðx; 0Þ ¼ 0;
0 6 x 6 a;
Wða; tÞ ¼ f ðtÞ; W x ða; tÞ ¼ gðtÞ;
ð1:3Þ
0 6 t 6 a;
ð1:4Þ
0 6 t 6 a;
ð1:5Þ 2
1
2
where qðxÞ; Lðx; tÞ; f ðtÞ; gðtÞ are known functions in the spaces L ð0; aÞ; CðD0 Þ; H ð0; aÞ; L ð0; aÞ, respectively, and f ð0Þ ¼ 0 ¼ gð0Þ. They are assumed real-valued. Here, E-mail address: [email protected] 0096-3003/$ - see front matter Ó 2013 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.amc.2013.05.041
124
M.-C. Drignei / Applied Mathematics and Computation 220 (2013) 123–141
Fig. 1. The domain of the problem (1.1)–(1.5).
CðD0 Þ is the space of continuous functions on D0 ; L2 ð0; aÞ :¼
u:
Z
a
u2 ðsÞ ds < 1 ;
0
Z H1 ð0; aÞ :¼ u 2 C½0; a : uðxÞ ¼ a þ
x
v ðsÞ ds;
for some a 2 R and some
v 2 L2 ð0; aÞ real-valued
:
0
Note that the requirement f ð0Þ ¼ 0 ¼ gð0Þ is needed in order to make (1.4) and (1.5) compatible with (1.3). The existence and uniqueness of the solution-pair fdqðxÞ; Wðx; tÞg to the problem (1.1)–(1.5) were proved in [[20] Theorem 4.17, pp. 152–154]. Our search for the numerical solution to the above boundary value problem was motivated by the fact that such a problem is to be solved at each stage of the Newton iteration process in order to obtain the numerical solution qðxÞ to the nonlinear equation:
FðqÞ ¼ f~f ðtÞ; g~ðtÞg;
ð1:6Þ
for a given pair of functions f~f ðtÞ; g~ðtÞg in H1 ð0; aÞ L2 ð0; aÞ. Here the map F is defined by
FðqÞ :¼ fKða; t; qÞ; Kx ða; t; qÞg;
ð1:7Þ
where Kðx; t; qÞ is the week solution to the Goursat problem in the triangle D0 :
8 ðx; tÞ 2 D0 ; > tt ðx; tÞ qðxÞKðx; tÞ ¼ 0; < Kxx ðx; tÞ K Rx Kðx; xÞ ¼ 12 0 qðsÞ ds; 0 6 x 6 a; > : Kðx; 0Þ ¼ 0; 0 6 x 6 a:
ð1:8Þ
In the theory of inverse spectral problems, solving (1.6) is equivalent to finding the potential function qðxÞ of the canonical ~ðtÞg obtained from two known spectral data of this operSturm–Liouville operator from knowledge of the Cauchy data f~f ðtÞ; g ator (see Section Appendix A). Then the Newton iterates are:
q0 ¼ initial guess; qnþ1 ¼ qn þ dqn ;
for n ¼ 0; 1; 2; . . . ;
ð1:9Þ
where fdqn ðxÞ; W n ðx; tÞg is the solution-pair to (1.1)–(1.5) with
8 qðxÞ :¼ qn ðxÞ; > > > < Lðx; tÞ :¼ Kðx; t; q Þ; n ~f ðtÞ Kða; t; q Þ; > f ðtÞ :¼ > n > : gðtÞ :¼ g~ðtÞ Kx ða; t; qn Þ:
ð1:10Þ
Complete details about the application of the Newton method to (1.6) will be provided in a different paper. We associated the name Goursat–Cauchy to the boundary value problem (1.1)–(1.5) because of its resemblance to the problem
8 ðx; tÞ 2 D0 ; > tt ðx; tÞ qðxÞWðx; tÞ ¼ 0; < W xx ðx; tÞ W R 1 x Wðx; xÞ ¼ 2 0 qðsÞ ds; 0 6 x 6 a; > : Wðx; 0Þ ¼ 0; 0 6 x 6 a; known in literature as the Goursat problem in the triangle D0 (see the statement right before Theorem 4.15 in [20]), and to the problem
M.-C. Drignei / Applied Mathematics and Computation 220 (2013) 123–141
125
8 > < W xx ðx; tÞ W tt ðx; tÞ qðxÞWðx; tÞ ¼ 0; ðx; tÞ 2 D0 ; Wða; tÞ ¼ f ðtÞ; 0 6 t 6 a; > : W x ða; tÞ ¼ gðtÞ; 0 6 t 6 a; known in literature as the Cauchy problem in the triangle D0 (see the statement right before Theorem 4.16 in [20]). For the past several decades the following boundary value problem, known as the Goursat problem on a rectangular domain has been abundantly studied (e.g. [1–4,9,12,16,15,17,22,23,25,26,28–31]):
8 > < uxy ¼ f ðx; y; u; ux ; uy Þ; 0 6 x 6 a and 0 6 y 6 b; uðx; 0Þ ¼ /ðxÞ; 0 6 x 6 a; > : uð0; yÞ ¼ wðyÞ; 0 6 y 6 b;
ð1:11Þ
where /ð0Þ ¼ wð0Þ, and various kinds of regularity on the known functions f ; /; w were assumed. In [14], the authors study a Goursat problem on a triangular domain. Specifically, they investigate the numerical resolvability of the problem
8 0 6 x 6 1 and x 6 y 6 1; > < uxy ¼ f ðx; y; u; ux ; uy Þ; uð0; yÞ ¼ sðyÞ; 0 6 y 6 1; > : H1 ux ðx; xÞ þ H2 uy ðx; xÞ þ uðx; xÞ ¼ rðxÞ; 0 6 x 6 1;
ð1:12Þ
where the functions f ; s; r are known and the real constants H1 ; H2 are also known, and point out that their method ’’can be used to approximately solve some inverse problems, which are further applied in areas such as geoscience and biomedical engineering’’. In the appendix of [14] the authors describe the connection between their problem (1.12) and an inverse spectral problem: recovery of the potential function of the canonical Sturm–Liouville operator from knowledge of the potential function only on half of its interval (i.e. ’partially known potential’). One difference between (1.12) and our current problem consists in the fact that they are applicable to different inverse spectral problems: ours is used for the recovery of the potential function from two spectral data via the Cauchy data, and (1.12) is applicable to the ’partially known potential’ case. Another difference is that, even with a change of variables, our problem (1.1)–(1.5) cannot be put into the form of (1.12), nor of (1.11). Indeed, by the change of the independent variables
ðx; tÞ $ ðn; gÞ; given by
(
n ¼ xþt ; 2 g ¼ xt ; 2
our domain D0 is transformed into a new domain D0 :¼ f0 6 g 6 n 6 a; and n þ g 6 ag, which is also a triangle (see Fig. 2), the PDE we study, namely (1.1) is transformed into a new PDE:
wng ðn; gÞ qðn þ gÞwðn; gÞ ¼ dqðn þ gÞlðn; gÞ; ðn; gÞ 2 D0 ;
ð1:13Þ
where wðn; gÞ :¼ Wðx; tÞ and lðn; gÞ :¼ Lðx; tÞ, and the boundary conditions (1.2)–(1.5) are transformed respectively into the following boundary conditions:
wðn; 0Þ ¼
1 2
Z
n
dqðsÞ ds;
0 6 n 6 a;
0
Fig. 2. The domain of the problem (1.13)–(1.17).
ð1:14Þ
126
M.-C. Drignei / Applied Mathematics and Computation 220 (2013) 123–141
wðn; nÞ ¼ 0;
06n6
a ; 2
wðn; a nÞ ¼ f ð2n aÞ;
ð1:15Þ a 6 n 6 a; 2
wn ðn; a nÞ þ wg ðn; a nÞ ¼ 2gð2n aÞ;
ð1:16Þ a 6 n 6 a; 2
ð1:17Þ
Thus, our original problem (1.1)–(1.5) is equivalent to the problem (1.13)–(1.17) for the unknown pair of functions fdqðn þ gÞ; wðn; gÞg. Note that (1.13) is not a particular case of the PDE in (1.11), nor in (1.12). Indeed, the PDE (1.13) is to be solved for two unknowns wðn; gÞ and dqðn þ gÞ, while those in (1.11) and in (1.12) are to be solved for only one, namely uðx; yÞ. The boundary conditions (1.14)–(1.17) refer to all boundaries of the domain in our case, whereas in all of the above mentioned papers they refer only to two boundaries of the rectangular or triangular domains. Also, the methods of numerical resolvability of the proposed problems are different in our paper and in Ref. [14]. In the present paper we make a change of dependent variables
Wðx; tÞ $ fZ 1 ðx; tÞ; Z 2 ðx; tÞg and place a mesh on the triangular domain D0 that consists of lines that are the characteristic curves of the hyperbolic PDE (1.1). We integrate along these characteristic curves to produce numerical values of Z 1 ðx; tÞ; Z 2 ðx; tÞ (from which Wðx; tÞ will Rx be obtained), dqðxÞ and rðxÞ :¼ 0 dqðsÞ ds at the grid points. The integration is performed in such a way that we progress from one vertical line to another, starting from the vertical boundary of the domain D0 and going backwards to the origin O of the Oxt coordinate system. Specifically, our each working cell is a triangle with two vertexes on one vertical line, and the third vertex on the neighboring vertical to the left (details and graphical illustrations are provided in Section 2). We utilize the previously calculated values of Z 1 ; Z 2 ; dq; r at the two vertical nodes of the triangular cell to produce their values at the third vertex of the cell. In Ref. [14], the mesh consists of lines parallel to the coordinate axes (x-axis and y-axis), the working cell is a square, the computations proceed from one horizontal line to another, starting from the origin O and going upwards, and from left to right on each horizontal line. Gou and Sun in [14] first find approximate values of u; ux ; uy at all nodes on the vertical boundary x ¼ 0, then at the mesh node ðx1 ; y1 Þ on the slanted boundary y ¼ x. Then for each j P 2 the computations move from the horizontal line y ¼ yj1 to the horizontal line y ¼ yj : for each square cell spanning between the lines y ¼ yj1 and y ¼ yj , approximate values of u; ux ; uy at the upper-right corner are produced from the previously calculated values at the other three corners. The values of u; ux ; uy at the grid point ðxj ; yj Þ on the slanted boundary are obtained either from the values at the grid points ðxi ; yj Þ, with i ¼ 0; 1; . . . ; j 1, or else from the values at the grid points ðxj1 ; yj Þ and ðxj1 ; yj1 Þ, depending on how the impedance boundary parameters H1 ; H2 are chosen. Theoretical approaches to the existence of global solutions to some non-linear Cauchy–Goursat problems for the wave equation have been developed for example in [18,19]. The equations considered are respectively utt uxx juja ¼ f ðx; tÞ, and utt uxx juja1 u ¼ f ðx; tÞ, in fðx; tÞ : 0 < x < t; 0 < t < Tg, with the conditions ux ð0; tÞ ¼ 0 and uðt; tÞ ¼ 0 for 0 6 t 6 T. The function f ðx; tÞ and the parameter a > 0 are known. 2. The method First, we note that the characteristic curves of the hyperbolic PDE (1.1) are the curves
x þ t ¼ constant; x t ¼ constant:
They are obtained by making a change of independent variables
ðx; tÞ $ ðn; gÞ that puts (1.1) into an equivalent, but easier to integrate form - at least for theoretical purposes. Setting wðn; gÞ :¼ Wðx; tÞ, straightforward calculations based on the change of variables n ¼ nðx; tÞ; g ¼ gðx; tÞ lead to the equivalent PDE:
~ðn; gÞwðn; gÞ ¼ dq ~ðn; gÞlðn; gÞ; ðn2x n2t Þwnn þ 2ðnx gx nt gt Þwng þ ðg2x g2t Þwgg þ ðnxx ntt Þwn þ ðgxx gtt Þwg q
ð2:1Þ
~ðn; gÞ :¼ qðxÞ; dq ~ðn; gÞ :¼ dqðxÞ; lðn; gÞ :¼ Lðx; tÞ. Choosing where q
(
n ¼ xþt ; 2
g ¼ xt ; 2 Eq. (2.1) takes the simpler form (1.13). Hence, we have justified the equations of the characteristic curves. However, we shall not pursue Eq. (1.13), but rather (1.1). We make the change of dependent variables
Wðx; tÞ $ fZ 1 ðx; tÞ; Z 2 ðx; tÞg given by:
ð2:2Þ
M.-C. Drignei / Applied Mathematics and Computation 220 (2013) 123–141
Z 1 ðx; tÞ :¼ W x ðx; tÞ þ aðxÞWðx; tÞ þ W t ðx; tÞ;
127
ð2:3Þ
Z 2 ðx; tÞ :¼ W x ðx; tÞ þ aðxÞWðx; tÞ W t ðx; tÞ; where aðxÞ is to be chosen conveniently later, and set
rðxÞ :¼
Z
x
dqðsÞ ds:
ð2:4Þ
0
From (2.3) we can write based on (1.1) that:
(
ðZ 1;x Z 1;t Þðx; tÞ ¼ ða0 ðxÞ a2 ðxÞ þ qðxÞÞWðx; tÞ þ aðxÞZ 2 ðx; tÞ þ dqðxÞLðx; tÞ; ðZ 2;x þ Z 2;t Þðx; tÞ ¼ ða0 ðxÞ a2 ðxÞ þ qðxÞÞWðx; tÞ þ aðxÞZ 1 ðx; tÞ þ dqðxÞLðx; tÞ:
ð2:5Þ
We shall choose aðxÞ such that the Riccati equation is satisfied:
a0 ðxÞ a2 ðxÞ þ qðxÞ ¼ 0; for 0 6 x 6 a:
ð2:6Þ
Inserting (2.6) into (2.5) and taking t ¼ x þ c in the first equation of (2.5), and t ¼ x c in the second equation of (2.5), we obtain:
d Z 1 ðx; x þ cÞ ¼ aðxÞZ 2 ðx; x þ cÞ þ dqðxÞLðx; x þ cÞ; dx
for 0 6 x 6 a;
ð2:7Þ
and respectively,
d Z 2 ðx; x cÞ ¼ aðxÞZ 1 ðx; x cÞ þ dqðxÞLðx; x cÞ; dx
for 0 6 x 6 a:
ð2:8Þ
Utilizing (1.4) and (1.5) in (2.3) we obtain two new boundary conditions for Z 1 and Z 2 on the vertical boundary of the triangular domain D0 :
Z 1 ða; tÞ ¼ gðtÞ þ aðaÞf ðtÞ þ f 0 ðtÞ;
for 0 6 t 6 a;
ð2:9Þ
for 0 6 t 6 a:
ð2:10Þ
and respectively
Z 2 ða; tÞ ¼ gðtÞ þ aðaÞf ðtÞ f 0 ðtÞ; Utilizing (1.3) in (2.3) we write:
Z 1 ðx; 0Þ þ Z 2 ðx; 0Þ ¼ 0;
for 0 6 x 6 a:
ð2:11Þ
Finally, using (1.2) and (2.4) in the first equation of (2.3) we obtain:
Z 1 ðx; xÞ ¼
1 ðdqðxÞ þ aðxÞrðxÞÞ; 2
for 0 6 x 6 a:
ð2:12Þ
Note that (2.9)–(2.11) do not lead to a value for aðaÞ. They rather show a compatibility relationship, for if we make t ¼ 0 in (2.9) and (2.10), and add them together we obtain
Z 1 ða; 0Þ þ Z 2 ða; 0Þ ¼ 2ðgð0Þ þ aðaÞf ð0ÞÞ; which on grounds of (2.11) with x ¼ a leads to
gð0Þ þ aðaÞf ð0Þ ¼ 0:
ð2:13Þ
The hypothesis f ð0Þ ¼ 0 ¼ gð0Þ (see (1.1)–(1.5)) makes (2.13) true for any value of aðaÞ. We do need a value for aðxÞ either at x ¼ 0, or else at x ¼ a in order to solve uniquely (2.6). We shall prescribe a final condition, aðaÞ. For a discussion on this matter see Subsection 6.1.2 in [10]. Also note that Eqs. (2.7) and (2.8) clearly show the necessity of a mesh that consists of characteristic curves for the hyperbolic PDE (1.1). Specifically, the mesh lines are:
x t ¼ c; where c ¼ x1 ; x3 ; x5 ; . . . ; xN ; x þ t ¼ c; where c ¼ x1 ; x3 ; x5 ; . . . ; xN ; xNþ2 ; . . . ; x2N1 :
ð2:14Þ
Here 0 ¼ x1 < x2 < x3 < < xN ¼ a < xNþ1 < < x2N1 ¼ 2a are equally spaced points, and N is chosen to be an odd number. See Fig. 3. For each fixed c > 0, Eq. (2.7) integrates along the characteristic curve t ¼ x þ c that is parallel to the secondary diagonal (i.e. the calculation of Z 1 progresses parallel to the secondary diagonal of the coordinate system Oxt), and Eq. (2.8) integrates along the characteristic curve t ¼ x c that is parallel to the main diagonal (i.e. the calculation of Z 2 progresses parallel to the main diagonal). See Fig. 4. In the following subsections we shall use (2.4) and (2.7)–(2.12) to find numerical values of Z 1 ; Z 2 ; r; dq at the mesh nodes. But first we shall solve (2.6).
128
M.-C. Drignei / Applied Mathematics and Computation 220 (2013) 123–141
Fig. 3. The mesh on the domain D0 .
2.1. The Riccati equation Eq. (2.6) along with a prescribed value of aðaÞ is a final value problem for the unknown function aðxÞ. We can call a MATLAB ode solver (e.g. ODE45) to solve it. Another way is as follows. Let
uðxÞ :¼
Z
x
qðsÞ ds:
ð2:15Þ
0
Then, using the notations uj :¼ uðxj Þ, for j ¼ 1; 2; . . . ; N we have:
(
u1 ¼ 0; because x1 ¼ 0; ðqj þ qjþ1 Þ; for j ¼ 1; 2; . . . N 1ðtrapezoid formulaÞ: ujþ1 uj ¼ dx 2
ð2:16Þ
The elementary trapezoid formula can be found in [[13], pp 153], or in any textbook on Calculus. Next, using (2.15) into (2.6) we write:
ða þ uÞ0 ðxÞ ¼ a2 ðxÞ;
for 0 6 x 6 a:
Letting bðxÞ :¼ ða þ uÞðxÞ, the above ODE becomes: 0
2
b ðxÞ ¼ ðb uÞ ðxÞ;
for 0 6 x 6 a:
ð2:17Þ
The final value bðaÞ ¼ aðaÞ þ uðaÞ ¼ aðaÞ þ uN is known since aðaÞ is prescribed and uN is found from (2.16). Hence, numerical values of bðxÞ at x ¼ x1 ; x2 ; . . . ; xN can be obtained via (2.17):
Fig. 4. How the calculations of Z 1 and Z 2 at the mesh nodes are to be performed.
M.-C. Drignei / Applied Mathematics and Computation 220 (2013) 123–141
8 b ¼ aðaÞ þ uN ; > < N 0 2 bN1 ¼ bN ðdxÞbN ¼ bN ðdxÞðbN uN Þ ðbackward difference schemeÞ; > : 0 2 bj2 ¼ bj ð2dxÞbj1 ¼ bj ð2dxÞðbj1 uj1 Þ ; for j ¼ N; N 1; . . . ; 3 ðcentered difference schemeÞ:
129
ð2:18Þ
Finally, from (2.16) and (2.18), and the definition of bðxÞ the numerical values of aðxÞ at x ¼ x1 ; x2 ; . . . ; xN will be obtained:
aN ¼ aðaÞ ¼ prescribed; aj ¼ bj uj ; for j ¼ N 1; N 2; . . . ; 2; 1:
ð2:19Þ
2.2. The data on the vertical boundary x ¼ xN ¼ a Numerical values for Z 1 and Z 2 at the grid points on the vertical boundary are obtained via (2.9) and (2.10), because aðaÞ is prescribed by us and f ðtÞ; gðtÞ are given functions. For the calculations of f 0 ðt j Þ; j ¼ 1; 2; . . . ; N t we shall use a forward, centered or backward difference scheme derivable from the Taylor’s series expansion:
8 0 f2 f1 f ¼ dt ; > > < 1 f f j ¼ 2; 3; . . . ; Nt 1; fj0 ¼ jþ12dtj1 ; > > : 0 fNt fNt 1 fNt ¼ dt ; where we used the abbreviations fj :¼ f ðt j Þ and fj0 ¼ f 0 ðtj Þ, and 0 ¼ t1 < t2 < . . . < t Nt ¼ a is the partition in t-direction with equally spaced points. Note that dt ¼ 2dx, and so N t ¼ Nþ1 (hence the necessity of N being an odd number). The forward 2 2 and backward difference scheme are of order OðdtÞ, whereas the centered difference scheme is of order Oðdt Þ. The numerical values of r and dq at the grid points on this boundary x ¼ xN ¼ a (i.e. r N :¼ rðxN Þ; dqN :¼ dqðxN Þ) will be calculated in Subsection 2.3. 2.3. The data on the vertical line x ¼ xN1 We integrate (2.7) with c ¼ x2N3 and (2.8) with c ¼ x1 ¼ 0 over ½xN1 ; xN and use the simple trapezoid rule (see Fig. 5). We obtain:
Z 1 ðGÞ Z 1 ðHÞ ¼
dx ðaN1 Z 2 ðHÞ þ aN Z 2 ðGÞ þ dqN1 LðHÞ þ dqN LðGÞÞ; 2
ð2:20Þ
Z 2 ðFÞ Z 2 ðHÞ ¼
dx ðaN1 Z 1 ðHÞ þ aN Z 1 ðFÞ þ dqN1 LðHÞ þ dqN LðFÞÞ: 2
ð2:21Þ
From (2.12) with x ¼ xN1 and respectively x ¼ xN we obtain:
Z 1 ðHÞ ¼
1 ðdqN1 þ aN1 r N1 Þ; 2
ð2:22Þ
Z 1 ðFÞ ¼
1 ðdqN þ aN rN Þ: 2
ð2:23Þ
Fig. 5. Z 1 ; Z 2 ; dq; r for the highest triangular cell between the lines x ¼ xN and x ¼ xN1 .
130
M.-C. Drignei / Applied Mathematics and Computation 220 (2013) 123–141
From (2.4) we have
rðxN Þ rðxN1 Þ ¼
Z
xN
dqðsÞ ds; xN1
which after the application of the simple trapezoid rule we obtain:
rN r N1 ¼
dx ðdqN1 þ dqN Þ: 2
ð2:24Þ
Taking x ¼ xN ¼ a in (2.4) and using (1.2) with x ¼ a and (1.4) with t ¼ a, we obtain:
rN ¼ 2f N :
ð2:25Þ
We used the abbreviations dqi :¼ dqðxi Þ; r i :¼ rðxi Þ; ai :¼ aðxi Þ. We note that (2.20)–(2.25) is a linear system with six equations and six unknowns: Z 1 ðHÞ; Z 2 ðHÞ; dqN1 ; rN1 ; dqN ; r N , because Z 1 ðGÞ; Z 2 ðGÞ; Z 1 ðFÞ; Z 2 ðFÞ are known from the previous step (see Subsection 2.2). We can solve this system easily. Next, we look at an arbitrary triangular cell between the vertical lines x ¼ xN and x ¼ xN1 . See Fig. 6. For each 2 6 i 6 N1 , 2 we integrate (2.7) with c ¼ x2N2i1 and (2.8) with c ¼ x2i1 over ½xN1 ; xN and use the simple trapezoid rule. We obtain respectively:
Z 1 ðBÞ Z 1 ðCÞ ¼
dx ðaN1 Z 2 ðCÞ þ aN Z 2 ðBÞ þ dqN1 LðCÞ þ dqN LðBÞÞ; 2
ð2:26Þ
Z 2 ðAÞ Z 2 ðCÞ ¼
dx ðaN1 Z 1 ðCÞ þ aN Z 1 ðAÞ þ dqN1 LðCÞ þ dqN LðAÞÞ: 2
ð2:27Þ
We solve the linear system (2.26) - (2.27) for the unknowns Z 1 ðCÞ; Z 2 ðCÞ, because Z 1 ðAÞ; Z 2 ðAÞ; Z 1 ðBÞ; Z 2 ðBÞ were obtained in Subsection 2.2 and dqN1 ; dqN were obtained from (2.20)–(2.25). 2.4. The data on the vertical line x ¼ x2j1 We shall obtain numerical values of Z 1 ; Z 2 ; r; dq at the grid points located on the vertical line x ¼ x2j1 from the values of Z 1 ; Z 2 ; r; dq at the grid points on the vertical line x ¼ x2j we calculated previously. We start with the top triangular cell spanned between the vertical lines x ¼ x2j and x ¼ x2j1 . See Fig. 7. We integrate (2.7) with c ¼ x4j3 , and (2.8) with c ¼ x1 ¼ 0 over ½x2j1 ; x2j and use the simple trapezoid rule. We obtain respectively:
Z 1 ðPÞ Z 1 ðMÞ ¼
dx ða2j1 Z 2 ðMÞ þ a2j Z 2 ðPÞ þ dq2j1 LðMÞ þ dq2j LðPÞÞ; 2
ð2:28Þ
Z 2 ðKÞ Z 2 ðMÞ ¼
dx ða2j1 Z 1 ðMÞ þ a2j Z 1 ðKÞ þ dq2j1 LðMÞ þ dq2j LðKÞÞ: 2
ð2:29Þ
From (2.12) with x ¼ x2j1 we write:
Z 1 ðMÞ ¼
1 ðdq2j1 þ a2j1 r2j1 Þ: 2
From (2.4) we have
Fig. 6. Z 1 ; Z 2 for an arbitrary triangular cell between the lines x ¼ xN and x ¼ xN1 .
ð2:30Þ
M.-C. Drignei / Applied Mathematics and Computation 220 (2013) 123–141
131
Fig. 7. Z 1 ; Z 2 ; dq; r for the highest triangular cell between the lines x ¼ x2j and x ¼ x2j1 .
rðx2j Þ rðx2j1 Þ ¼
Z
x2j
dqðsÞ ds;
x2j1
which after the application of the simple trapezoid rule becomes:
r2j r2j1 ¼
dx ðdq2j1 þ dq2j Þ: 2
ð2:31Þ
We solve the linear system (2.28)–(2.31) for the unknowns Z 1 ðMÞ; Z 2 ðMÞ; r 2j1 ; dq2j1 . Next, we look at an arbitrary cell that spans between the vertical lines x ¼ x2j and x ¼ x2j1 . See Fig. 8. For each 2 6 i 6 j 1, we integrate (2.7) with c ¼ x4j2i1 , and (2.8) with c ¼ x2i1 over ½x2j1 ; x2j and use the simple trapezoid rule. We obtain respectively:
Z 1 ðSÞ Z 1 ðTÞ ¼
dx ða2j1 Z 2 ðTÞ þ a2j Z 2 ðSÞ þ dq2j1 LðTÞ þ dq2j LðSÞÞ; 2
ð2:32Þ
Z 2 ðRÞ Z 2 ðTÞ ¼
dx ða2j1 Z 1 ðTÞ þ a2j Z 1 ðRÞ þ dq2j1 LðTÞ þ dq2j LðRÞÞ: 2
ð2:33Þ
We solve the linear system (2.32) - (2.33) for the unknowns Z 1 ðTÞ; Z 2 ðTÞ. Note that dq2j1 was obtained from (2.28)–(2.31), and dq2j from the calculations done for the previous vertical line x ¼ x2j . Finally, we calculate Z 1 and Z 2 at the lowest grid point on the vertical line x ¼ x2j1 . See Fig. 9. From (2.11) with x ¼ x2j1 we write:
Z 1 ðUÞ þ Z 2 ðUÞ ¼ 0:
ð2:34Þ
Taking c ¼ x2j1 in (2.8) and integrating over ½x2j1 ; x2j we obtain:
Fig. 8. Z 1 ; Z 2 for an arbitrary triangular cell between the lines x ¼ x2j and x ¼ x2j1 .
132
M.-C. Drignei / Applied Mathematics and Computation 220 (2013) 123–141
Fig. 9. Z 1 ; Z 2 at the lowest grid point on the vertical line x ¼ x2j1 .
Z 2 ðVÞ Z 2 ðUÞ ¼
dx ða2j1 Z 1 ðUÞ þ a2j Z 1 ðVÞ þ dq2j1 LðUÞ þ dq2j LðVÞÞ: 2
ð2:35Þ
We solve the linear system (2.34) - (2.35) for Z2(U) for Z 1 ðUÞ; Z 2 ðUÞ. 2.5. The data on the vertical line x ¼ x2j2 We shall perform our calculations from the highest triangular cell to the lowest triangular cell spanned between the vertical lines x ¼ x2j1 and x ¼ x2j2 . The highest cell is illustrated in Fig. 10. Integrating (2.7) with c ¼ x4j5 , and (2.8) with c ¼ x1 ¼ 0 over ½x2j2 ; x2j1 and using the simple trapezoid rule, we obtain respectively:
Z 1 ðZÞ Z 1 ðXÞ ¼
dx ða2j2 Z 2 ðXÞ þ a2j1 Z 2 ðZÞ þ dq2j2 LðXÞ þ dq2j1 LðZÞÞ; 2
ð2:36Þ
Z 2 ðYÞ Z 2 ðXÞ ¼
dx ða2j2 Z 1 ðXÞ þ a2j1 Z 1 ðYÞ þ dq2j2 LðXÞ þ dq2j1 LðYÞÞ: 2
ð2:37Þ
From (2.12) with x ¼ x2j2 we obtain:
Z 1 ðXÞ ¼
1 ðdq2j2 þ a2j2 r 2j2 Þ: 2
From (2.4) we have:
rðx2j1 Þ rðx2j2 Þ ¼
Z
x2j1
dqðsÞ ds;
x2j2
which after applying the simple trapezoid rule gives:
Fig. 10. Z 1 ; Z 2 ; dq; r for the highest triangular cell between the lines x ¼ x2j1 and x ¼ x2j2 .
ð2:38Þ
M.-C. Drignei / Applied Mathematics and Computation 220 (2013) 123–141
r2j1 r 2j2 ¼
dx ðdq2j2 þ dq2j1 Þ: 2
133
ð2:39Þ
We solve the linear system (2.36)–(2.39) for the unknowns Z 1 ðXÞ; Z 2 ðXÞ; r 2j2 and dq2j2 . Note that Z 1 ðZÞ; Z 2 ðZÞ; Z 1 ðYÞ; Z 2 ðYÞ; r 2j1 ; dq2j1 were obtained in Subsection 2.4. Next we look at an arbitrary triangular cell spanned between the vertical lines x ¼ x2j1 and x ¼ x2j2 . See Fig. 11. For each 2 6 i 6 j 1, we integrate (2.7) with c ¼ x4j2i3 , and (2.8) with c ¼ x2i1 over ½x2j2 ; x2j1 . We obtain:
Z 1 ðEÞ Z 1 ðIÞ ¼
dx ða2j2 Z 2 ðIÞ þ a2j1 Z 2 ðEÞ þ dq2j2 LðIÞ þ dq2j1 LðEÞÞ; 2
ð2:40Þ
Z 2 ðDÞ Z 2 ðIÞ ¼
dx ða2j2 Z 1 ðIÞ þ a2j1 Z 1 ðDÞ þ dq2j2 LðIÞ þ dq2j1 LðDÞÞ: 2
ð2:41Þ
We solve the linear system (2.40) - (2.41) for the unknowns Z 1 ðIÞ; Z 2 ðIÞ. Note that dq2j2 was obtained from (2.36)–(2.39)), and Z 1 ðDÞ; Z 2 ðDÞ; Z 1 ðEÞ; Z 2 ðEÞ; dq2j1 were obtained in Subsection 2.4. 2.6. The data on the vertical line x ¼ x1 There is only one grid point on the vertical line x ¼ x1 , namely O – the origin of the coordinate system. See Fig. 12. Taking x ¼ x1 ¼ 0 in (2.11) we have:
Z 1 ðOÞ þ Z 2 ðOÞ ¼ 0:
ð2:42Þ
Integrating (2.8) with c ¼ x1 ¼ 0 over ½x1 ; x2 and using the simple trapezoid formula, we obtain:
Z 2 ðQ Þ Z 2 ðOÞ ¼
dx ða1 Z 1 ðOÞ þ a2 Z 1 ðQ Þ þ dq1 LðOÞ þ dq2 LðQ ÞÞ: 2
ð2:43Þ
From (2.12) with x ¼ x1 ¼ 0 we have:
Z 1 ðOÞ ¼
1 ðdq1 þ a1 r 1 Þ: 2
ð2:44Þ
Due to (2.4) we write:
rðx2 Þ rðx1 Þ ¼
Z
x2
dqðsÞ ds; x1
which after applying the simple trapezoid rule gives:
r2 r1 ¼
dx ðdq1 þ dq2 Þ: 2
We solve the linear system (2.42)–(2.45) for the unknowns Z 1 ðOÞ; Z 2 ðOÞ; r 1 and dq1 . 2.7. The calculation of W; W x ; W t at the mesh nodes From (2.3) we can obtain W t ðx; tÞ and W x ðx; tÞ þ aðxÞWðx; tÞ at all mesh nodes ðx; tÞ:
Fig. 11. Z 1 ; Z 2 for an arbitrary triangular cell between the lines x ¼ x2j1 and x ¼ x2j2 .
ð2:45Þ
134
M.-C. Drignei / Applied Mathematics and Computation 220 (2013) 123–141
Fig. 12. Z 1 ; Z 2 ; dq; r at the lowest grid point O.
W t ðx; tÞ ¼
Z 1 ðx; tÞ Z 2 ðx; tÞ ; 2
W x ðx; tÞ þ aðxÞWðx; tÞ ¼
Z 1 ðx; tÞ þ Z 2 ðx; tÞ : 2
ð2:46Þ
ð2:47Þ
Since ri :¼ rðxi Þ, for i ¼ 1; 2; . . . ; N were calculated in Subsections 2.3 - 2.6 we will also be able to calculate Wðx; xÞ via (1.2) and (2.4) at all mesh nodes on the slanted boundary. Letting Wðk; lÞ :¼ Wðxk ; xl Þ; W x ðk; lÞ :¼ W x ðxk ; xl Þ; W t ðk; lÞ :¼ W t ðxk ; xl Þ and using (2.46), (2.47), one way to obtain W t ; W; W x at all mesh nodes is as follows. For each k ¼ N; N 1; . . . ; 6; 5, perform calculations on the vertical line x ¼ xk going from top to bottom:
8 Z 1 ðk;kÞZ 2 ðk;kÞ ; > 2 < W t ðk; kÞ ¼ 1 Wðk; kÞ ¼ 2 r k ; > : 2 ðk;kÞ W x ðk; kÞ ¼ Z1 ðk;kÞþZ ak Wðk; kÞ; 2 8 Z 1 ðk;k2ÞZ 2 ðk;k2Þ ; > 2 < W t ðk; k 2Þ ¼ Wðk; k 2Þ ¼ Wðk; kÞ ð2dxÞW t ðk; kÞ; ðbackward difference schemeÞ; > : 2 ðk;k2Þ ak Wðk; k 2Þ; W x ðk; k 2Þ ¼ Z1 ðk;k2ÞþZ 2 and for each i ¼ k; k 2; k 4; . . . ; 7; 5 (if k is odd), or i ¼ k; k 2; k 4; . . . ; 8; 6 (if k is even)
8 Z 1 ðk;i4ÞZ 2 ðk;i4Þ ; > 2 < W t ðk; i 4Þ ¼ Wðk; i 4Þ ¼ Wðk; iÞ ð4dxÞW t ðk; i 2Þ; ðcentered difference schemeÞ; > : 2 ðk;i4Þ ak Wðk; i 4Þ: W x ðk; i 4Þ ¼ Z1 ðk;i4ÞþZ 2 For k ¼ 4 and k ¼ 3, perform calculations on the vertical line x ¼ xk going from top to bottom:
8 Z 1 ðk;kÞZ 2 ðk;kÞ ; > 2 < W t ðk; kÞ ¼ 1 Wðk; kÞ ¼ 2 r k ; > : 2 ðk;kÞ W x ðk; kÞ ¼ Z1 ðk;kÞþZ ak Wðk; kÞ; 2 8 Z 1 ðk;k2ÞZ 2 ðk;k2Þ ; > 2 < W t ðk; k 2Þ ¼ Wðk; k 2Þ ¼ Wðk; kÞ ð2dxÞW t ðk; kÞ; ðbackward difference schemeÞ; > : 2 ðk;k2Þ ak Wðk; k 2Þ; W x ðk; k 2Þ ¼ Z1 ðk;k2ÞþZ 2 For k ¼ 2, set
8 Z 1 ðk;kÞZ 2 ðk;kÞ ; > 2 < W t ðk; kÞ ¼ 1 Wðk; kÞ ¼ 2 r k ; > : 2 ðk;kÞ W x ðk; kÞ ¼ Z1 ðk;kÞþZ ak Wðk; kÞ; 2 For k ¼ 1, set
M.-C. Drignei / Applied Mathematics and Computation 220 (2013) 123–141
135
8 Z 1 ðk;kÞZ 2 ðk;kÞ > ; < W t ðk; kÞ ¼ 2 Wðk; kÞ ¼ 0; ðdue to ð1:3Þ with x ¼ x1 ¼ 0Þ; > : W x ðk; kÞ ¼ 0; ðdue to ð1:3Þ with x ¼ x1 ¼ 0Þ: 3. The algorithm Let N be an odd number. Let 0 ¼ x1 < x2 < < xN ¼ a < xNþ1 < xNþ2 < < x2N1 ¼ 2a be equally spaced points. Place the mesh (2.14) on the domain D0 . Find Z 1 and Z 2 at the grid points on the vertical boundary x ¼ xN ¼ a. See Subsection 2.2. Find Z 1 ; Z 2 at the grid points on the vertical line x ¼ xN1 and rN ; dqN ; rN1 ; dqN1 . See Subsection 2.3. For each j ¼ N1 ; N3 ; . . . ; 4; 3; 2 2 2 – find Z 1 ; Z 2 ; dq; r at the grid points on the vertical line x ¼ x2j1 (see Subsection 2.4), – find Z 1 ; Z 2 ; dq; r at the grid points on the vertical line x ¼ x2j2 (see Subsection 2.5). Find Z 1 ; Z 2 ; dq1 ; r1 at the grid points on the vertical line x ¼ x1 ¼ 0. (See Subsection 2.6), Find W; W x ; W t at the mesh nodes. (See Subsection 2.7). 4. Remarks
Remark 1. The ODE (2.8) with c ¼ x1 ¼ 0 can be further manipulated due to (2.12), (2.4) and (2.6):
d 1 Z 2 ðx; xÞ ¼ aðxÞZ 1 ðx; xÞ þ dqðxÞLðx; xÞ ¼ aðxÞ ðdqðxÞ þ aðxÞrðxÞÞ þ dqðxÞLðx; xÞ dx 2 1 1 0 2 ¼ ðaðxÞr ðxÞ þ a ðxÞrðxÞÞ þ dqðxÞLðx; xÞ ¼ ðaðxÞr 0 ðxÞ þ a0 ðxÞrðxÞ þ qðxÞrðxÞÞ þ dqðxÞLðx; xÞ; 2 2
for 0 6 x 6 a;
which is equivalent to
d 1 1 Z 2 ðx; xÞ aðxÞrðxÞ ¼ qðxÞrðxÞ þ dqðxÞLðx; xÞ; dx 2 2
for 0 6 x 6 a:
ð4:1Þ
Integrating (4.1) over appropriate intervals and using the simple trapezoid formula we produce alternative equations to (2.21), (2.29), (2.37), and (2.43), respectively:
1 1 dx dx Z 2 ðFÞ Z 2 ðHÞ ðaN r N aN1 r N1 Þ ¼ ðqN1 r N1 þ qN rN Þ þ ðdqN1 LðHÞ þ dqN LðFÞÞ; 2 2 2 2
ð4:2Þ
1 1 dx dx Z 2 ðKÞ Z 2 ðMÞ ða2j r 2j a2j1 r 2j1 Þ ¼ ðq2j1 r 2j1 þ q2j r 2j Þ þ ðdq2j1 LðMÞ þ dq2j LðKÞÞ; 2 2 2 2
ð4:3Þ
1 1 dx dx Z 2 ðYÞ Z 2 ðXÞ ða2j1 r2j1 a2j2 r2j2 Þ ¼ ðq2j2 r 2j2 þ q2j1 r2j1 Þ þ ðdq2j2 LðXÞ þ dq2j1 LðYÞÞ; 2 2 2 2
ð4:4Þ
1 1 dx dx Z 2 ðQ Þ Z 2 ðOÞ ða2 r 2 a1 r1 Þ ¼ ðq1 r1 þ q2 r2 Þ þ ðdq1 LðOÞ þ dq2 LðQ ÞÞ: 2 2 2 2
ð4:5Þ
There is however not a significant difference in the algorithm’s accuracy when using formulas (4.2)–(4.5) in place of (2.21), (2.29), (2.37), (2.43), respectively.
5. Numerical examples We illustrate the effectiveness of our numerical method/algorithm on the following two examples where analytical solutions are available, thus making the comparison with the numerical solutions possible. To obtain these examples, we tackled the boundary value problem (1.1)–(1.5) backwards: we rather started with a given formula for Wðx; tÞ and produced dqðxÞ; qðxÞ; Lðx; tÞ; f ðtÞ; gðtÞ such that (1.1)–(1.5) are met. qðxÞ t ð1þ2x2 2 Þ 2 2 2 2 Example 1. For Wðx; tÞ ¼ ex t, we obtain dqðxÞ ¼ 2ex ð1 þ 2x2 Þ; Lðx; tÞ ¼ ; f ðtÞ ¼ ea t; gðtÞ ¼ 2aea t. Thus, the 1þ2x2 2 hypothesis condition f ð0Þ ¼ 0 ¼ gð0Þ is satisfied. We chose qðxÞ ¼ 1þx3 and a ¼ 1, for the numerical illustration.
136
M.-C. Drignei / Applied Mathematics and Computation 220 (2013) 123–141
Fig. 13. Example 1 – the comparison for dqðxÞ, when N ¼ 15 x-partition points have been chosen.
Fig. 14. Example 1 – illustration of
Fig. 15. Example 1 – illustration of
W exact ðx;tÞW numeric ðx;tÞ , kW exact k2
when N ¼ 15 x-partition points have been chosen.
W exact x ðx;tÞW numeric ðx;tÞ x Þ, k2 kW exact x
when N ¼ 15 x-partition points have been chosen.
½1þqðxÞðsin xþcos xÞt Example 2. For Wðx; tÞ ¼ t ðsin x þ cos xÞ, we obtain dqðxÞ ¼ 2½ð1 xÞ sin x þ ð1 þ xÞ cos x; Lðx; tÞ ¼ 2½ð1xÞ ; f ðtÞ ¼ sin xþð1þxÞ cos x ðsin a þ cos aÞ t; gðtÞ ¼ ðcos a sin aÞ t, so that f ð0Þ ¼ 0 ¼ gð0Þ, as required. For the numerical illustration we took a ¼ 1 and 1 qðxÞ ¼ sin xþcos . x MATLAB figure-outputs are provided for our two examples with N ¼ 15 partition points in x-direction. See Figs. 13–20, where the exact solution dqðxÞ (in blue) and the numerical solution dqðxÞ (in red) are displayed for comparison, as well as the divided differences:
W exact ðx; tÞ W numeric ðx; tÞ ; kW exact k2
M.-C. Drignei / Applied Mathematics and Computation 220 (2013) 123–141
Fig. 16. Example 1 – illustration of
W exact ðx;tÞW numeric ðx;tÞ t t , k2 kW exact t
when N ¼ 15 x-partition points have been chosen.
Fig. 17. Example 2 – the comparison for dqðxÞ, when N ¼ 15 x-partition points have been chosen.
Fig. 18. Example 2 – illustration of
W exact ðx; tÞ W numeric ðx; tÞ x x ; k2 kW exact x W exact ðx; tÞ W numeric ðx; tÞ t t : k kW exact 2 t
W exact ðx;tÞW numeric ðx;tÞ , kW exact k2
when N ¼ 15 x-partition points have been chosen.
137
138
M.-C. Drignei / Applied Mathematics and Computation 220 (2013) 123–141
Fig. 19. Example 2 – illustration of
Fig. 20. Example 2 – illustration of
W exact ðx;tÞW numeric ðx;tÞ x x Þ, kW exact k2 x
W exact ðx;tÞW numeric ðx;tÞ t t , k2 kW exact t
when N ¼ 15 x-partition points have been chosen.
when N ¼ 15 x-partition points have been chosen.
The reason we chose to display the derivatives W x ðx; tÞ and W t ðx; tÞ in addition to the solution fdqðxÞ; Wðx; tÞg to (1.1)–(1.5) is because they are directly involved in the computation of Wðx; tÞ (see Subsection 2.7).
6. Conclusions We solved a boundary value problem on a triangular domain for two unknown functions. The numerical method consisted of placing a mesh on the triangular domain with lines that are the characteristic curves of the hyperbolic PDE in the boundary value problem. Then a change of dependent variables was made and the numerical values for the needed unknown functions were obtained at the mesh nodes by going from one vertical line to another, starting from the vertical boundary. Integration was performed along the characteristic curves of the original PDE, and the simple trapezoid formula was utilized to carry out these calculations. Numerical examples were provided to illustrate the effectiveness of the method/ algorithm. This algorithm can be further embedded in a Newton type algorithm through which the numerical solution to an inverse Sturm–Liouville problem can be obtained, provided two spectral data are known (i.e. an inverse two-spectra problem can be solved numerically). Appendix A In this section we present the connection with the inverse two-spectra problems for the canonical Sturm–Liouville operator. The canonical Sturm–Liouville equation on a finite interval ½0; a is
u00 ðxÞ þ qðxÞuðxÞ ¼ kuðxÞ;
0 6 x 6 a;
ð7:1Þ
where qðxÞ is the coefficient function (also called the potential function of the canonical Sturm–Liouville operator). Two boundary conditions are associated to (7.1), one per end-point. The type of boundary conditions gives the name of the boundary value problem so formed. To fix the ideas, let us consider the following two situations:
139
M.-C. Drignei / Applied Mathematics and Computation 220 (2013) 123–141
uð0Þ ¼ 0 ¼ uðaÞ
ð7:2Þ
uð0Þ ¼ 0 ¼ u0 ðaÞ:
ð7:3Þ
and
Thus, (7.1) & (7.2) is called the Dirichlet eigenvalue problem for the canonical Sturm–Liouville operator, because both boundary conditions are Dirichlet type, and (7.1) & (7.3) is refereed to as the Dirichlet–Neumann eigenvalue problem for the canonical Sturm–Liouville operator, because the boundary condition at the left end point x ¼ 0 is Dirichlet type, and the boundary condition at the right end point x ¼ a is Neumann type. The direct problem is formulated as follows: given qðxÞ in a certain space of functions (e.g. L2 ð0; aÞ; H1 ð0; aÞ), find the eigenpairs fk; uðxÞg of either one of the BVP mentioned above, or any other BVP that consists of (7.1) and two boundary conditions (one per end-point). The inverse problem formulates as follows: given information about the spectral data, find qðxÞ. Note that the spectral data do not mean the eigenpairs, because if so, then the problem would become trivial (i.e. even with the knowledge of only one eigenpair you can get qðxÞ from (7.1) as if you would solve a linear equation). Spectral data mean for instance, one sequence of eigenvalues and one sequence of norming constants. Or they can be two sequences of eigenvalues. We discuss the two-spectra case. This means, we know two sequences of eigenvalues. And since we decided to focus on the boundary conditions (7.2) and (7.3), we shall go ahead with the sequence of Dirichlet eigenvalues fkn gnP1 and the sequence of Dirichlet–Neumann eigenvalues fmn gnP1 . That means, fkn ; un ðxÞgnP1 satisfy (7.1) & (7.2), and fmn ; v n ðxÞgnP1 satisfy (7.1) & (7.3), for the same potential function qðxÞ. The facts that the eigenvalues are real (provided qðxÞ is a real-valued function), countably many, and satisfy specific asymptotic formulas were rigorously established in the literature devoted to the direct Sturm–Liouville problems. See for example [[20] Section 4.3]:
kn ¼
np2 þ ½q þ cn ; a
for n ¼ 1; 2; 3; . . . ;
where fcn gnP1 is an l2 sequence of real numbers (i.e.
mn
2 ðn 12Þp ¼ þ ½q þ dn ; a
ð7:4Þ P1
2
n¼1 ðcn Þ
< 1), and
for n ¼ 1; 2; 3; . . . ;
ð7:5Þ
Ra with fdn gnP1 2 l2 . Here ½q :¼ 1a 0 qðxÞ dx is the mean of the function q 2 L2R ð0; aÞ. Also, one can prove that k is a Dirichlet eigenvalue of the canonical Sturm–Liouville operator with potential function qðxÞ, if and only if k is a zero of the function
k 2 C ! Sða; q; kÞ;
ð7:6Þ
called the characteristic function of the BVP (7.1) & (7.2), and m is a Dirichlet–Neumann eigenvalue of the canonical Sturm– Liouville operator with potential function qðxÞ, if and only if m is a zero of the function
k 2 C ! S0 ða; q; kÞ;
ð7:7Þ 2
called the characteristic function of the BVP (7.1) & (7.3). Here, Sðx; q; kÞ stands for the weak solution (i.e. in H ð0; aÞ) of the initial value problem
8 00 > < u ðxÞ þ qðxÞuðxÞ ¼ kuðxÞ; uð0Þ ¼ 0; > : 0 u ð0Þ ¼ 1:
0 6 x 6 a; ð7:8Þ
pffiffi The notation SðÞ was chosen to resemble the function sine, specifically sinðpffiffikkxÞ which is the weak solution to (7.8) with qðxÞ ¼ 0. We choose to detail for the second case, the equivalence between being an eigenvalue and being a zero of the characteristic function. For the other case, the reasons are similar. So, let us assume that m is a Dirichlet–Neumann eigenvalue. Let v ðxÞ be an associated eigenfunction. That is, fm ; v ðxÞg satisfies (7.1) & (7.3). It follows that ðv Þ0 ð0Þ – 0, as otherwise (7.1) and the first boundary condition of (7.3) would imply that v ðxÞ ¼ 0, for all x 2 ½0; a, because the initial value problem (7.1) & uð0Þ ¼ 0 ¼ u0 ð0Þ has only one solution, and the identically zero function is a solution. We cannot have v ðxÞ ¼ 0, for all ~ ðxÞ :¼ ðvv ÞðxÞ ~ ðxÞ satisfies (7.8) x 2 ½0; a, because v ðxÞ is an eigenfunction. Thus, v is well-defined. It follows immediately that v 0 ð0Þ 0 ðaÞ 0 ~ ðxÞ ¼ Sðx; q; m Þ. Hence, S ða; q; m Þ ¼ v ~ 0 ðaÞ ¼ ððvv ÞÞ0 ð0Þ with k ¼ m , and therefore v ¼ 0. We showed that a Dirichlet–Neumann eigenvalue m is a zero of (7.7). Conversely, let m be a zero of (7.7). Then, due to the definition of Sðx; q; kÞ with k ¼ m , we have that fm ; Sðx; q; m Þg satisfies (7.1) & (7.3). This means nothing but m is a Dirichlet–Neumann eigenvalue of the canonical Sturm–Liouville operator with potential function qðxÞ. Hence, we have
Sða; q; kn Þ ¼ 0;
for n ¼ 1; 2; 3; . . . ;
ð7:9Þ
S0 ða; q; mn Þ ¼ 0;
for n ¼ 1; 2; 3; . . . :
ð7:10Þ
and
140
M.-C. Drignei / Applied Mathematics and Computation 220 (2013) 123–141
It is also known from the theory of direct Sturm–Liouville problems that Sðx; q; kÞ has an integral representation (see for instance [[20] Example 4.19 and Theorem 4.18]):
Sðx; q; kÞ ¼
pffiffiffi pffiffiffi Z x sinð kxÞ sinð ktÞ pffiffiffi þ dt; Kðx; t; qÞ pffiffiffi k k 0
for 0 6 x 6 a:
ð7:11Þ
Here Kðx; t; qÞ goes under the name of Gelfand-Levitan-Marchenko kernel. As seen in (7.11), the solution to one initial value problem, namely (7.8) with qðxÞ ¼ 0, is connected to the solution to another initial value problem (7.8) by means of a second kind Volterra integral operator with kernel Kðx; t; qÞ. The kernel satisfies the Goursat problem (1.8) in the triangle D0 :¼ fðx; tÞ 2 R2 : 0 6 t 6 x 6 ag. From (7.11) we obtain:
pffiffiffi pffiffiffi Z x pffiffiffi sinð kxÞ sinð ktÞ þ dt; S0 ðx; q; kÞ ¼ cosð kxÞ þ Kðx; x; qÞ pffiffiffi Kx ðx; t; qÞ pffiffiffi k k 0
for 0 6 x 6 a:
ð7:12Þ
From (7.9), (7.10), (7.11), (7.12), and the first boundary condition of (1.8) we obtain that Kða; t; qÞ is the solution ~f ðtÞ to
Z pffiffiffiffiffi sinð kn aÞ þ
a
pffiffiffiffiffi ~f ðtÞ sinð kn tÞ dt ¼ 0;
for n ¼ 1; 2; 3; . . . ;
ð7:13Þ
0
and Kx ða; t; qÞ is the solution g~ðtÞ to
pffiffiffiffiffi
pffiffiffiffiffi
mn cosð mn aÞ þ
pffiffiffiffiffi a½q sinð mn aÞ þ 2
Z
a
g~ðtÞ sinð
pffiffiffiffiffi
mn tÞ dt ¼ 0; for n ¼ 1; 2; 3; . . . :
ð7:14Þ
0
Note that in (7.14), ½q - the mean of the function qðxÞ can be estimated from the eigenvalues. Due to (7.4) and (7.5), we can write
kN1
2 2 ðN2 12Þp N1 p ½q mN2 ; a a
for N 1 ; N 2 large enough. Hence, everything in (7.13) and (7.14) depends only on the known quantities fkn gnP1 and fmn gnP1 . After solving the systems of integral equations (7.13) and (7.14), we can proceed to the determination of the function qðxÞ – the unknown in the inverse two-spectra problem, by applying the Newton method to the non-linear equation (1.6), where the map FðqÞ is defined by (1.7). ~ðtÞ with Note that, numerically we can solve (7.13) and (7.14) by writing a Fourier series expansion for each of ~f ðtÞ and g 1 2 appropriate bases of functions. Since in practice only the first several eigenvalues are available, say fkn gNn¼1 and fmn gNn¼1 , we truncate these Fourier series after the first N 1 , and respectively N 2 terms. Thus, we will need to keep the equations that correspond to n ¼ 1; 2; . . . N 1 from (7.13), and the equations that correspond to n ¼ 1; 2; . . . N 2 from (7.14). The two linear systems so obtained will be square N 1 N 1 and respectively N 2 N 2 linear systems for the first N 1 and N 2 Fourier coefficients of ~f ðtÞ and g~ðtÞ, respectively. Once these coefficients are calculated, we have estimates of ~f ðtÞ and g~ðtÞ and we can continue as outlined above. We close this section by saying that a canonical Sturm–Liouville equation models the infinitesimal, vertical vibration of an elastic string of negligible mass, and the boundary conditions associated to the equation are derivable from the constraints the ends of the string are subjected to (see for example [[10] Chapter 2 and Section 5.2]). The reader may find of interest the papers [5–8,11,14,21,24,27], which are either on inverse problems or on numerical partial differential equations. Appendix B. Supplementary data Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/ j.amc.2013.05.041. References [1] R.T. Dames, Stability and convergence for a numerical solution of the Goursat problem, Journal of Mathematical Physics 38 (1959–1960) 42–67. [2] J.T. Day, A Gaussian quadrature method for the numerical solution of the characteristic initial value problem uxy ¼ f ðx; y; uÞ, Mathematics of Computation 17 (1963) 438–441. [3] J.T. Day, A Runge-Kutta method for the numerical solution of the Goursat problem in hyperbolic partial differential equations, The Computer Journal 9 (1966) 81–83. [4] J.T. Day, On the numerical solution of the Goursat problem, Nordisk Tidskrift for Informationsbehandling (BIT) 11 (1971) 450–454. [5] M. Dehghan, On the solution of an initial-boundary value problem that combines Neumann and integral condition for the wave equation, Numerical Methods for Partial Differential Eqs. 21 (2005) 24–40. [6] M. Dehghan, A. Ghesmati, Combination of meshless local weak and strong (MLWS) forms to solve the two dimensional hyperbolic telegraph equation, Engineering Analysis with Boundary Elements 34 (2010) 324–336. [7] M. Dehghan, A. Mohebbi, High order implicit collocation method for the solution of two-dimensional linear hyperbolic equation, Numerical Methods for Partial Differential E qs. 25 (2009) 232–243. [8] M. Dehghan, R. Salehi, A method based on meshless approach for the numerical solution of the two-space dimensional hyperbolic telegraph equation, Mathematical Methods in the Applied Science 35 (2012) 1220–1233.
M.-C. Drignei / Applied Mathematics and Computation 220 (2013) 123–141
141
[9] J. Diaz, On an analogue of the Euler–Cauchy polygon method for the numerical solution of uxy ¼ f ðx; y; u; ux ; uy Þ, Archive for Rational Mechanics and Analysis 1 (1957) 154–180. [10] M.C. Drignei, Inverse Sturm–Liouville problems using multiple spectra, Ph.D Thesis, Iowa State University, 2008. [11] M.C. Drignei, Constructibility of an L2R ð0; aÞ solution to an inverse Sturm–Liouville problem using three Dirichlet spectra, Inverse Problems 26 (2010) 025003 (29pp). [12] D.J. Evans, B.B. Sanugi, Numerical solution of the Goursat problem by a nonlinear trapezoidal formula, Applied Mathematics Letters 1 (3) (1988) 221– 223. [13] W. Gautschi, Numerical Analysis an Introduction, Birkhauser, Boston-Besel-Berlin, 1997. [14] K. Gou, B. Sun, Numerical solution of the Goursat problem on a triangular domain with mixed boundary conditions, Applied Mathematics and Computation 217 (2011) 8765–8777. [15] N. Gradinaru, G. Teodoru, On certain numerical methods for the Goursat problem, (Romanian) Buletinul Institutului Politehnic din Iasi, Serie Noua 20 (24) (1974) 43–46. [16] E. Goursat, A Course in Mathematical Analysis, vol. 3, Dover, New York, 2006. [17] M.K. Jain, K.D. Sharma, Cubature method for the numerical solution of the characteristic initial value problem, Journal of the Australian Mathematical Society 8 (1968) 355–368. [18] O. Jokhadze, Cauchy–Goursat problem for one-dimensional semilinear wave equations, Communications in Partial Differential Equations 34 (4–6) (2009) 367–382. [19] O. Jokhadze, On existence and nonexistence of global solutions of Cauchy–Goursat problem for nonlinear wave equations, Journal of Mathematical Analysis and Applications 2 (2008) 1033–1045. [20] A. Kirsch, An introduction to the mathematical theory of inverse problems, Applied Mathematical Sciences 120, Springer-Verlag, New York, 1996. [21] M. Lakestani, M. Dehghan, The use of Chebyshev cardinal functions for the solution of a partial differential equation with an unknown time-dependent coefficient subject to an extra measurement, Journal of Computational and Applied Mathematics 235 (2010) 669–678. [22] R.H. Moore, A Runge-Kutta procedure for the Goursat problem in hyperbolic partial differential equations, Archive for Rational Mechanics and Analysis 7 (1961) 37–63. [23] M.A.S. Nasir, A.I.M. Ismail, Numerical solution of a linear Goursat problem: stability, consistency and convergence, WSEAS Transactions on Mathematics 3 (2004) 616–619. [24] A. Saadatmandi, M. Dehghan, Computation of two time-dependent coefficients in a parabolic partial differential equation subject to additional specifications, International Journal of Computer Mathematics 87 (2010) 997–1008. [25] N.P. Salihov, Solution of the Goursat problem by multipoint difference methods. I, (Russian) Vestnik Moskovskogo Universiteta. Serija I. Matematika, Mehanika 4 (1961) 25–34. [26] N.P. Salihov, Solution of the Goursat problem by multipoint difference methods. II, (Russian) Vestnik Moskovskogo Universiteta. Serija I. Matematika, Mehanika 6 (1961) 3–16. [27] M. Shamsi, M. Dehghan, Determination of a control function in three-dimensional parabolic equations by Legendre pseudospectral method, Numerical Methods for Partial Differential E qs. 28 (2012) 74–93. [28] H.J. Stetter, W. Torning, General Multistep Finite Difference Methods for the Solution of uxy ¼ f ðx; y; u; ux ; uy Þ, Rendiconti del Circolo Matematico di Palermo 12 (1963) 281–298. [29] A.M. Wazwaz, On the numerical solution of the Goursat problem, Applied Mathematics and Computation 59 (1993) 89–95. [30] A.M. Wazwaz, The decomposition method for approximate solution of the Goursat problem, Applied Mathematics and Computation 69 (1995) 299– 311. [31] A.M. Wazwaz, The variational iteration method for a reliable treatment of the linear and the nonlinear Goursat problem, Applied Mathematics and Computation 193 (2007) 455–462.