An equation, which contains the derivatives of one or more dependent variables with respect to one or more independent variables, is called a differential equation (DE).

There are two types of differential equations:

• • Ordinary differential equation (ODE)

• • Partial differential equation (PDE)

Based on the numerical methods formulated and used to solve the ODEs, we further classify the ODE problems into two categories:

• • Initial value problem (IVP)

• • Boundary value problem (BVP)

In Chapters 4 and 5, we will cover the numerical solution of the ODE-IVP and ODE-BVP systems, respectively. The numerical schemes for the PDE systems will be discussed subsequently in Chapter 6.

Key Learning Objectives

• • Learning how to solve systems of non-linear algebraic equations

• • Knowing the bracketing, open, and polynomial-based numerical methods

• • Learning how to develop the numerical algorithm for different systems

• • Simulating a couple of industry-relevant processes with choosing a suitable algorithm

INTRODUCTION

This chapter covers the numerical solution of the systems of non-linear algebraic equations. For this, there are several numerical techniques available, which are often called the iterative convergence methods. Prior to discussing these methods, let us start with a non-linear function of a single variable represented by:

y = f(x) 0 (3.1)

One can solve it graphically by producing a plot of y versus x. As shown in Figure 3.1a, the curve intersects the x-axis at point x*, which is called the root (or zero) of the function f(x), since it satisfies Equation (3.1). It is true that a non-linear algebraic equation may have multiple solutions or roots (e.g., x1*, x2* and x3* ) as shown in Figure 3.1b. In fact, we all know that an nth order polynomial has n solutions.

This graphical method is useful for rough estimates of the root and because of the lack of precision, it has very limited application. However, one can use this rough estimate as an initial guess. Alternatively, one can use the trial-and-error approach, which starts with assuming the initial guess value of x and evaluates whether f(x) = 0. If it is not satisfied (that is almost always true), re-assume another guess for x and again evaluate f(x) to check the improvement in result, if any. This is repeated until f(x) reaches sufficiently close to zero.

It is somewhat obvious that finding root through the graphical or trial-and-error method is insufficient and inadequate for scientific and engineering problems. Thus the necessity of numerical methods arises. Next, we will pay our attention to know various numerical methods and how to solve the non-linear systems by judiciously using them.

Key Learning Objectives

• • Learning some basics of partial differential equation (PDE)

• • Knowing several finite difference methods for PDE

• • Solving parabolic, hyperbolic and elliptic PDEs

• • Developing solution methodology with the use of global index

• • Numerically solving various linear and non-linear process examples

INTRODUCTION

There are many systems which are typically modeled with more than one independent variable [e.g., space (x) and time (t)] to characterize their states. These models thus contain partial derivatives and are constituted with partial differential equations (PDEs). Indeed, the PDEs are used to describe a wide variety of phenomena, such as heat, diffusion, fluid dynamics, gravitation, quantum mechanics, sound, elasticity, among others.

Let us consider a general form of a linear PDE, which has two independent variables, namely x and y, as:

Here, u is the dependent variable. Obviously, this is a second-order PDE.

Rewriting,

Note that here, A, B, C, D, E, F, G and are all functions of x and y, not of u or its derivatives.

BASIC DIFFERENCES BETWEEN ODE AND PDE

A couple of basic differences between the ODE and PDE are highlighted in Table 6.1.

Note that in this chapter, we will confine ourselves mostly to two-dimensional (2D) differential systems, which have two independent variables (i.e., either time and space or two spatial dimensions). However, we will also have a few cases with three independent variables (see, for example, Section 6.10), namely time and two spatial dimensions.

INTIAL AND BOUNDARY CONDITIONS

The partial differential equation systems often include derivative terms on spatial variables. In many cases, time derivative term is also additionally contained. For example, along with a second-order derivative with respect to space variable, the PDE may contain a first-order time derivative as:

For such a case, along with two boundary conditions (BCs), an initial condition (IC) needs to be prescribed. It thus becomes clear that we specify one IC for a first-order time derivative contained in a PDE, two ICs for a second-order time derivative, and so on.

Algebraic equation is an equation in the form of a polynomial containing a finite number of terms and equated to zero (dictionary meaning). Sometimes, the algebraic equations may include transcendental functions (e.g., trigonometric, exponential and logarithmic functions).

We will cover two types of the systems of algebraic equations:

• • Linear algebraic equation system

• • Non-linear algebraic equation system

For solving the linear algebraic equations, various methods are to be discussed in Chapter 2, and in Chapter 3, the numerical methods are covered for solving the non-linear algebraic equations.

This is the introductory part of this text book and it consists of just one chapter, Chapter 1, titled Introduction to Numerical Method and Process Simulation. Here, our goal is to mainly cover some basic concepts related to:

• • Numerical method

• • Process simulation

Broadly speaking, these are the two main topics of this book, along with engineering applications of the numerical algorithms. They are further detailed in subsequent parts of the book, having a total of eight chapters.

Key Learning Objectives

• • Understanding the basics of interpolating polynomial

• • Knowing a couple of commonly used polynomial interpolation approaches

• • Learning how to estimate the errors (truncation errors)

• • Knowing the piecewise polynomial interpolation

INTRODUCTION

In the last chapter, we have learned how to fit an analytical function, f based on the available y versus x data set. The least-squares curve fitting method was used there for this purpose to know the overall trend of our experimental observations. It reveals that the function is not necessarily going through any of the data points and it is expected so, since the data are not precise enough. Now here we would like to cover various techniques of interpolation to find the function that should pass through all of the data points with the presumption that the data are known to be very precise.

For the given (tabulated) data points, (x0, y0), (x1, y1), ….,(xm, ym), the value of y at x is called:

Basically, we need to approximate y = f(x) at non-tabulated abscissas over the interval [x0, xm]for interpolation and outside that interval for extrapolation.

Let us consider a non-linear equation (linear in the parameters, a’s) that we would like to use to predict our data points:

where, for example:

Depending on the given data points, one can choose the polynomial, periodic (sine) or exponential function. However, here our interest is confined to polynomial function.

Now to make the above model Equation (8.1) ready for use, one needs to find the values of all coefficients/parameters, a0, a1…..am. Accordingly, we fix our objective as the construction of the function f(x) as a polynomial of minimum degree. It is fairly true that as m in the polynomial form of Equation (8.1) increases, the difficulty level in finding the associated parameters values gets intensified.

Before leaving the introduction section, let us briefly highlight the basic differences between the least-squares and polynomial interpolation in Table 8.1.

Analytical methods are used to find exact solutions of various physical and engineering problems. But these methods have certain limitations in that they are generally able to solve linear (or linearized) systems or those having low dimensionality and simple geometry. It is a fact that the natural systems (i.e., practical problems) are inherently non-linear and complex. For them, thus numerical methods are inevitable for finding approximate solution. In the pre-computer era, the application of numerical methods was limited as they involved repeated arithmetic operations. Today, with the advent of fast and cheaper digital computers, the role of numerical techniques in solving scientific and engineering problems has exploded.

It is fairly correct that there are several excellent books available on numerical methods in the market. Among them, many of which are specific to only science stream, and a few are relevant to specific engineering discipline. Usually, the theoretical part is covered more in the science streams, whereas application part gets more importance in the engineering disciplines. At this juncture, we felt the necessity of a book for readers who are keen to investigate the wide applicability of various numerical approaches to interesting examples and case studies from both the science and engineering streams. Hence attempt is made here to bridge this gap to some extent.

This text book covers various elegant and advanced numerical techniques that are commonly employed to solve the algebraic and differential equations (ordinary and partial differential equations). To easily understand the theory of numerical methods, they are presented in a simple manner in that the Lemma, Corollary etc. are not included. Furthermore, to gain insights into these methods, all involved calculations in them are shown in details. The issues related to their usability, limitations, and conservativeness (if any) are concisely discussed with beautiful worked-out examples from both science and engineering fields. This apart, there are a couple of advanced and effective numerical techniques covered with their relevant applications.

An appealing part of this book is that a wide range of process units are used to illustrate the numerical methods. At first, the numerical solutions of some mathematically worked-out examples are discussed that are typically solved by the students of all disciplines at the early stage of their undergraduate (UG) course.

You have the hardware and understand its architecture. You have a large problem to solve. You suspect that a parallel program may be helpful. Where do you begin? Before we can answer that question, an understanding of the software infrastructure is required. In this chapter, we will discuss general organization of parallel programs, that is, typical software architecture. Chapter 5 elaborates this further and discusses how to design solutions to different types of problems.

As we have noted, truly sequential processors hardly exist, but they execute sequential programs fully well. Some parts of the sequential program may even be executed in parallel, either directly by the hardware's design, or with the help of a parallelizing compiler. On the other hand, we are likely to achieve severely sub-par performance by relying solely on the hardware and the compiler. With only a little more thought, it is often possible to simply organize a sequential program into multiple components and turn it into a truly parallel program.

This chapter introduces parallel programming models. Parallel programming models characterize the anatomy or structure of parallel programs. This structure is somewhat more complex than that of a sequential program, and one must understand this structure to develop parallel programs. These programming models will also provide the context for the performance analysis methodology discussed in Chapter 3 as well as the parallel design techniques described in Chapter 5.

