NUMERICAL SOLUTIONS FOR PARTIAL DIFFERENTIAL EQUATIONS OF CONSERVATION LAWS

 

The world we live in is in nature a very complicated one. Quit often we can not explain what is happening or perhaps how it happened. But by the help of mathematics we can predict the future occurrences, one is because mathematics is cheaper since instead of running experiments everyday we use mathematical equations to predict the future. Secondly this problems are so complex in nature such that when we relate them to partial difference equations (PDEs) we more often than not have non-linearities or complications coming from the pdes themselves or the from the boundaries that is when we have complicated boundaries [1]. For that reason its is evident we can not solve them analytically and that’s why we need numerical solutions, that means we try to get approximate solutions as accurate as possible.

In Figure (1), the first image shows the Solai tragedy that occurred in Kenya. The second one shows a traffic situation in an intersection. The third one simulates a bullet in motion, while the last shows a linear accelerator used the to treat cancer.

Numerical methods are procedures or techniques that are used to approximate mathematical processes (an example of a mathematical process is an integral). They are schemes that are applied to find an approximation when analytical answers are either absent or impossible to find. Numerical methods for all partial differential equations (PDEs) is made up of two parts: First a space discretisation that transforms the system of PDEs into a system of ODEs and the second is a time discretisation that transforms the system of ODEs into a series of algebraic equations that can be solved using techniques of linear algebra [2].

When modeling the physical world, conservation becomes a fundamental principle to be applied which stretches from quantum mechanics, continuum mechanics to gravitational physics, having prominent examples that include the famous Navier’s equations of elasticity, Euler equations of gas dynamics, Einstein’s equations of gravitation and Maxwell’s equations of electromagnetics.

It therefore does not come as a surprise why there is deep research in the area of mathematical analysis of conservation laws whose stretches back to the Leibniz. The process that involve developing accurate and efficient computational methods used to solve such problems have taken the center stage of scientific pioneers including von Neumann, Lax, and Lions. However, even if a problem may appear very simple, a detail investigation reveals serious hurdles when we attempt to understand the most basic conservation law together with the nature of the solutions.

When we find an exact solution to a problem and in a shorter time than forever, then we say that the problem is "analytically solvable”. For instance, “Victor has 2 mangoes while James has 3 mangoes, how many mangoes do they have together?” is solvable analytically. It’s 5. Exactly 5.

Precisely solving such problems is what we often think that mathematicians are busy doing, however more often than not scientists stumble on problems that can’t be solved analytically. In real sense “mathematics” widely involves finding an answer rather than the answer. While it could be hard to find the exact answer, we can try to find answers which are “arbitrary precise”. This is like saying, you can determine an approximate solution and the more computer time you are willing to invest, the closer that approximate answer will be to the correct answer. This trick is called the “numerical method”. Numerical methods are rather strange: they determine solutions close to the exact answer without ever knowing what that answer is in the first place. So that there is actually an answer: we require numerical methods because majority of problems are not solvable analytically and we know that these methods work because each one of them method comes with a proof that it works.

The choice of a particular numerical method to use in solving a certain problem, depends on the ease in applying it to that particular problem, the ability of the method to produce accurate results when compared to other numerical methods and the consistency of the numerical methods.

Systems of hyperbolic conservation consists of nonlinear and time dependent systems of pdes with simple structure. In one-dimensional space the equations are of the form

Here u is a vector of m-dimension and represent quantities which are conserved state variables for example momentum, energy dynamics and mass problem. Even more precisely, u j represents the state variables density function at position j. The explanation is that at time t the integral ∫ x2 u j(x, t) dx equals to all the quantity in the interval (x1, x2) of state variable.

When we say that we are conserving the state variables we are trying to say that the integral ∫ ∞ u j(x, t)dx is supposed to be fixed respecting t. These functions u j themselves at time t which stand for how the state variables are distributed changes with time. We

are assuming that for equation (1) when we know u(x, t) at a certain points and time enables determination of the flux, of every state variable (x, t).

The functions of the flux are typically functions of u, which are nonlinear that leads to systems (nonlinear) of PDEs. Generally it’s impossible to determine the actual solutions to the equations, therefore the necessity to formulate and discuss numerical schemes for getting those solutions which are approximate.

As a example of equations we are looking at in our entire investigations, we look at a conservation law, typically written as

 

together wit initial condition

subject to some periodic boundary conditions. Suppose that the partition of (0, 1) is given  with

to find the element’s center we use

if (2) and (15) are multiplied by an arbitrary function say v(x) then integrating by parts over Ij, we obtain the weak problem statement, that is

               (4)

 

                               (5)

