119 

 

Ibn Al-Haitham Jour. for Pure & Appl. Sci. 34 (2) 2021 

 

 

  

 

 

The Galer kin -I mplicit Methods for Solving Nonlinear Hyperbolic  

Boundary  Value Proble m  
 

 

 

 

Department of  Mathematics, College of  Science, Mustansiriyah University, Baghdad, Iraq 

 

 

 

 

Abstract 

This paper is concerned with finding the approximation solution (APPS) of a certain 

type of nonlinear hyperbolic boundary value problem (NOLHYBVP).  The given BVP is 

written in its discrete (DI) weak form (WEF), and is proved that  it has a unique APPS, which 

is obtained via the mixed Galerkin finite element method (GFE) with implicit method 

(MGFEIM) that reduces the problem to solve the Galerkin nonlinear algebraic system  

(GNAS).  In this part, the predictor and the corrector technique (PT and CT) are proved 

convergent and are used to transform the obtained GNAS to linear (GLAS), then the GLAS is 

solved using the Cholesky method (ChMe). The stability and the convergence of the method 

are studied. The results are given by figures and shown the efficiency and accuracy for the 

method. 

 

Keywords:  nonlinear hyperbolic boundary value problem, Galekin finite element method, 

implicit method, convergence, stability. 

 

1. Introduction 

          Hyperbolic partial differential equations play a very important role as real life problems 

in many fields of sciences as in technology, fluid dynamics, optics, science and many others. 

In the past few decades, there have been many researchers interested in their study to solve 

boundary value problems in general and in particular NLHBVE. Many researchers have used 

different methods to solve the NLHBVE, Smiley studied in 1987, was used Eigen function 

methods to solve problems of nonlinear hyperbolic value at resonance [1]. In 1989, Chi, 

Wiener, and Shah used in the exponential growth of solutions of nonlinear hyperbolic 

equations [2], while in 2001 Minamoto used the existence and demonstration of the 

uniqueness of solutions [3]. In 2004, Krylovas, and Čiegis, used the numerical asymptotic 

averaging for weakly nonlinear hyperbolic waves [4]. 

In 2018, Ashyralyev and Agirseven solved NLHBVE with a time delay [5]. 

Ibn Al Haitham Journal for Pure and Applied Science 

Journal homepage: http://jih.uobaghdad.edu.iq/index.php/j/index 

 

Doi: 10.30526/34.2.2618 

Article history: Received 21, April ,2020, Accepted,23,June ,2020, Published in April 2021 

 

Jamil A. Ali Al-Hawasy   
jhawassy17@mustansiriyah.edu.iq 

Nuha Farhan Mansour 

hawasy20@yahoo.com 

mailto:jhawassy17@mustansiriyah.edu.iq
mailto:hawasy20@yahoo.com


  

120 

 

Ibn Al-Haitham Jour. for Pure & Appl. Sci. 34 (2) 2021 

 

      The specific element method has been studied by several researchers interested in this 

field, for example, in 2010 Bangerth and Rannacher touched on Galerkin's specific adaptation 

techniques for wave equation [6]. Whereas, in 2017, Al-Haq and Muhammad discussed 

numerical methods to solve LHYBVP by difference method and the method of the specified 

elements [7].  

In this paper, we are concerned the study of the APPS of the NOLHYBVP. The given 

BVP is written in its WEF, and in its discrete equation (DI) type. It is proved to have unique 

APPS. The APPS is obtained via the MGFEIM. The problem then reduces to solve the 

GNAS, then the PT and CT are proved convergent and are used to transform the GNAS to a 

GLAS. This GLAS is solved by using the ChMe. The stability and the convergence of the 

method are studied. A computer program is codding to find the numerical solution for the 

problem. The results are given by figures, and are shown the efficiency and accuracy for the 

method which is highly considered in this work. 

 

2. Description of The NOLHYBVP  

Let  𝐼 = [0, T], with  0 < 𝑇 < ∞ ,  𝜓 ⊂ ℝ2  be a bounded and open region with smooth 

boundary ∂𝜓 ,   𝜑 = 𝜓 × 𝐼 , Σ = 𝜕𝜓 × 𝐼 ,   then the NLHYPVP is given by: 

𝑤𝑡𝑡 −  Δ𝑤 + 𝑤 = ℎ(�⃗�, 𝑡, 𝑤),  in     𝜑                                                                                           

(1) 

𝑤(�⃗�, 0) =  𝑤 0 (�⃗�),    in    𝜓                                                                                                         

(2)   

𝑤𝑡 (�⃗�, 0)  =  𝑤
1 (�⃗�),   in     𝜓                                                                                                       

(3)  

𝑤(�⃗�, 𝑡) =  0 ,  on  Σ                                                                                                                     

(4) where  𝑤 = 𝑤(�⃗�, 𝑡) ∈  𝐻0
2(𝜓), ,  Δ𝑤 =  ∑

𝜕2𝑤

𝜕𝑥𝑖
2

2
𝑖=1   and  ℎ ∈  𝐿

2(𝜓)  is a given  function.   

Now, let 𝑉 = 𝐻0
1(𝜓)={ 𝜂:𝜂 ∈ 𝐻1(𝜓), 𝜂 = 0 on  ∂𝜓}, 𝑤𝑡 = 𝑝, then the WEFM of  (1-4) is:  

⟨ 𝑤𝑡𝑡  , 𝜂 ⟩  +  (∇𝑤, ∇ 𝜂)  +  (𝑤, 𝜂)  =  (ℎ(𝑤) , 𝜂 ) ,  ∀ 𝜂 ∈ 𝑉  are on  𝐼,                                    

(5)  

(𝑤(0), 𝜂)  =  (𝑤 0, 𝜂)    in  𝜓 ,  𝑤 0∈ 𝑉                                                                                       

(6) 

 (𝑝(0), 𝜂)  =  (𝑤1, 𝜂)    in  𝜓 , 𝑤1 ∈  𝐿2(𝜑)  ,                                                                            

(7) 

 

Definition (1),[8]: A point 𝑠∗ ∈ 𝐷 ⊂  ℝ2 is called a fixed point of the function 𝑦 ∶ 𝐷 ⟶  ℝ2 

, 

if 𝑦(𝑠∗) =  𝑠∗ . 

Definition 2,[8]: A function 𝑦: 𝐷 ⊂ ℝ2 ⟶ ℝ2 is called contractive on 𝐷 if for each 𝑑1, 𝑑2 ∈

𝐷: 

  ‖𝑦(𝑑2) − 𝑦(𝑑1)‖ ≤ 𝑎‖𝑑2 − 𝑑1‖, where 𝑎 ∈ (0,1).  

