Fast Convergence Methods for Hyperbolic Systems of Balance Laws with Riemann Conditions

In this paper, we develop an accurate technique via the use of the Adomian decomposition method (ADM) to solve analytically a 2 × 2 systems of partial differential equation that represent balance laws of hyperbolic-elliptic type. We prove that the sequence of iteration obtained by ADM converges strongly to the exact solution by establishing a construction of fixed points. For comparison purposes, we also use the Sinc function methodology to establish a new procedure to solve numerically the same system. It is shown that approximation by Sinc function converges to the exact solution exponentially, also handles changes in type. A numerical example is presented to demonstrate the theoretical results. It is noted that the two methods show the symmetry in the approximate solution. The results obtained by both methods reveal that they are reliable and convenient for solving balance laws where the initial conditions are of the Riemann type.


Introduction
Many mathematical and practical problems in physics can be modeled by balance laws, which are systems of nonlinear partial differential equations, usually of hyperbolic type. Balance laws describe very large branches of scientific fields, specifically in fluid and gas dynamics, quantum mechanics, and astrophysics.
In this paper, we present the mathematical formulation for finding numerical solution using Sinc basis functions for 2 × 2 systems of balance laws, while for comparison purposes, we also use the Adomian series technique to find approximate analytical solutions for the same model system that which can be written as The solution U, which is a vector of two components u(x, t) that might represent the length of a wave, and another component, say v(x, t) which is the velocity of that wave. i.e., U(x, t) = [u(x, t), v(x, t)] T . On the other hand the flux function can be represented as F(U) = [ f (u, v), g(u, v)] T , and finally, the source term has the form G(U) = [h 1 (x, t), h 2 (x, t)] T . where F is a smooth function on IR 2 of U, and the derivative of F possesses strictly nonlinear or linearly degenerate eigenfunctions and eigenvalues. The system (1) is of mixed type, it is of elliptic type, for all values of u, v that satisfy ε = {(u, v) ∈ R 2 : ( f u + g v ) 2 < ( f u + g v )} and of hyperbolic type for data that lying in H = {(u, v) ∈ R 2 : ( f u + g v ) 2 > ( f u + g v )}. We consider in Equation (1) the simplest piece-wise constant initial data, i.e., for certain dispositions of the Riemann initial data For the sinc method, the basis functions on the interval (−∞, ∞) for z ∈ D E , and for the time interval (0, T) are derived from the composite translated sinc functions and, π((z−kh)/h) , z = kh 1, z = kh, k = 0, ∓1, ∓2, . . .
Let f be a function defined on IR, then for h > 0 we define the Whittaker cardinal series The properties of (6) were studied and surveyed in depth in Stenger and in Lund [15][16][17]. They are based on the infinite strip D d defined in C as To approximate by the Sinc method, we will refer to an important class of functions, see [17] called Play-Wiener class, denoted by W(π/h) which is the family of all analytic functions that are in C, and satisfy some decaying conditions. For a function f and its derivative f to be in the class W(π/h), the kth derivative of the function f can be approximated by i−n denotes the kth derivative of the cardinal function. Now we present a fundamental definition that has an important role in the development of Sinc methods on arc Γ in C. For the simply connected domain D, let φ be a map that take D onto D d with φ(a) = −∞ and φ(b) = ∞. If the inverse map of φ is denoted by ψ, define When we evaluate a function and its derivative with respect to φ, we obtain Since we will replace the first derivative by it matrix approximations, here we define the m × m matrix I (1) whose in − th entry is given by δ (1) i−n . An important class of functions denoted by L α (D) play crucial rule in approximating derivatives will be needed in next formula (see, [17]). Let F ∈ L α (D), α is a positive constant, then for h = πd (αN) we have An application to the above formula, and in order to approximate the first derivative with respect to x for the function u(x, t) in the domain (−∞, ∞) leads to the approximation where With m x = 2N + 1 and the skew-symmetric matrix I (1) m x is defined by (9). Finally, we shall indicate a general formula for approximating the integral ν a F(t)dt, ν ∈ Γ. If F ∈ W(π/h), then for all k ∈ Z, To simplify the situation, we define the new matrix I (−1) , as, for N positive integer, define the (2N + 1) × (2N + 1) matrix I (−1) by