where  represents the limit from the right and the le respectively.

The problem statement is that we find a solution uh to u for every t ∈ (0, T ), with v(x) and uh(x, t) belonging to semi-finite dimension space

    (6)

This means that, in the space dimension  and v are polynomials of degree K. Since uh and v are both discontinuous at  , then the ambiguity in the terms of (4) which involves the non-linear fuxes  should be substituted by numerical fuxes which entirely depends on  i.e

 

         (7)

 

which is yet to be chosen. Therefore if we use Legendre’s polynomials Pm to expand the approximation uh in terms of basis functions

             (8)

Then the test functions vh are considered to be equal to the functions of the basis, that is, ,[3], whereas at the same time we invoke the orthogonality property of Legendre’s polynomials

 

                               (9)

 

where

 

 

 

 

Notice that we apply the property

There are multiple reasons as to why we consider this particular group of equations alone:

One of the earliest development in an attempt to solve conservation laws was the Dis- continuous Galerkin (DG) method. DG methods were originally developed around 1970s [5] as a way of solving partial differential equations numerically [4] . In 1973 Reed and Hill [6] discovered a DG method to solve the NEUTRON TRANSPORT EQUATION (NTE) equation.

where σ, a(x) ∈ R and u is to be found

This equation has roots going back more than a century ago to the Boltzamann’s equation [7].

which was formulated initially to facilitate the study of kinetic theory of gases. Further study of radiation transport in atmosphere lead to several analytical Solutions the trans- port problems in early 1930s. The physics surrounding these problems, however, caused a confinement in interest to one dimensional semi-infinite medium geometries. During the advent of nuclear chain reactors in 1940s there became interest involving neutral particle transport problems in a wider range of configurations in geometry which were found in the applications of nuclear reactor and radiation shielding. A number of good analytical methods to solve the transport problems who pursued since the 1940s. The Weiner-Hopf method [8], singular eigenfunction expansions among other analytical methods provided a great deal of input into the general nature of the transport processes through the study of greatly idealized configurations, which included the Milne problem.

At the same time there was an increase in sophisticated numerical methods which began being developed, concurrent with the rapidly increasing use of digital computers as a computational power. Among these numerical methods the discontinuous galerkin was conceived.

Until their most recent development these method made it’s way into the heart of computational fluid dynamics where they found use in a wide variety of applications including but has not been limited to nonlinear conservation laws, the compressible Navierstokes equation and Hamilton Jacobi-like equations.

This quest, however, is far from trivial simply because of two main reasons

Thus, when constructing numerical methods for these problems, it must be guaranteed that the discontinuities of the approximate solution are those which are physically relevant. In addition, it must be ensured that the appearance of a discontinuity in the approximated solution does not bring about sparious oscillations that tamper with the quality of the approximation; on the other hand, when ensuring this, the method must remain accurate enough near that discontinuity in order to capture the possibly rich structure of the exact solution.

Owing to their finite element nature, the DG methods have the following main advantages over other classical finite volume and finite difference methods:

The original finite element method was introduced in 1973 by Reed and Hill for solving the neutron transport equation (10).The usefulness of the method was recognised by LeSaint and Raviart who in 1974 published its first mathematical analysis. To display the method, we multiply the equation by a test function v and integrate over an arbitrary subset of Ω say, K. A er a formal integration by parts, we get

σ (u, v)K − (u, a.∇v)K + (a.nku, v)dK = ( f , v)K

where nk denotes the outward unit normal of dK, and

Next, we construct a triangulation τh = K of Ω, and take our approximate solution uh to be a polynomial of degree at most k on each elements of the triangulation. The approximate solution uh is then determined as the unique solution of the following weak formulation:

∀k ∈ τh :

σ (uh, v)K − (uh, a.∇v)K) + (h, v)dK = ( f , v)K, ∀v ∈ Pk(K),

where Pk(K) denotes the space of polynomials of degree at most k on the element Ka and h on the numerical flux given by

Note that the value lim(x − sa) is nothing but the value of uh upstream the characteristic direction a. As a consequence the degree of freedom of the approximate solution h the element K can be computed in terms of the values of uh upstream the characteristics hitting dK. In other words, the approximate solution uh can be computed element by element when the elements are suitably ordered according to the characteristic direction a.

 

The first analysis of the DG method as applied to ordinary differential equations, was performed in 1974 by LeSaint and Raviart [9] who showed that the method is strongly stable of order 2k + 1 at mesh points, and that the Gauss-Radau discretization of the DG method is also of order 2k + 1 when polynomials (piecewise) of degree k are applied. It is rather interesting that about a year later before the introduction of the DG method by Reed and Hill, another mathematician, Hulme had studied a method for ordinary differential equations which word is same week formulation as the DG method but employed a continuous approximate solution uh; this method is, however, of order k at mesh points.