Theorem (3),[8]: A contractive function 𝑦 on a complete normed space 𝐷 has a unique fixed 

point 𝑠∗ in 𝐷. 

Theorem (4),[9]: Let  {𝑣𝑛}  be a bounded sequence in the space in 𝐿
∞(𝜓). Then, there exists 

a subsequence {𝑛′} and a function 𝑣0 ∈ 𝐿
∞(𝜓) such that, in 𝐿∞(𝜓) then  𝑣𝑛′  ⟶ 𝑣0. 

 



  

121 

 

Ibn Al-Haitham Jour. for Pure & Appl. Sci. 34 (2) 2021 

 

Assumptions 2.1:                                                                                                                       

 (i)  Let  κ1 and κ2  be two positive constants  such that the following are satisfied: 

    a)⎹(∇μ1 , ∇μ2)⎸≤  κ1 ∥ ∇μ1 ∥1 ∥ ∇μ2 ∥1  , ∀ μ1 , μ2 ∈ 𝑉  

    b) (∇μ , ∇μ ) ≥  κ2 ∥ ∇μ ∥1
2    ,    ∀ μ∈ 𝑉                                                                                     

(ii) The function ℎ is defined on  𝜑 × ℝ , continuous with respect to 𝑤𝑗
𝑛

 satisfies the 

following: 

    a)|ℎ(�⃗�, 𝑡, 𝑤)|  ≤  𝛽 (�⃗�, 𝑡) +  𝛿 |𝑤| where δ > 0, 𝑤 ∈  𝜑  and 𝛽 ∈  𝐿2(𝜑).  

    b) |ℎ(�⃗�, 𝑡, 𝑤1) − ℎ(�⃗�, 𝑡, 𝑤2)|  ≤ 𝐿|𝑤1 − 𝑤2|, where 𝐿 is a Lipchitz constant and 𝑤1, 𝑤2 ∈

 ℝ. 

 

3. Discretization of the Continuous Equation (COE): 

The WEF of ((5)- (7)) is discretized  by using the GFEME , let 𝜑 be divided  into sub 

regions  𝜑𝑖𝑗  = 𝜓𝑖
𝑛 × 𝐼𝑗

𝑛
 , let {𝜓𝑖

𝑛 }𝑖=1
𝑁(𝑛)

  be a triangulation of  �̅�   and  {𝐼𝑗
𝑛}𝑗=0  be a subdivision 

of the interval 𝐼 ̅into 𝑌(n) intervals, then  𝐼𝑗 =  𝐼𝑗
𝑛 ≔ [𝑡𝑗

𝑛 , 𝑡𝑗+1
𝑛 ] has the same length ∆𝑡 =  

𝑇

𝑌
 , 

also, let   𝑉𝑛 ⊂ V = 𝐻0
1(𝜓) be the space of piecewise affine functions in 𝜓. 

Now, the discrete equations (DES),  where ∀ 𝜂 ∈ 𝑉𝑛 are written  as follows:  

⟨𝑝𝑗+1
𝑛 − 𝑝𝑗

𝑛  , 𝜂 ⟩  +  Δt (∇𝑤𝑗+1
𝑛 , ∇ 𝜂) +  Δt (𝑤𝑗+1

𝑛 , 𝜂)  =  Δt (ℎ(𝑤𝑗+1
𝑛 ) , 𝜂)                                  

(8)  

 𝑤𝑗+1
𝑛 − 𝑤𝑗

𝑛  =  Δt 𝑝𝑗+1
𝑛

                                                                                                                

(9) (𝑤(0), 𝜂)  =  (𝑤 0, 𝜂)    in   𝜓                                                                                                     

(10) 

 (𝑝(0), 𝜂)  =  (𝑤1, 𝜂)    in  𝜓                                                                                                      

(11)       

where, 𝑤 0 ∈ 𝑉, 𝑤1 ∈ 𝐿2(𝜓), and  𝑤𝑗  
𝑛 = 𝑤 

𝑛(𝑥, 𝑡𝑗
𝑛 ), 𝑝𝑗

𝑛 = 𝑝 
𝑛(�⃗� , 𝑡𝑗

𝑛 ) ∈ 𝑉𝑛 , ∀  𝑗 = 0,1, … , 𝑌 −

1.   

4. The APPS of the NLHYBVP: 

To find the APPS �̅� 𝑛 = (𝑤0  
𝑛 , 𝑤1  

𝑛 , … , 𝑤𝑌  
𝑛 )for the DES (8-11), the MGFEIM is used through 

the following steps: 

(1) Let  { 𝜂𝑖 ∶  𝑖 = 1,2, … . 𝑁, 𝑤𝑖𝑡ℎ 𝜂𝑖 (�⃗�) = 0 , 𝑜𝑛 ∂𝜓 } be a basis of 𝑉𝑛, and by using the 

GFEME, let   �̅� 𝑛 (�⃗� , 𝑡𝑗
𝑛 ) (with   �̅�𝑡

𝑛(�⃗� , 𝑡𝑗
𝑛 ) = �̅�(�⃗� , 𝑡𝑗

𝑛 ))be an APPS    of (8-11) such that  

 �̅� 𝑛 (�⃗� , 𝑡𝑗
𝑛 ) =  ∑ 𝑟𝑘

𝑗𝑁
𝑘=1 𝜂𝑖  and �̅�

𝑛(�⃗� , 𝑡𝑗
𝑛 ) =  ∑ 𝑢𝑘

𝑗𝑁
𝑘=1 𝜂𝑖   ∀ 𝜂𝑖  ∈  𝑉𝑛,   

        where  𝑟𝑘
𝑗

= 𝑟𝑘 (𝑡𝑗
𝑛 ) ,and   𝑢𝑘

𝑗
= 𝑢𝑘 (𝑡𝑗

𝑛 )  are unknown constants  ∀𝑗 = 0,1, … , 𝑌 − 1.  

(2)  Using the APPs in (8-11) to get ,   ∀𝑗 = 0,1, … , 𝑌 − 1:   

     (𝑀 + (Δt)2𝑄 )𝑅𝑗+1 = 𝑀𝑅𝑗 + (Δt) 𝑀𝑈𝑗 + (Δt)2 �⃗⃗� (𝑡𝑗
𝑛 ,  �⃗� 𝑇 𝑅𝑗+1)                                 (12) 

     𝑈𝑗+1 =
1