The Sinc-Galerkin Method: Balance Laws
Among the benefits of the Sinc method, is that the approximate solution is approaching the exact solution exponentially, as well as the Sinc method, it carefully handles systems in which the state changes from elliptic to hyperbolic and vis versa, as this depends on the change of the eigen-values of the system. To find an approximate solution for the balance law appeared in (1), without loss of generality, we begin with a 1-dimensional balance law defined in the domain R × (0, T 0 ) where the initial condition is taken to be Suppose k(x) ∈ L α (D), and for the step-size h = πd αN . Integrate Equation (12) with respect to time, and collocation with respect to the variable x will produce a system of Volterra integral equations of the form where the column vector u(t) represent u(t) = (u −N x (t), . . . , u N x (t)) T , here we use the notation u i (t) = u i (x i , t), k = (k(z −N x ), . . . , k(z N x )) T and H = (H(z −N x ), . . . , H(z N x )) T , while the square matrix A = −1 m x represents the discrete Sinc approximation for the first derivative. Subsequently an application of the use of the indefinite integration formula with respect to time variable t, and since our t ∈ (0, T 0 , then the domain . . , N t , to be the basis function for the interval (0, T 0 ). Now, we define a matrix E = h t I as defined in (11), with m t = 2N t + 1. Then the solution of Equation (12) is in matrix form is given by the rectangular m x × m t matrix U = [u ij ]: we use • to represents the Hadamard multiplication. For convergence purposes for the obtained solution we state two Theorems, in which the proof can be done by just resemble the theorems in [16].

Theorem 1 ([16]).
Let the matrix U be as defined in (15). Then for the number of points in both space and time N x , N t to be bigger than some constant that depends on d, α, we have To guarantee that our approximate solution as obtained in the discrete system (15) converges uniquely to the true solution, we state the following theorem.
Theorem 2 ( [16,18]). For any positive constant R, we can find another positive constant T 0 so that whenever U 1 − U 0 < R 2 , then the solution of (15) is unique. Also, the iteration system leads to a solution that converges to the unique solution.
Next, we state and prove a fact regarding the behavior of the solution, whenever the initial condition satisfies a decay condition. Consider the initial value problem for a balance law of the form with initial condition u(x, 0) = u 0 (x), x ∈ IR. We need to assume some decaying conditions on the right hand side of Equation (16).

Definition 1.
We say a function u(x, t) decays exponentially in x and uniformly for t ∈ (0, T) if |u(x, t)| ≤ C exp(−α|x|) for some positive constants al pha, C independent of t.
, and can be written as u 0 (x) = sech (αx)g(x), where the function g(x) is analytic in the strip D d and grows no more than a polynomial in the strip D d , and suppose the source term h(x, t, u) and its derivatives with respect to x, y, u are all exponentially decaying in x and uniformly for small t. Then the solution u(x, t) of (16) is in the class L α (D).
Proof. In order to find the solution of (16), we consider the system of first order ODEs.
such that The system (17) with condition (18) implies that However, h(x, t, u) and its derivatives with respect to the variables x, and t are all exponentially decaying in x and uniformly for small t. This implies that is Lipschitz in the two variables u, x. In addition, so, by the Picard theorem, the solution u(x, t) is unique.

Treatment of Non-Zero Boundary Conditions
The obtained solution in (15) was valid if the initial condition satisfy some certain decaying conditions, i.e, the initial condition has to belong to the class L α . In Section 5, we will discuss problems which have non-zero boundary conditions with the Riemann-type initial condition, which means that the initial condition is no longer in L α . To illustrate the situation, we consider the problem in (1) with the initial condition (3), and since the Sinc functions composed with the various conformal maps, S(m, h x ) • φ, are zero at the end-points of the domain of the problem, and since the boundary conditions are non-homogeneous Dirichlet conditions, then the change of variables (similarly for v(x, t)), will convert system (1) to another one with homogeneous Dirichlet conditions, with smooth initial condition. We substitute the transformation (19) into Equation (1) that leads to a new balance law for the new unknownũ(x, t) with an initial condition that is in the class L α for |x| large. Therefore, to find the solution forũ(x, t) we can proceed as mentioned above.

The Adomian Decomposition Method (ADM)
A lot of published research papers talked about the method of Adomian, and presented it in a detailed way. Here, for the sake of non-repetition, we will suffice to mention the references that contributed to the development of the method. In addition, because of its abundance and where we cannot limit it, we mention the following [19][20][21][22]. We start by rewriting the system in (1) as follows: (20) or, we may express Equations (20) using linear differential operators as In the above, the operators L t , L x are defined to be L t = ∂ ∂t and L x = ∂ ∂x . Performing the inverse operator L −1 t = t 0 . . . dt PDE system in (21) yields where The ADM [19,20] suggests that our approximate solution for unknown functions u(x, t) and v(x, t) are given by the infinite series while the four nonlinear operators φ i (u, v), i = 1, 2, 3, 4 can be expressed as an infinite series of Adomian polynomials as φ ( u, v) = ∑ ∞ n=0 A ni , where A ni are what is called Adomian's polynomials that can be computed according to specific formula set in [20]. Here, the general formulas for Adomian's polynomials A ni , i = 1, 2, 3, 4, has the form These formulas are smoothly found using coding programs to get as many terms as we want. Plug the nonlinear system (23) into Equation (22) and we arrive at the recursive relation The construction of the solutions of u(x, t) and v(x, t) are given by Regarding the convergence of the above two series, there are previous studies that dealt with the topic in detail and in general, we mention [21].

Convergence of the ADM Approximation
Here we give a proof to show that our obtained solution using ADM converges to the exacts solution via the use of the fixed-point theory. The partial differential equation that we are encountered to solve is which can be written as where H is the Hilbert space, and Comparing terms, we end up with the algorithm u 0 = F(x, t), u n+1 = A n (u 0 , u 1 , . . . , u n ). Then the approximate solution is given by a finite sum of the series S n = ∑ n n=1 u n = u 1 + u 2 + · · · + u n . If S 0 = 0, S n+1 = N(S n + u 0 ) = ∑ n k=0 A k . If the limit exist, then S = lim n→S n is the solution of the fixed-point functional equation S = N(u 0 + S) ∈ H. Need to show that the sequence {S n } is Cauchy. The distance between two consecutive iterations S n and S n+1 is a ball S r of radius r is given by To prove that {S n } is Cauchy sequence, then for all positive integers m, n where m > n, we have Choose r ∈ (0, 1) such that, S m − S n can be made arbitrarily small provided n is large enough for all m, n ∈ IN. With this choice of r and u 0 , u 1 , . . . we see that all iterations will remain in the ball There is an integer N such that u m − u n < , ∀n > N, for any m. So the sequence is Cauchy in the Hilbert space. Hence converges to some S ∈ H., i.e., lim n→∞ S n = S, where S = ∑ ∞ n=0 u n . Solving Equation (26) is the same as solving the functional equation N(S + u 0 ) = S, by assuming that N is a continuous operator, we get When viewing numerical calculations in the last section, and for purposes of comparison we will calculate what is called the order of convergence, we will provide the following definition Definition 2. A sequence S n is said to converge to S of order p ∈ IN, if there is a positive constant α such that Let us examine the order of convergence for the Cauchy sequence S n . We write Taylor's series for N(S n + u 0 ) about (S + u 0 ) as However, it is known that N(S + u 0 ) = S and N(S n + u 0 ) = S n+1 , we may rewrite the above equation as Assume our operator N ∈ C p [a, b] with the property that N (m) (S + u 0 ) = 0, for all m = 1, 2, . . . , p − 1, while it satisfy N p (S + u 0 ) = 0, therefore we have divide by (S n − S) p and take limit as n → ∞, we obtain Using the fact that S n is Cauchy that converges to S, the above equation reduces to Therefore, the order of convergence of S n is p.
For the purposes of calculating the order of convergence, we assume that the first four iterations of the computed solution are given by x k−2 , x k−1 , x k and x k+1 so the following formula will be used to calculate the order of convergence As mentioned in [23], the first study comparing ADM and Picard's method was in 1987 by Rach [24]. The study was done by reviewing several examples. Rach showed that solving differential equations of linear type using ADM was equivalent to the classical method of successive approximations, which is known as Picard iteration. Also it was shown that the Picard method gives more accurate solutions, but it needs more time than ADM; this happens because the main advantage of Adomian's method over the Picard method is the ease of computation of successive terms.

Applications: Riemann Type
To investigate the accuracy of the ADM compared with the Sinc method, we choose examples with known solutions that allows for a more complete error analysis. So a 2 × 2 system of balance laws with Riemann type conditions is tested numerically by using the Sinc function methodology, and for comparison purposes, we also solve the same model by ADM. The example reported here is selected to show the convergence of the two schemes. Consider the system where T is some small constant. In addition, with flux function F; IR 2 → IR 2 , such that with the approximations to Riemann type condition given by The left and right boundary conditions are .
In Section 3 when we construct the approximate solution by means of the Sinc methodology, it was necessary that the initial condition belonged to the class of functions L α (D) previously mentioned in Section 3, but we know that the Riemann condition does not belong to this family, so a transformation had to be made to make the primitive condition belong to the family. Since the boundary conditions are non-homogeneous, then the transformatioñ will convert the given boundary conditions to homogeneous conditions, provided that |x| → ∞. Now, after substituting the transformation (29) into Equation (28), we get a new system with the unknowñ U(x, t) = (ũ(x, t),ṽ(x, t)) T given by where W(x, 0) = (w 1 (x, 0), w 2 (x, 0)) T . Equation (30) can be solved forŨ(x, t), and then, using the transformation in (29), for U(x, t). To proceed, since the space and the time domains are (−∞, ∞) and (0, T), respectively, choose the conformal maps φ(x) = x and Υ(t) = log( t T−t ). We solve the above system using ADM subject to the initial conditions u(x, 0), v(x, 0). Integrate the system in (28) with respect to t and using the initial conditions u(x, 0), v(x, 0), we arrive at and, v( In the above two integrals, set u(x, t) = ∑ ∞ n=0 u n (x, t), v(x, t) = ∑ ∞ n=0 v n (x, t), and for the nonlinear terms, we find Adomian polynomials as The first few terms of Adomian polynomials are given by Substitute the above infinite sum for Adomian polynomials, together with the approximate solutions into Equations (32) and (33), exactly as mentioned in Section 5. Balancing terms in Equations we obtain an approximate solution of the form Table 1 shows that the method converges for T = 2 using Sinc methodology, where the second column reports the supremum norm of the error between the exact solution and the Sinc approximate u−solution, while the third column reports the error between the exact solution and the Sinc approximate v−solution. The error in the approximate u-solution using ADM is reported in Figure 1. Figure 2 reported the v(x)-solution using ADM for different values of time t. The error in approximating u, v show that most of the error concentrates at origin, which caused by approximating Riemann type condition by a smooth function. Through our study, and upon consulting the Tables 1 and 2, we are able to find the values for the order of convergence using Equation (27), it was found that the order of convergence for the ADM is q = 1.0007, while it is approximately q = 2.843 using the Sinc method. This is in accordance with the theories formulated in Section 4, which show that the Sinc methodology is faster, better, and more comprehensive in finding an approximate solution. Because of the nature of the issues under discussion, it is noted that the two methods show the symmetry of the approximate solution. A quick look at in calculated approximate solutions for both u(x, t) and v(x, t) shows the symmetry of the solutions where the error value at the point (x, t) is the same as at (−x, t).
The main goal of this research paper is to compare the two methods together, but if the reader wants more ways to solve systems from the same context, we cite the following [22]. It was shown in [25] that traditional methods, like finite difference method is unstable when the system is of mixed-type when Riemann conditions are involved, which means that the two methods used in this research are an improvement over the previous traditional methods. To present a physical and engineering application for the purposes of explaining the importance of the topic, please review the research paper [18] where the p−system was solved. Finally, what we intend to refer to in a future research in the hope of discussing it from all sides is the importance of the Telegraph system. As an illustration, we shall consider conductors (such as telephone wires or submarine cables) in which the current I = I(x, t) may leak to ground. The resultant decrease in current is governed by I x = −GV − CV t , where V = V(x, t) is voltage, G = G(x) is the conductance to ground, and C = C(x) is the capacitance with the ground. The change in voltage is governed by V x = −RI − LI t , where R = R(x) is the resistance and L = L(x) is the inductance of the cable.