A study of global error control for ordinary differential equations for DG method was done in 1994 by Estep and French. Another work on DG methods for ODEs was carried out in 1981 by Delfour, Hager and Trochu; they introduced a class of DG methods which are proven to give an order of accuracy up to 2k + 2 at the mesh points. Recently, Schotzau and Schwab have obtained a new estimate on the size of time step needed to solve the implicit system of equations determined by the DG method by way of a simple fixed point iteration method.

In 1988,Johnson [10] gave an analysis of error control for the DG method for harder ODEs. And lastly, in 1996, Bottcher and Rannacher introduced a new global error control method for ODEs by using the DG method.

With a variety of well-tested and successfully applied methods, one cannot help but ask why there is any need to consider another method. To embrace this, let us take off by trying to understand the strengths and weaknesses of the standard methods. we will consider how one-dimensional conservation law for the solution u(x, t)

corresponding to a set of boundary and initial condition on, ∂ Ω. Where f = flux and h(x, t) is a function which has been prescribed. To come up with any numerical method that solves a pde calls for consideration of two factors

This draws the line between separate methods giving different properties.

Consider a very simple method used over the years and perhaps for the oldest time, called the FDM. Here, xk, is given in space and using difference methods we approximate derivatives; i.e,

                     (12)

 

where ug is solution and fg is flux numerically approximated, whereas

represents size of grid. To construct a FDM we need, within the surrounding of every point xk, the flux and its solution are therefore considered as approximated numerically using polynomials.

whereas ai(t) and bi(t) are calculated ensuring the approximated function exists at points, xk. putting this approximation in equation (11), gives the residual

Clearly,

 

 

In the table below we give the summary of methods. By a look at it one should remember that this comparison only look at some basic problems and that many of those problems which have been addressed and restrictions can be rectified and overcome in a variety of ways. In addition, the comparison gives an insight of which shortcoming one should struggle to resolve when trying to come up with a new method.

Table 1. Generic properties of methods used in discretizing partial differential equations

On the table a C shows success, while × shows that there is a shortcoming in that method. finally, (C) shows that the method, when subjected to modifications, can be used to solve search problems but it should not the recommended for as a natural choice.

We will briefly look into the standard Galerkin finite element method for the approximate solution of initial-boundary value problem for the heat equation (note that the method is equally applicable to similar equations),

the Laplacian. Before we proceed with the discussion of this problem we need to revise on some basic material around the finite element method for its corresponding stationary problem, the Dirichlet problem for Poisson’s equation,

−∆u = f

in Ω, with u = 0, on dΩ. using variational formulation of this problem, we shall define an approximation of the solution u.as a function uh which belongs to a finite dimensional linear space Sh of functions of x with certain properties. This function, is the simplest case a continuous, piecewise linear function on some partition of Ω, will be a solution of finance system of linear algebra equation. We show basic error estimates for this approximate solution in energy at least square norms. Looking back at the parabolic problem which we first write in a weak form, we continue to discretize this problem, first in the special variable x, which gate is result in an approximate solution uh(., t) in the finite element space Sh, for t ≥ 0, as a solution of an initial value problem for a finite- dimensional system of ordinary differential equations. We then define the fully discrete approximation by application of some finite difference time stepping method to this finite dimensional initial value problem. This yields an approximate solution U = Uh which belongs to Sh at discrete-time levels. Errors estimates will be derived for both the special and fully discrete solutions.

The outline of the thesis is as follows:

Chapter 2: Having research on the earlier numerical methods for PDEs, this chapter will focus on the strides made towards achieving the ultimate goal. The discontinuous galerkin method will be considered in particular where we will look at how the method has been applied in attempts to solve nonlinear equations like the Parabolic, Maxwell’s and Visco-elastic flows in the recent past. We will briefly present other original numerical methods and describe their theoretical and computational developments in the frame- work of linear hyperbolic systems.

In addition we will review the application of the Finite Difference method in various equations including but not limited to the Viscid Shallow Water equations and the Isothermal Euler equations.

 

Chapter 3: In this chapter we will discuss the finite difference method in details, it’s development will be considered. We will derive the and shallow water equations in non- conservative and conservative form and then using finite difference method we will solve it.

Chapter 4: For better understanding of how the finite difference method works we will run simulations of the results in chapter 3 using the MATLAB so ware before concluding our study in chapter 5.

 

No Review Found.