Δt
 (𝑅𝑗+1 −  𝑅𝑗)                                                                                                    (13) 

     𝑀𝑅0 =  𝑠0                                                                                                                        (14) 

     𝑀𝑈0 =  𝑠1                                                                                                                         (15) 

     where,  𝑀 = (𝑚𝑖𝑘 )𝑁×𝑁  , 𝑚𝑖𝑘 = (𝜂𝑘 , 𝜂𝑖 ), 𝑄 = (𝑞𝑖𝑘 )𝑁×𝑁,   𝑞𝑖𝑘 =  (∇𝜂𝑘 , ∇𝜂𝑖 ), �⃗⃗�  = (𝐿𝑖 )𝑁×1 ,      

      𝐿𝑖 = (ℎ( �⃗�
𝑇 𝑅𝑗+1), 𝜂𝑖 ) ,  𝑅𝑁×1

𝑗
= (𝑟1

𝑗
, 𝑟2

𝑗
 , … , 𝑟𝑁

𝑗
)𝑇 ,𝑈𝑁×1

𝑗
= (𝑢1

𝑗
, 𝑢2

𝑗
 , … , 𝑢𝑁

𝑗
)𝑇,   𝑠0 =

(𝑠𝑖
0)𝑁×1,     

       𝑠𝑖
0 = ( 𝑤 0, 𝜂𝑖 ),  𝑠

1 = (𝑠𝑖
1)𝑁×1  and  𝑠𝑖

1 = ( 𝑤1, 𝜂𝑖 ) , for each 𝑖 , 𝑘 = 1,2, … , 𝑁. 



  

122 

 

Ibn Al-Haitham Jour. for Pure & Appl. Sci. 34 (2) 2021 

 

(3)  System (12-15), is GNAS and has a unique solution. To solve it, we find at first 𝑅0 and 

𝑈0 from solving (14) and (15) respectively, then, the PT and the CT are utilized to solve 

(12) for each 𝑗( 𝑗 = 0,1, … , 𝑌 − 1) as follows:  

In the PT we suppose 𝑅𝑗+1 = 𝑅𝑗 in the components of  �⃗⃗�  in the R.H.S of (12), then it turn 

to a GLAS, which is solved to get the predictor solution 𝑅𝑗+1, then in the CT we resolve 

(12) with setting �̅�𝑗+1 = 𝑅𝑗+1 (in the components of  �⃗⃗�  of the R.H.S of it)  to get the 

corrector solution  𝑅𝑗+1, finally substituting  𝑅𝑗+1 in (13) to get 𝑈𝑗+1, we can  repeat this 

procedure if we want more than one time. This reputation can be expressed as follows: 

 

(𝑤
𝑗+1
(𝑙+1)

, 𝜂𝑖 ) + (Δt)
2(∇𝑤𝑗+1

(𝑙+1)
, ∇𝜂𝑖 ) + (Δt)

2 (𝑤𝑗+1
(𝑙+1)

, 𝜂𝑖 ) = (𝑤𝑗
𝑛, 𝜂𝑖 ) +  Δt (𝑝𝑗

𝑛, 𝜂𝑖 )  

  + (Δt)2ℎ(𝑡𝑗
𝑛 , 𝑤𝑗+1

(𝑙)
) , 𝜂𝑖 )                                                                                                 (16)  

𝑝
𝑗+1
(𝑙+1)

=  
(𝑤

𝑗+1
(𝑙+1)

−𝑤𝑗
𝑛)  

Δt
                                                                                                              (17)    

Equation (17) tells us the iterative method depending on just 𝑤
𝑗+1
(𝑙+1)

. Thus, equation (16) is 

reformulated as  𝑤 (𝑙+1) =  𝛿(𝑤 (𝑙+1)) , where  𝑙 is the number of the iterations. And this led us 

to the following theorem. 

Theorem (5): For any fixed point, the DES (8)-(11), and for Δ  sufficiently small, has a 

unique solution 𝑤 𝑛 = (𝑤0
𝑛, 𝑤1

𝑛, … . . , 𝑤𝑁
𝑛) and the sequence of the corrector solutions 

converges on ℝ. 

proof: Let  𝑤 (𝑙+1) = (𝑤0
(𝑙+1)

, 𝑤1
(𝑙+1)

, … . . , 𝑤𝑁
(𝑙+1)

 ) and  

 𝑤
(𝑙+1)

= (𝑤0
(𝑙+1)

, 𝑤1
(𝑙+1)

, … . . , 𝑤𝑁
(𝑙+1)

 ) where  𝑤 (𝑙+1)  and  𝑤
(𝑙+1)

 are   

solutions of equation (16). This means, 

(𝑤
𝑗+1
(𝑙+1)

, 𝜂𝑖 )   +  (Δt)
2(∇𝑤𝑗+1

(𝑙+1)
, ∇  𝜂𝑖 ) + (Δt)

2  (𝑤
𝑗+1
(𝑙+1)

, 𝜂𝑖 )  =  (𝑤𝑗
𝑛, 𝜂𝑖 )  +    Δt (𝑝𝑗

𝑛, 𝜂𝑖 )   

+ (Δt)2 (ℎ(𝑡𝑗
𝑛 , 𝑤𝑗+1

(𝑙)
) , 𝜂𝑖  )                                                                                                   (18)      

and 

(𝑤
𝑗+1

(𝑙+1)
, 𝜂𝑖 )   + (Δt)

2(∇ 𝑤
𝑗+1

(𝑙+1)
, ∇ 𝜂𝑖 ) + (Δt)

2 (𝑤
𝑗+1

(𝑙+1)
, 𝜂𝑖 ) =    (𝑤𝑗

𝑛, 𝜂𝑖 )  +   Δt (𝑝𝑗
𝑛, 𝜂𝑖 )     

+ (Δt)2(ℎ (𝑡𝑗
𝑛 , 𝑤

𝑗+1

(𝑙)
) , 𝜂𝑖  )                                                                                                    (19)                

subtracting )19) from (18) then setting  𝜂𝑖 = (𝑤𝑗+1
(𝑙+1)

−  𝑤𝑗+1
(𝑙+1)

) in the obtained equation , we 

get that                                                      

 (𝑤
𝑗+1

(𝑙+1)
− 𝑤𝑗+1

(𝑙+1)
,  𝑤𝑗+1

(𝑙+1)
− 𝑤𝑗+1

(𝑙+1)
) + (Δt)2(∇ 𝑤

𝑗+1

(𝑙+1)
− ∇𝑤𝑗+1

(𝑙+1)
,  ∇ 𝑤𝑗+1

(𝑙+1)
− ∇𝑤𝑗+1

(𝑙+1)
) + 