We will see in Chapter 7 that many efficient sequential algorithms are not so efficient if trivially parallelized. Many problems instead require specially designed parallel algorithms suitable for the underlying system architecture. These parallel algorithms are often designed directly in terms of these programming models.

A program broadly consists of executable parts and memory where data is held, in addition to input and output. A large parallel program usually performs input and output through a parallel file system. We will discuss parallel file systems in Section 5.4, but in the context of the current discussion they behave much like memory – data of some size can be fetched from an address or written to an address by executable parts.

We are now ready to start implementing parallel programs. This requires us to know:

• • How to create and manage fragments (and tasks).

• • How to provide the code for the fragments.

• • How to organize, initialize, and access shared memory.

• • How to cause tasks to communicate.

• • How to synchronize among tasks.

This chapter discusses popular software tools that provide answers to these questions. It offers a broad overview of these tools in order to familiarize the reader with the core concepts employed in tools like these, and their relative strengths. This discussion must be supplemented with detailed documentation and manuals that are available for these tools before one starts to program.

The minimal requirement from a parallel programming platform is that it supports the creation of multiple tasks or threads and allows data communication and synchronization among them. Modern programming languages, Java, Python, and so on, usually have these facilities – either as a part of language constructs or through standard library functions. We start with OpenMP, which is designed for parallel programming on a single computing system with memory shared across threads of a processor. It is supported by many C/C++ and Fortran compilers. We will use the C-style.

OpenMP

Language-based support for parallel programming is popular, especially for single node computing systems. Compiling such programs produces a single executable, which can be loaded into a process for execution, similar to sequential programs. The process then generates multiple threads for parallel execution. OpenMP is a compiler-directive-based shared-memory programming model, which allows sequential programmers to quickly graduate to parallel programming. In fact, an OpenMP program stripped off its directives is nothing but a sequential program. A compiler that does not support the directives could just ignore them. (For some things, OpenMP provides library functions – these are not ignored by the compiler.) Some compilers that support OpenMP pragmas still require a compile time flag to enable that support.

Preliminaries

C/C++ employs #pragma directives to provide instructions to the compiler. OpenMP directives all are prefixed with #pragma omp followed by the name of the directive and possible further options for the directive as a sequence of clauses, as shown in Listing 6.1.

This chapter introduces some general principles of parallel algorithm design. We will consider a few case studies to illustrate broad approaches to parallel algorithms. As already discussed in Chapter 5, the underlying goal for these algorithms is to pose the solution into parcels of relatively independent computation, with occasional interaction. In order to abstract the details of synchronization, we will assume the parallel RAM (PRAM) or the bulk-synchronous parallel (BSP) model to describe and analyze these algorithms. It is a good time for the reminder that going from, say, a PRAM algorithm to one that is efficient on a particular architecture requires refinement and careful design for a particular platform. This is particularly true when “constant time” concurrent read and write operations are assumed. Concurrent reads and writes are particularly inefficient for distributed-memory platforms, and are inefficient for shared-memory platforms as well. It requires synchronization of the processors’ views of the shared memory, which can be expensive.

Recall that PRAM models focus mainly on the computational aspect of algorithm, whereas practical algorithms also require close attention to memory, communication, and synchronization overheads. PRAM algorithms may not always be practical, but they are easier to design than those for more general models. In reality, PRAM algorithms are only the first step toward more practical algorithms, particularly on distributed-memory systems.

Parallel algorithm design often seeks to maximize parallelism and minimize the time complexity. Even if the number of actually available processors is limited, higher parallelism translates to higher scalability in practice. Nonetheless, the work-time scheduling principle (Section 3.5) indicates that low work complexity is paramount for fast execution in practice. In general, if the best sequential complexity of solving the given problem is, say To(n), we would like the parallel work complexity to be O(To(n)). It is a common algorithm design pattern to assume up to To(n) processors and then try to minimize the time complexity. With maximal parallelism, the target time complexity using To(n) processors is O(1). This is not always achievable, and there is often a trade-off between time and work complexity. We then try to reduce the work complexity to O(To(n)), without significantly increasing the time complexity.

Parallel programming is challenging. There are many parts interacting in a complex manner: algorithm-imposed dependency, scheduling on multiple execution units, synchronization, data communication capacity, network topology, memory bandwidth limit, cache performance in the presence of multiple independent threads accessing memory, program scalability, heterogeneity of hardware, and so on. It is useful to understand each of these aspects separately. We discuss general parallel design principles in this chapter. These ideas largely apply to both shared-memory style and message-passing style programming, as well as task-centric programs.

At first cut, there are two approaches to start designing parallel applications:

• 1. Given a problem, design and implement a sequential algorithm, and then turn it into a parallel program based on the type of available parallel architecture.