(Δt)2(𝑤
𝑗+1

(𝑙+1)
− 𝑤𝑗+1

(𝑙+1)
,  𝑤𝑗+1

(𝑙+1)
− 𝑤𝑗+1

(𝑙+1)
) = (Δt)2(ℎ ( 𝑤𝑗+1 

(𝑙)
) − ℎ( 𝑤𝑗+1 

(𝑙)
),  𝑤𝑗+1

(𝑙+1)
− 𝑤𝑗+1

(𝑙+1) 
) 

(20)                       

From Assumptions 2.1 (ib) the 2𝑛𝑑  and 3𝑟𝑑  terms in the L.H.S of equation (20) are positive, 

and applying Assumption 2.1 (iib) on  ℎ in R.H.S of equation (20), and by using the Cauchy 

Schwarz inequality on this side, we get 

 ‖𝛿 ( 𝑤𝑗+1 
(𝑙)

) −  𝛿( 𝑤𝑗+1 
(𝑙)

) ‖
0

= ‖𝑤𝑗+1
(𝑙+1)

−  𝑤𝑗+1
(𝑙+1)

‖
0

≤  𝜆 ‖𝑤𝑗+1
(𝑙)

−  𝑤𝑗+1
(𝑙)

‖
0
                     (21) 

where  𝜆 =  (Δt)2𝐿 < 1, for  sufficiently small  Δ𝑡.                             

which implies  that 𝛿 is contractive, also since { 𝑤 (𝑙)} ∈ ℝ ∀ 𝑙,  that  𝛿(𝑤 (𝑙+1)) =  𝑤 (𝑙+1)  ∈

ℝ ∀ 𝑙,  𝑖. 𝑒  𝛿(𝑤) ∈ ℝ ,  hence ,by theorem 3  the sequence  { 𝑤 (𝑙)}  converges to a point in ℝ.  



  

123 

 

Ibn Al-Haitham Jour. for Pure & Appl. Sci. 34 (2) 2021 

 

 

5. Stability:  

Lemma (6): If Δ  is sufficiently small, then ∀ 𝑗 = 0,1, … , 𝑌 

  ‖𝑤𝑗 
𝑛 ‖1

2  ≤ �̅�, ‖𝑝𝑗 
𝑛 ‖0

2  ≤ �̅� ,   ∑  ‖ 𝑤𝑗+1
𝑛 − 𝑤𝑗

𝑛 ‖1 
2𝑌−1

𝑗=0  ≤ �̅� , and  ∑  ‖ 𝑝𝑗+1
𝑛 − 𝑝𝑗

𝑛 ‖0 
2𝑌−1

𝑗=0  ≤ �̅�   

where �̅� refers to a various constants. 

proof: Let  𝜂 = 𝑝𝑗+1
𝑛 

  substituting in equation (8), and rewriting the first term in the L.H.S of 

the obtained equation, we get 

‖ 𝑝𝑗+1
𝑛 ‖0

2  −  ‖ 𝑝𝑗
𝑛 ‖0

2  +  ‖ 𝑝𝑗+1
𝑛 − 𝑝𝑗

𝑛 ‖0 
2  +  Δt (∇𝑤𝑗+1

𝑛 , ∇𝑝𝑗+1
𝑛 )  +  Δt (𝑤𝑗+1

𝑛 , 𝑝𝑗+1
𝑛 )  =   

Δt (ℎ(𝑤𝑗+1
𝑛 ) , 𝑝𝑗+1

𝑛 )                                                                                                                (22)  

Since,   

 𝛥𝑡 [(∇𝑤𝑗+1
𝑛  ,   ∇𝑝𝑗+1

𝑛 )  +  (𝑤𝑗+1
𝑛 , 𝑝𝑗+1

𝑛 )]  =   
1

2
 [ (∇𝑤𝑗+1

𝑛 −  ∇𝑤𝑗
𝑛 , ∇𝑤𝑗+1

𝑛 −  ∇𝑤𝑗
𝑛)  +                    

𝑤𝑗+1
𝑛 −  𝑤𝑗

𝑛  , 𝑤𝑗+1
𝑛 −  𝑤𝑗

𝑛 )  +  (∇𝑤𝑗+1
𝑛  , ∇𝑤𝑗+1

𝑛 )  + (𝑤𝑗+1
𝑛  , 𝑤𝑗+1

𝑛 )    −  (∇𝑤𝑗
𝑛 , ∇𝑤𝑗

𝑛) −

 (𝑤𝑗
𝑛 , 𝑤𝑗

𝑛) ]   

By substituting above equality in the L.H.S in equation (22), summing both sides of the 

obtained equality, for  𝑗 = 0 𝑡𝑜 𝑗 = 𝑙 − 1,  then set  c = max  (1,
𝑘2

2
 ),  we get  

 𝑐‖ 𝑝𝑙
𝑛 ‖0

2 +  𝑐 ∑  ‖ 𝑝𝑗+1
𝑛 − 𝑝𝑗

𝑛 ‖0 
2𝑙−1

𝑗=0 +  𝑐‖𝑤𝑙
𝑛 ‖1

2  + 𝑐 ∑  ‖  𝑤𝑗+1
𝑛 − 𝑤𝑗

𝑛 ‖1 
2𝑙−1

𝑗=0   ≤   ‖ 𝑝0
𝑛 ‖0

2 +   
𝑘2

2
 

‖ 𝑤𝑙
𝑛 ‖1

2   +  ∑  𝑙−1𝑗=0 Δt (ℎ(𝑤𝑗+1
𝑛 ) , 𝑝𝑗+1

𝑛 )                                                                                  (23)  

 

Now, using the assumptions on ℎ and then by the Cauchy Schwarz inequality, to get    

 │(ℎ(𝑤𝑗+1
𝑛 ) , 𝑝𝑗+1

𝑛 )│ ≤ ‖ 𝛽𝑗
 ‖0

2 + 𝛿‖ 𝑤𝑗+1
𝑛 ‖1

2 + 𝛿̅‖ 𝑝𝑗+1
𝑛 ‖0

2 , 𝛿̅ =  𝛿 + 1                               (24) 

since ‖ 𝑤𝑗+1
𝑛 ‖1

2  = 2‖ 𝑤𝑗+1
𝑛 − 𝑤𝑗

𝑛 ‖1 
2 + 2‖ 𝑤𝑗

𝑛 ‖1
2                                                                   (25) 

and    ‖ 𝑝𝑗+1
𝑛 ‖0

2  = 2‖ 𝑝𝑗+1
𝑛 − 𝑝𝑗

𝑛 ‖0
2 + 2‖ 𝑝𝑗

𝑛 ‖0
2                                                                      (26) 

Substituting (25) and (26) in inequality (23), and assume that 𝑑 =  max(2𝛿, 2𝛿̅), to get  

 𝑐‖ 𝑝𝑙
𝑛 ‖0

2  +  (𝑐 − 𝑑Δt) ∑  ‖ 𝑝𝑗+1
𝑛 − 𝑝𝑗

𝑛 ‖0
2𝑙−1

𝑗=0  +  𝑐‖𝑤𝑙
𝑛 ‖1

2 + (𝑐 −

𝑑Δt) ∑  ‖ 𝑤𝑗+1
𝑛 − 𝑤𝑗

𝑛 ‖1 
2𝑙−1

𝑗=0 ≤    

‖ 𝑝0
𝑛 ‖0

2 +   
𝑘2

2
 ‖𝑤𝑙

𝑛 ‖1
2 +  ‖𝛽‖𝑄

2  +  𝑑(Δt) ∑ ‖𝑤𝑗
𝑛 ‖1

2 𝑙−1𝑗=0 + 𝑑(Δt) ∑ ‖𝑝𝑗
𝑛 ‖0

2 .𝑙−1𝑗=0                    (27) 

 Now, let ∆t < 𝑐 𝑑⁄   then the  2𝑛𝑑   and  4𝑡ℎ   terms in the R.H.S of (27) are positives, by using 

the discrete Gronwall’s (DGs) inequality [10], one obtains    

c(‖ 𝑝𝑙
𝑛 ‖0

2 +  ‖ 𝑤𝑙
𝑛 ‖1

2) ≤  𝑎𝑒
∑ 𝑑(Δt) 𝑙−1𝑗=0 = 𝑎𝑒𝑙𝑑(Δt) ≤ 𝑏  ,  

which gives that 

‖ 𝑤𝑙
𝑛 ‖1

2  ≤ 𝑑1 =  
𝑏

𝑐
 , and  ‖ 𝑝𝑙

𝑛 ‖0
2  ≤ 𝑑1,  for any arbitrary index 𝑙.  

Hence, ‖ 𝑤𝑗
𝑛 ‖1

2  ≤  𝑑1 and ‖ 𝑝𝑗
𝑛 ‖0

2  ≤  𝑑1, for each   𝑗 = 0,1, … . . , 𝑌 − 1. 

Therefore 

  (Δt)𝑑 ∑ ‖𝑤𝑗
𝑛 ‖1

2 𝑌−1𝑗=0 +  (Δt)𝑑 ∑ ‖𝑝𝑗
𝑛 ‖1

2 𝑌−1𝑗=0  ≤ 2𝑑1 𝑑 Δt Y =  2cT =  �̅� .   

We back to (27) substituting 𝑙 = 𝑌, the 1st and the 3rd term in the L.H.S are positives, then we 

use the above results in the R.H.S. of it , keeping in mind the first three terms in this side that  

are bounded (from the above steps), to obtain    

∑ ‖ 𝑤𝑗+1
𝑛 − 𝑤𝑗

𝑛 ‖1
2 𝑌−1𝑗=0 ≤ �̅�                                                                                          (28a) 

∑ ‖ 𝑝𝑗+1
𝑛 − 𝑝𝑗

𝑛 ‖0
2 𝑌−1𝑗=0  ≤ �̅�                                                                                                   (28b) 

6. Convergence: 



  

124 

 

Ibn Al-Haitham Jour. for Pure & Appl. Sci. 34 (2) 2021 

 

The following definitions for the functions "almost everywhere on I " are useful in the proof 

of next theorem, so let  

𝑤−
𝑛 (𝑡)  ∶=  𝑤𝑗

𝑛 , 𝑡 ∈  𝐼𝑗
𝑛 , ∀ 𝑗 = 0,1, … . , 𝑌,   

 𝑤+
𝑛 (𝑡)  ∶=  𝑤𝑗+1

𝑛   ,   𝑡 ∈  𝐼𝑗
𝑛 , ∀ 𝑗 = 0,1, … . , 𝑌 − 1,     

 𝑝+
𝑛 (𝑡)  ∶=  𝑝𝑗+1

𝑛 , 𝑡 ∈  𝐼𝑗
𝑛 , ∀ 𝑗 = 0,1, … . , 𝑌 − 1,   

𝑝−
𝑛 (𝑡𝑗

𝑛 )  ∶=  𝑝𝑗
𝑛 , 𝑡𝑗

𝑛 ∈  𝐼𝑗
𝑛 , ∀ 𝑗 = 0,1, … . , 𝑌,   

Also, Let  𝑤^
𝑛 (𝑡)  ∶=  𝑤𝑗

𝑛 
 be an affine function on each 𝐼𝑗

𝑛 
,   ∀ , 𝑗 = 0,1, … . , 𝑌, and 

  𝑝^
𝑛 (𝑡)  ∶=  𝑝𝑗

𝑛  
, be an affine function on each 𝐼𝑗

𝑛 
,   ∀ , 𝑗 = 0,1, … . , 𝑌. 

Theorem (7): The discrete solutions 𝑤−
𝑛 (𝑡) , 𝑤+

𝑛 (𝑡) , and 𝑤^
𝑛 (𝑡)  are converges strongly in 

𝐿2(𝜑), where  𝑛 ⟶  ∞. 

proof:  we start with using lemma (6) we have for any  𝑗 = 0,1, … . , 𝑌 that 

 

‖𝑤𝑗 
𝑛 ‖1

2  ≤ �̅�   and ‖𝑝𝑗 
𝑛 ‖0

2  ≤ �̅�,  then   

‖𝑤− 
𝑛 ‖

𝐿2(𝐼,𝑉)
2  , ‖ 𝑤+ 

𝑛 ‖
𝐿2(𝐼,𝑉)
2 , ‖𝑤^ 

𝑛 ‖
𝐿2(𝐼,𝑉)
2 , ‖𝑝− 

𝑛 ‖
𝐿2(𝜑)
2  , ‖𝑝+ 

𝑛 ‖
𝐿2(𝜑)
2  , and ‖𝑝^

𝑛‖
𝐿2(𝜑)
2

  are 

bounded. 

From (28a), we have  

 Δt ∑ ‖ 𝑤𝑗+1
𝑛 − 𝑤𝑗

𝑛 ‖0
2 𝑌−1𝑗=0 ≤ Δt�̅� ⟶ 0, as Δt ⟶ 0,  

gives 

�̅�+
𝑛 ⟶  �̅�−

𝑛   strongly (ST) in  𝐿2(𝐼, 𝑉) and then in 𝐿2(𝜑).                                                      

(29a) 

by the same way from (28b), we get that , 

 𝑝+
𝑛 ⟶ 𝑝− 

𝑛  ST in 𝐿2(𝜑)                                                                                                        (29b) 