• 2. Start ab initio. Design a parallel algorithm suitable for the underlying architecture and then implement it.

In either case, performance, correctness, reusability, and maintainability are important goals. We will see that for many problems, starting with a sequential algorithm and then dividing it into independent tasks that can execute in parallel leads to a poor parallel algorithm. Instead, another algorithm that is designed to maximize independent parts, may yield better performance. If a good parallel solution cannot be found – and there do exist inherently sequential problems, for which parallel solutions are not sufficiently faster than sequential ones – it may not be a problem worth solving in parallel.

Once a parallel algorithm is designed, it may yet contain parts that are sequential. Further, the parallel parts can also be executed on a sequential machine in an arbitrary sequence. Such “sequentialization” allows the developer to test parts of a parallel program. If a purely sequential version is already available, or can be implemented with only small effort, it can also serve as a starting point for parallel design. The sequential version can be exploited to develop the parallel application incrementally, gradually replacing sequential parts with their parallel versions. The sequential version also provides performance targets for the parallel version and allows debugging by comparing partial results.

Programs need to be correct. Programs also need to be fast. In order to write efficient programs, one surely must know how to evaluate efficiency. One might take recourse to our prior understanding of efficiency in the sequential context and compare observed parallel performance to observed sequential performance. Or, we can define parallel efficiency independent of sequential performance. We may yet draw inspiration from the way efficiency is evaluated in a sequential context. Into that scheme, we would need to incorporate the impact of an increasing number of processors deployed to solve the given problem.

Efficiency has two metrics. The first is in an abstract setting, for example, the asymptotic analysis of the underlying algorithm. The second is concrete – how well does the algorithm's implementation behave in practice on the available hardware and on data sizes of interest. Both are important.

There is no substitute for measuring the performance of the real implementation on real data. On the other hand, developing and testing iteratively on large parallel systems is prohibitively expensive. Most development occurs on a small scale: using only a few processors, p, on small input of size n. The extrapolation of these tests to a much larger scale is deceptively hard, and we often must resort to simplified models and analysis tools.

Asymptotic analysis on simple models is sometimes criticized because it oversimplifies several complex dynamics (like cache behavior, out-of-order execution on multiple execution engines, instruction dependencies, etc.) and conceals constant multipliers. Nonetheless, with large input sizes that are common in parallel applications, asymptotic measures do have value. They can be computed somewhat easily, in a standardized setting and without requiring iterations on large supercomputers. And, concealing constants is a choice to some degree. Useful constants can and should be retained. Nonetheless, the abstract part of our analysis will employ the big-O notation to describe the number of steps an algorithm takes. It is a function of the input size n and the number of processors p.

Asymptotic notation or not, the time t(n, p) to solve a problem in parallel is a function of n and p. For this purpose, we will generally count in p the number of sequential processors – they complete their program instructions in sequence.

This chapter is not designed for a detailed study of computer architecture. Rather, it is a cursory review of concepts that are useful for understanding the performance issues in parallel programs. Readers may well need to refer to a more detailed treatise on architecture to delve deeper into some of the concepts.

There are two distinct facets of parallel architecture: the structure of the processors, that is, the hardware architecture, and the structure of the programs, that is, the software architecture. The hardware architecture has three major components:

• 1. Computation engine: it carries out program instructions.

• 2. Memory system: it provides ways to store values and recall them later.

• 3. Network: it forms the connections among processors and memory.

An understanding of the organization of each architecture and their interaction with each other is important to write efficient parallel programs. This chapter is an introduction to this topic. Some of these hardware architecture details can be hidden from application programs by well-designed programming frameworks and compilers. Nonetheless, a better understanding of these generally leads to more efficient programs. One must similarly understand the components of the program along with the programming environment. In other words, a programmer must ask:

• 1. How do the multiple processing units operate and interact with each other?

• 2. How is the program organized so it can start and control all processing units? How is it split into cooperating parts and how do parts merge? How do parts cooperate with other parts (or programs)?

One way to view the organization of hardware as well as software is as graphs (see Sections 1.6 and 2.3). Vertices in these graphs represent processors or program components, and edges represent network connection or program communication. Often, implementation simplicity, higher performance, and cost-effectiveness can be achieved with restrictions on the structure of these graphs. The hardware and software architectures are, in principle, independent of each other. In practice, however, certain software organizations are more suited to certain hardware organizations. We will discuss these graphs and their relationship starting in section 2.3.

Another way to categorize the hardware organization was proposed by Flynn and is based on the relationship between the instructions different processors execute at a time. This is popularly known as Flynn’s taxonomy.

SISD: Single Instruction, Single Data

A processor executes program instructions, operating on some input to produce some output. An SISD processor is a serial processor.