Then by theorem 3, there exist subsequences of {𝑤− 
𝑛 }, {𝑤+

𝑛}, {𝑤^
𝑛 }), and of ({𝑝− 

𝑛 }, {𝑝+
𝑛}, 

{𝑝^
𝑛 },) use again the same notations which converge weakly to some 𝑤 in  𝐿2(𝐼, 𝑉), to some 𝑝 

in 𝐿2(𝜑), i.e.  

 𝑤−
𝑛 ⟶ 𝑤 , 𝑤+

𝑛 ⟶ 𝑤, 𝑤^
𝑛 ⟶ 𝑤  weakly  in  𝐿2(𝐼, 𝑉)   

            𝑝−
𝑛 ⟶ 𝑝 , 𝑝+

𝑛 ⟶ 𝑝, 𝑝^
𝑛 ⟶ 𝑝  weakly  in  𝐿2(𝜑), 

N this point the first compactness theorem[9] is used , to get that  𝑤^
𝑛 ⟶ 𝑤 ST in 𝐿2(𝜑),then 

 𝑤+
𝑛 ⟶ 𝑤 and 𝑤−

𝑛 ⟶ 𝑤 ST  in 𝐿2(𝜑). 

 

Now, let   {𝑉𝑛 }𝑛=1
∞  be a sequence of subspaces of 𝑉, where  𝑉𝑛 is as defined above. Then by 

using the Galerkin approach, for each 𝜂 ∈ 𝑉, there exists a sequence {𝜂𝑛 }, with  𝜂𝑛 ∈  𝑉𝑛    

for each 𝑛,  such that  𝜂𝑛 ⟶ 𝜂 ST in 𝐿
2(𝜑).  

Consider that 𝜉(𝑡) ∈  𝐶2[0, 𝑇], for which   𝜉(𝑇) =  𝜉′(𝑇) = 0 and  𝜉(0) =  𝜉′(0) ≠

0 , let  𝜉𝑛(𝑡)  continuous piecewise(CP) interpolation  of    𝜉(𝑡) with  respect to 𝐼𝑗
𝑛

, and let  

𝜁 = 𝜂 𝜉(𝑡), with  𝜁𝑛 =  𝜂𝑛 𝜉
𝑛 (𝑡), with  

𝜁− 
𝑛  ∶=  𝜂𝑛 𝜉−

𝑛(𝑡) , 𝑡 ∈  𝐼𝑗
𝑛 , 𝑗 = 0,1, … . , 𝑌 − 1 , 𝜂𝑛 ∈  𝑉𝑛 ,  

𝜁+ 
𝑛  ∶=  𝜂𝑛 𝜉+

𝑛(𝑡), 𝑡 ∈  𝐼𝑗
𝑛 ,    𝑗 = 0,1, … . , 𝑌 − 1 , 𝜂𝑛 ∈  𝑉𝑛  , 

𝜁^ 
𝑛  ∶=  𝜂𝑛 𝜉

𝑛(𝑡), 𝑡 ∈  𝐼 ,    𝜂𝑛 ∈  𝑉𝑛 ,     

Setting  𝜂 =  𝜁𝑗+1 
𝑛 

 in equation (8), and summing both sides of the obtained equation for 𝑗 = 0, 

to 𝑗 = 𝑌 − 1 , then using discrete integrating by parts (DIBP) for the 1st term in the L.H.S.,  

once can get  that     



  

125 

 

Ibn Al-Haitham Jour. for Pure & Appl. Sci. 34 (2) 2021 

 

− ∫ (𝑝− 
𝑛  , (𝜁^

𝑛 𝑇

0
)′) 𝑑𝑡  +  ∫ [(

𝑇

0
∇ 𝑤+ 

𝑛 , ∇ 𝜁+ 
𝑛 ) + ( 𝑤+ 

𝑛 , 𝜁+ 
𝑛 )] 𝑑𝑡 = ∫ (

𝑇

0
ℎ(𝑡−  

𝑛 , 𝑤+ 
𝑛 ), 𝜁+ 

𝑛 ) 𝑑𝑡 +

 (𝑝0 
𝑛 , 𝜂𝑛 ) 𝜉(0)                                                                                                                             

(30) 

On the other hand, from (9), once has 

   ((𝑤^
𝑛 )′, 𝜂𝑛 )(𝜉

𝑛)′ = (𝑝+ 
𝑛 , 𝜂𝑛 )(𝜉

𝑛)′ 

Integrating both sides on [0, 𝑇], then using DEBP for the 1st term in the L.H.S., to obtain  

− ∫ (𝑤+
𝑛 𝑇

0
, 𝜂𝑛 ) (𝜉

𝑛(𝑡))′′𝑑𝑡 =  ∫ (𝑝+ 
𝑛  

𝑇

0
, 𝜂𝑛 )(𝜉

𝑛)′𝑑𝑡 + (𝑤0 
𝑛 , 𝜂𝑛 )(𝜉

𝑛(0))′                             

(31) 

Now, since 

 𝜉𝑛(𝑡)  ⟶  𝜉(𝑡)   in 𝐶(𝐼) ⊂ 𝐿2(𝐼), 𝜂𝑛  ⟶ 𝜂  ST  in 𝐿
2(𝐼, 𝑉) and in  𝐿2(𝜓), then, we have 

𝜁+ 
𝑛  = 𝜂𝑛 𝜉+

𝑛  ⟶ 𝜂 𝜉 = 𝜁    ST in 𝐿2(𝐼, 𝑉) and in 𝐿2(𝜑), 𝜂𝑛 𝜉
𝑛(0)   ⟶ 𝜂 𝜉(0)  ST  in 𝐿2(𝜑),  

 (𝜁^
𝑛 )′ =   𝜂𝑛 𝜉

′𝑛  ⟶ 𝜂 𝜉′ =  𝜂𝜁′ ST  in 𝐿2(𝐼, 𝑉). 

And since , 𝑡−
𝑛  ⟶ 𝑡   ST  in 𝐿∞(𝐼),  𝑤+

𝑛 , 𝑤−
𝑛,  𝑤^

𝑛 ⟶ 𝑤  ST   in 𝐿2(𝜑),  𝑤0
𝑛  ⟶  𝑤 0  ST  in  

𝑉 and  𝑝0
𝑛  ⟶  𝑤1 ST  in 𝐿2(𝜓). 

Now, from assumptions ℎ, and the above convergences, one  can passage to the limit in (30) 

and in (31), to obtain  

 − ∫ (𝑝
𝑇

0
, 𝜂) 𝜉′𝑑𝑡 + ∫ [(∇𝑤

𝑇

0
, ∇𝜂) + (𝑤 , 𝜂)] 𝜉 𝑑𝑡 = ∫ (

𝑇

0
ℎ(𝑡, 𝑤), 𝜂) 𝜉 𝑑𝑡 + (𝑤1, 𝜂) 𝜉(0)   (32)     

and  

− ∫ (𝑤
𝑇

0
, 𝜂)𝜉′′(𝑡) 𝑑𝑡 = ∫ (𝑝 

𝑇

0
, 𝜂)𝜉′(𝑡)𝑑𝑡 + (𝑤 

0 , 𝜂)𝜉′(0)                                                    (33) 

The following cases appear: 

Case (1) : Consider  𝜉(𝑡) ∈  𝐶2[0, 𝑇], such that 𝜉(𝑇) =  𝜉′(𝑇) = 𝜉(0) =  𝜉′(0) = 0, by 

setting  𝜉′(0) = 0 in equation (32) and 𝜉(0) = 0 in (33), then we use IBP for the 1st term of 

each one of the obtained equation, one gets respectively 

∫ (𝑤𝑡
𝑇

0
, 𝜂)𝜉′(𝑡) 𝑑𝑡 =  ∫ (𝑝

𝑇

0
, 𝜂)𝜉′(𝑡)𝑑𝑡 ⟹  𝑤𝑡 = 𝑝 ,   

∫ (𝑤𝑡𝑡
𝑇

0
, 𝜂) 𝜉  𝑑𝑡 + ∫ [(∇𝑤

𝑇

0
, ∇𝜂) + (𝑤 , 𝑣)] 𝜉 𝑑𝑡 = ∫ (

𝑇

0
ℎ(𝑡, 𝑤), 𝜂)𝜉 𝑑𝑡,                              (34) 

Thus  

(𝑤𝑡𝑡  , 𝜂) + (∇w, ∇𝜂) + (𝑤 , 𝜂)  =   (ℎ(𝑡, 𝑤), 𝜂) , 𝜂 ∈ 𝑉 a. e. on 𝐼.   

 

Case (2) : Consider 𝜉(𝑡) ∈  𝐷[0, 𝑇],  𝜉(0) ≠  0, 𝜉(𝑇) = 0 and use   IBP   the first term in the 

L.H.S of (34),  once get  that 

− ∫ (𝑤𝑡
𝑇

0
, 𝜂) 𝜉′𝑑𝑡 + ∫ [(∇𝑤

𝑇

0
, ∇𝜂) + (𝑤 , 𝜂)] 𝜉 𝑑𝑡 = ∫ (

𝑇

0
ℎ(𝑡, 𝑤), 𝜂)𝜉 𝑑𝑡 + 𝑤𝑡 (0), 𝜂)𝜉(0)  

(35) 

Setting 𝑝 =  𝑤𝑡 in (32), subtracting the resulting equation from (35), to get  

(𝑤𝑡 (0) , 𝜂)𝜉(0) = (𝑤
1, 𝜂)𝜉(0)  ⟹  (𝑤𝑡 (0) , 𝜂) = (𝑤

1, 𝜂) , for each 𝜂  

then  𝑤𝑡 (0)  =  𝑤
1(0). 

 

Case (3): Consider  𝜉(𝑡) ∈  𝐷[0, 𝑇], with 𝜉′(0) ≠  0, 𝜉(0) = 0, and  𝜉(𝑇) =  𝜉′(𝑇) = 0. 

Using twice the IBP for the 1st term in the L.H.S. of (34), to obtain  

∫ (𝑤
𝑇

0
, 𝜂)𝜉′′𝑑𝑡 + ∫ [(∇𝑤

𝑇

0
, ∇𝜂) + (𝑤 , 𝜂)]𝜉 𝑑𝑡 = ∫ (

𝑇

0
ℎ(𝑡, 𝑤), 𝜂)𝜉 𝑑𝑡 − (𝑤(0), 𝜂)𝜉′(0)    (36) 

Rewritten (33), in the following form 

− ∫ (𝑝
𝑇

0
, 𝜂)𝜉′(𝑡) 𝑑𝑡 =  ∫ (𝑤

𝑇

0
, 𝜂)𝜉′′(𝑡)𝑑𝑡 + (𝑤 

0 , 𝜂)𝜉′(0)                                                   (37) 



  

126 

 

Ibn Al-Haitham Jour. for Pure & Appl. Sci. 34 (2) 2021 

 

Substituting (37) in (32), and using 𝜉(0) = 0, then subtracting the resulting equation from 

(36) to get 

 (𝑤(0) , 𝜂)𝜉′(0) = (𝑤 0, 𝜂)𝜉′(0)  ⟹  (𝑤(0) , 𝜂) = (𝑤 0, 𝜂) for each 𝜂 , then  𝑤(0) = 𝑤 0(0) ,  

Thus limit point  𝑤 is a solution to the WEF for the COE. 

 

7. Cholesky factorization 

     Cholesky method is used using to solve  GLAS with conditions that the coefficient matrix 

𝐴 must be a symmetric and positive definite. In this method the matrix 𝐴 can be factorized 

into the product of an Upper triangular matrix 𝐿 and Lower triangular matrix   𝐿𝑇  [11],  and 𝐿 

calculates as follows: 

    𝑓𝑜𝑟 𝑖 = 1,2, … , 𝑛      𝑡ℎ𝑒𝑛      𝑙𝑖𝑖 =  (𝑎𝑖𝑖 − ∑ 𝑙𝑘𝑖
2  

𝑖−1

  𝑘=1

)
1
2     

  𝑓𝑜𝑟 𝑗 = 𝑖 + 1, … , 𝑛.    𝑡ℎ𝑒𝑛      𝑙𝑖𝑗 = ( 𝑎𝑖𝑗 − ∑ 𝑙𝑘𝑖 

𝑖−1

𝑘=1

𝑙𝑘𝑗  ) / 𝑙𝑖𝑖   

8. Numerical Examples:  

     The problems in the following examples are coded by Mat lap soft.  

 Example 1: Consider the following NLHBVP:     

 𝑤𝑡𝑡   −  Δ 𝑤 + 𝑤 = ℎ(�⃗�, 𝑡, 𝑤),   �⃗� = (𝑥, 𝑦) , 𝜑 =  𝜓 ×  𝐼 , 𝜓 = (0,1) × (0,1), 𝐼 = [0,1] 

 𝑤(�⃗�, 0)  =  𝑥𝑦(1 − 𝑥)(1 − 𝑦), in    𝜓                                                                                            

 𝑤𝑡 (�⃗�, 0)  =  𝑤
1 (�⃗�) ,  in     𝜓                                                                                                           

𝑤(�⃗�, 𝑡) = 0  , on   ∑ = 𝜕𝜓 × 𝐼 

where ℎ(�⃗�, 𝑡, 𝑤) =  
1

2
(𝑥𝑦 − 𝑥𝑦2 − 𝑦𝑥2 + 𝑥2𝑦2)√𝑐𝑜𝑠2 [1 − 2sin(𝑥𝑦 − 𝑥𝑦2 − 𝑦𝑥2 +

𝑥2𝑦2)√𝑐𝑜𝑠𝑡 ] + 2(𝑦 + 𝑥 − 𝑦2 − 𝑥2)√𝑐𝑜𝑠𝑡 + (𝑥𝑦2 − 𝑦𝑥2 − 𝑥𝑦 − 𝑥2𝑦2) 𝑠𝑖𝑛2𝑡 4 √cos (𝑡)
3

⁄ +

𝑤 𝑠𝑖𝑛𝑤 

and the exact solution (EXS)of the problem is 𝑤(�⃗�, 𝑡) =  𝑥𝑦(1 − 𝑥)(1 − 𝑦)√cos (𝑡) .   

The MGFEIM is utilized to solve this problem with = 9 ,  𝑌 = 20 and 𝑇 = 1, the results are 

shown in figure 1. (a) the APPS , and figure 1.(B) shown the EXS    at  �̂� = 0.5. 

                                         
                                            Figure1. (a) shows the APS and  (b) shows the   EXS 

 

Example 2: Consider the following NLHYBVP :     

 𝑤𝑡𝑡 −  Δ 𝑤 + 𝑤 = ℎ(�⃗�, 𝑡, 𝑤)    where   �⃗� = (𝑥, 𝑦)  

 𝑤(�⃗�, 0) =  (𝑥 − 1)(1 − 𝑦) sin(𝑥𝑦)   in    𝜓                                                  



  

127 

 

Ibn Al-Haitham Jour. for Pure & Appl. Sci. 34 (2) 2021 

 

  𝑤𝑡  (�⃗�, 0)  =  𝑤
1 (�⃗�) , in     𝜓                                                                           

 𝑤(�⃗�, 𝑡) = 0  , on   ∑ = 𝜕𝜓 × 𝐼                                                     

 ℎ(�⃗�, 𝑡, 𝑤) = 2(𝑦 − 𝑥 − 𝑦2 + 𝑥2)√𝑒𝑡
2

cos(𝑥𝑦) + (1 − 𝑥 − 𝑦 + 𝑥𝑦)√𝑒𝑡
2

sin(xy)  

[𝑥2 + 𝑦2 − sin ((1 − 𝑥 − 𝑦 + 𝑥𝑦)√𝑒𝑡
2

sin(xy))] + 𝑤 sin (𝑤). 

and the EXS is 𝑤(�⃗�, 𝑡) =  (𝑥 − 1)(1 − 𝑦)√𝑒𝑡
2

sin(𝑥𝑦) . 

The MGFEIM is utilized to solve this problem with = 9 ,  𝑌 = 20 and 𝑇 = 1, the results are 

shown in Figure 2. (a) the APPS , and Figure 2.(b) shown the EXS    at  �̂� = 0.5. 

         
                        Figure2. (a) shows the APS and  (b) shows the   EXS 

 

9. Conclusions 

The MGFEIM is used successfully to solve the DI of the WEF of a certain type of 

NOLHYBVP. The existence theorem of a unique convergent APP is proved. The convergent 

of the PT and CT which are used to solve the GNAS that is obtained from applying the 

MGFEIM, is proved and  the ChMe which is used inside these technique is highly efficient 

for solving large GAS. The DI of the WEF is proved itis stable and convergent. The results 

are given by figures and show the efficiency and accuracy for the method. 

 

References  

1. Smiley, M. W. Eigenfunction methods and nonlinear hyperbolic boundary value problems 

at resonance. Journal of Mathematical Analysis and Applications .1987, 122(1), 129-151.               

2.  Chi, H.; Poorkarimi, H.; Wiener, J;  Shah, S. M. On the exponential growth of solutions to 

non-linear hyperbolic equations. International Journal of Mathematics and Mathematical 

Sciences.1989,12(3),539-545,https://doi.org/10.1155/S0161171289000670.                                                                                                         

3. Minamoto, T. Numerical existence and uniqueness proof for solutions of nonlinear 

hyperbolic equations. Journal of Computational and Applied Mathematics 2001, 135(1), 79-

90.                     

 4.  Krylovas, A.;  Čiegis, R.  A review of numerical asymptotic averaging for weakly 

nonlinear hyperbolic waves. Mathematical Modelling and Analysis 2004, 9(3), 209-222.                           
DOI:10.1080/13926292.2004.9637254.                                                                                                  

5. Ashyralyev, A.; Agirseven, D. Bounded solutions of nonlinear hyperbolic equations with 

time delay.  Electron. J. Differential Equations 2018[Internet], 2018(21), 1-15, Available 

from:  http://ejde.math.txstate.edu . 

6.  Bangerth, W.; Geiger, M.;  Rannacher, R.; Adaptive Galerkin finite element methods for 

the wave equation. Computational Methods in Applied Mathematics 2010, 10(1), 3-48. DOI:  

10.2478/cmam-2010-0001. 



  

128 

 

Ibn Al-Haitham Jour. for Pure & Appl. Sci. 34 (2) 2021 

 

 7.  Al-Haq, A. A. J. M. Numerical Methods for Solving Hyperbolic Type Problems. M.Sc. 

thesis. An-Najah National University, Nablus. 2017.     
8.  Al-Hawasy,J. A.; Jawad, M. A. The Approximation Solution of a Nonlinear Parabolic 

Boundary Value Problem Via Galerkin Finite Elements Method with Crank-Nicolson. 

IHJPAS.  2018, 31(3), 126-134, https://doi.org/10.30526/31.3.2002. 

9.   Al-Hawasy, J. A.;  Jasim, D. K. The Continuous Classical Optimal Control Problems for 

Triple Elliptic Partial Differential Equations. IHJPAS. 2020, 33(1), 143-151, DOI: 

10.30526/33.1.2380.                                                               

10.  Holte, J. M. Discrete Gronwall lemma and applications. In MAA-NCS meeting at the 

University of North Dakota 2009, 24,1-7.                                                                                             

11.   Chen, J.; Jin. Z.; Shi. Q.; Qiu. J.; Liu. W.. Block Algorithm and Its Implementation for 

Cholesky Factorization. ICCGI (2013): The Eighth International Multi-Conference on 

Computing in the Global Information Technology China ICCGI20013.232-236. ISBN: 978-      

1-61208-283-7.