Study with the several resources on Docsity
Earn points by helping other students or get them with a premium plan
Prepare for your exams
Study with the several resources on Docsity
Earn points to download
Earn points by helping other students or get them with a premium plan
These notes will help you understand clearly most of the things needed in Numerical Methods, whether you are an engineering student, Mathematics student, or Statistics student.
Typology: Lecture notes
1 / 45
This page cannot be seen from the preview
Don't miss anything!
'''''' (^) $$$$$
&
$
%
Faculty of Natural Sciences Department of Mathematical Sciences & Computing
Notes for the course
Prepared by
Prof. W. Sinkala
These notes are tailored for the course Numerical Methods (APM37W1). In this course, our focus is on practical applications of Numerical Methods for solving Par- tial Differential Equations (PDEs), rather than delving into the theoretical aspects of Numerical Analysis, which include topics like convergence, stability, and error analysis.
The primary goal of these notes is to elaborate on various techniques used in numerically solving PDEs. By emphasising the practical implementation of these methods, we aim to provide students with a thorough understanding of how they operate on a computer. It is expected that students will actively engage with the material by experimenting with the provided numerical routines using calculators or computers. This hands-on approach is crucial for grasping concepts and routines.
For computational tasks, we will utilise MATLAB, a widely used programming lan- guage known for its simplicity and robust visualisation capabilities. MATLAB’s pop- ularity in academic and industrial settings, coupled with its user-friendly nature and powerful visualisation tools, makes it an ideal choice for our purposes.
During classes, straightforward MATLAB programmes will be provided, and students will be required to comprehend each line’s specific function. This practise aims to foster a deeper understanding of programming logic and prepares students for modifying these programmes to solve specific problems effectively.
Prof. W. Sinkala
February 2024
i
Partial differential equations occur widely in many scientific applications. To give a simple example, consider the temperature in the room you are sitting in. the temper- ature obviously varies from one point to another, that is it is a function of the three cartesian coordinates x, y and z. It also varies with time, t. So the temperature may be written as u = u(x, y, z, t). It can easily be shown using basic laws of physics that u has to satisfy an equation involving partial derivatives with respect to x, y, z and t. Such an equation is called a partial differential equation. In order to solve this equation we require a good deal more information. For instance we should need to about the sources of heat in the room, and the way in which heat is absorbed by the walls and furniture.
Many physical phenomena in such diverse domains as financial mathematics, fluid mechanics, electricity, magnetism, optics, or heat flow are conveniently described by partial differential equations. In general a partial differential equation is an equation that contains partial derivatives. In contrast to ordinary differential equations (ODEs), where the unknown function depends only one variable, in partial differential equations, the unknown function depends on several variables.
Notation. The following (usual) nation is used:
ut =
∂u ∂t
ux =
∂u ∂x
uxx =
∂^2 u ∂x^2
uxt =
∂^2 u ∂x∂t
Partial differential equations are classified according to many things. Classification is an important concept because the general theory and methods of solution usually apply only to a given class of equations. Six basic classifications are:
(a) ut = uxx B^2 − 4 AC = 0 (parabolic) (b) utt = uxx B^2 − 4 AC = 4 > 0 (hyperbolic) (c) uxx + uyy = 0 B^2 − 4 AC = − 4 < 0 (elliptic)
(d) yuxx + uyy = 0 B^2 − 4 AC = − 4 y
elliptic for y > 0 parabolic for y = 0 hyperbolic for y < 0
Remark 1.2.
A general classification diagram is given in Figure 1.
Linearity Linear Nonlinear
Order 1 2 3 4 5 · · · m
Kind of Coefficients (linear equations) Constant Variable
Homogeneity (linear equations) Homogeneous Nonhomogeneous
Number of variables 1 2 3 4 5 · · · n
Basic type (linear equations) Hyperbolic Parabolic Elliptic
Table 1.1: Classification diagram for partial differential equations.
1.3 Boundary and Initial Conditions
Typically, the number of independent solutions of a PDE is infinite. Moreover, the boundaries of the regions of the independent variables over which we desire to solve the PDE are usually not discrete; they are often continuous curves or surfaces. Thus, the complete formulation of a physical system in terms of PDEs requires careful attention not only to the equations that govern the system but also to the correct formulation of the boundary conditions (BCs). For many PDEs we also require (in addition to the BCs) to specify the initial conditions (ICs) of the system in order to obtain a unique solution. For second-order partial differential equations linear boundary conditions can take one of the following three forms:
Example 1.3.
(a) If u(x, t) is the displacement of a vibrating string and its ends are fixed at x = 0 and x = L, then the conditions u(0, t) = 0 and u(L, t) = 0 are Dirichlet conditions.
(b) Suppose u(x, t) is the temperature in a rod of length L. If the rod is perfectly insulated at x = 0 and x = L, the heat flux at these pints is zero. In this case the appropriate boundary conditions (from the Fourier law of heat conduction) are ∂u/∂x(0, t) = 0 and ∂u/∂x(L, t) = 0. These are Neumann are boundary conditions.
(c) Suppose in Example (b) above that poor insulation is used at the end of the rod. The boundary conditions might take the form
u(0, t) +
∂u ∂x
(0, t) = 0 and u(L, t) +
∂u ∂x
(L, t) = u 0
This is an example of Robin’s boundary conditions.
1.4 Problem Set 1
(a) ut = uxx + u (b) ut = uxx + e−t (c) uxx + 3uxy + uyy = sin x (d) utt = uuxxx + e−t
(a) uxx − uxy − 2 uyy = 0 (b) 2uxx + 4uxy + 3uyy + 7u = 0 (c) yuxx − 2 uxy + exuyy − u = 3
Mathematical modelling is the process whereby the evolution or the state of a real-life system is represented by a set of mathematical relations, after proper approximations and idealisations. In this chapter we present a number of partial differential equations which arise in various areas of Applied Mathematics.
Suppose that a string is tightly stretched between two fixed point x = 0 and x = L on the x-axis (Figure 2.1(a)). We assume that the string is perfectly flexible and elastic. At time t = 0 the string is picked up in the middle (Figure 2.1(b)), to a distance h (where h is assumed to be small compared with L. Then the string is released. The
x = 0 (^) x = L
x = L 2
x
y
(a)
(b)
x
y
x = 0 x = L
h
Figure 2.1: Elastic string
problem is to describe the motion which takes place. Suppose that at some instant t,
the string has a shape as shown in Figure 2.2. We shall call Y (x, t) the displacement of point x on the string (measured from the equilibrium position which we take as the x-axis) at the time t.
It turns out the Y (x, t) is governed by the partial differential equation
∂t^2
= a^2
∂x^2
where a^2 = τ /ρ, τ being the tension in the string (assumed constant) and ρ the density (mass per unit length) of the string. Equation (2.1.1) is called the vibrating string equation.
x
y
x
x
x
x
x
y y(x,t) (x+ ,t)
+
x+ x Point
Equilibrium position
Point
Figure 2.2: Vibrating elastic string
Let us now consider the boundary conditions. Since the string is fixed at points x = 0 and x = L, we have
Y (0, t) = 0, Y (L, T ) = 0, for t ≥ 0. (2.1.2)
These conditions state that the displacements at the ends of the string are always zero. Referring to Figure 2.1(b) it is seen that
Y (x, 0) =
2 hx L
, 0 ≤ x ≤
2 h L
(L − x),
≤ x ≤ L
This merely gives the equations of the two straight line portions of that figure. Y (x, 0) denotes the displacements of any point x at t = 0. Since the string is released from rest, its initial velocity everywhere is zero. Denoting by Yt the velocity ∂Y /∂t we may write Yt(x, 0) = 0 (2.1.4)
which says that the velocity at any place x at time t = 0 is zero.
There are many other boundary value problems which can be formulated using the same partial differential equation (2.1.1). For example, the string could be plucked at another point besides the midpoint or even at two or more points. We could also have a string or rope of which one end is fixed while the other is moved up and down according to some law of motion.
It is also possible to generalise the vibrating string equation (2.1.1). For example, suppose that we have a membrane or drumhead in the form of a square in the xy-plane
We call (2.1.9) Laplace’s equation (in 3-dimensions) and ∇^2 the Laplacian after the French mathematician Pierre Simon Laplace (1749–1827), who investigated many of its important and interesting properties.
Our concern with the solution of Laplace’s equation will be restricted to finite domains, and in particular to rectangular domains. The example below is typical of Laplace boundary value problems on rectangular domains. Note that time is not involved in Laplace’s equation.
Example 2.1.1 [Laplace boundary value problem]
uxx + uyy = 0 0 < x < 4 , 0 < y < 5 , uy(x, 0) − 2 u(x, 0) = 2x u(x, 5) = 3 u(0, y) = y^2 ux(4, y) + 3u(4, y) = 0
Suppose that a thin metal bar of length L is placed on the x-axis of an xy coordinate system (see Figure 2.4).
x
y
x = 0 x = L
Figure 2.4: Thin metal bar
Suppose that the bar is immersed in boiling water so that it is at a temperature 100◦C. Then it is removed and the ends x = 0 and x = L are kept in ice so that the ends are at a temperature 0◦C. We shall suppose that no heat escapes from the surface of the bar, ie the surface is insulated. We are interested in knowing what the temperature will be at any place in the bar at any time. Denoting the temperature of the bar by U it is easily seen that U depends on the position x of the bar, as well as the time t (measured from time zero when the bar is at 100◦C) of observation. We denote this dependency by U (x, t). Thus we have the dependent variable U depending on the two independent variables x and t. The mathematical formulation of this problem results in a partial differential equation that must be satisfied by the temperature function, U , of the metal bar, viz.
∂^2 U ∂t
= α^2
∂x^2
where α^2 is called the diffusivity of the material. Equation (2.1.11) is called the (one- dimensional) heat equation or diffusion equation.
Taking the special case where the ends of the rod are are kept at 0◦C and where the initial temperature of the bar is 100◦C, the following boundary conditions result:
U (0, t) = 0, U (L, t) = 0 for t > 0 U (x, 0) = 100, for 0 < x < L
The first two express the fact that the temperature at x = 0 and x = L is zero for all time. The third expresses the fact that the temperature at any place x between 0 and L at time zero is 100◦C. Thus the boundary value problem is the problem of determining the solution of the partial differential equation (2.1.11) which satisfies the conditions (2.1.12).
Remark 2.1.1 The same boundary value problem may have different interpretations, ie, it can arise from apparently different physical situations. For example, suppose that we have an infinite slab of conducting material initially at 100 ◦C where the bounding planes x = 0 and x = L are kept at 0 ◦C. Then the temperature U (x, t) at any location x of the slab at any time t > 0 is given by the boundary value problem defined by partial differential equation (2.1.11) and (2.1.12) above.
It is easy to generalise equation (2.1.11) to the case where heat may flow in more than one direction. For example, for heat conduction in two and in three dimensions we have ∂^2 U ∂t
= κ
∂x^2
∂y^2
and ∂^2 U ∂t
= κ
∂x^2
∂y^2
∂z^2
respectively. We call (2.1.13) the two-dimensional heat conduction equation, and (2.1.14) the three-dimensional heat conduction equation.
NB. The equations (2.1.11), (2.1.13) and (2.1.14) can be written more concisely in terms of the Laplacian operator. For example (2.1.14)can be written as
∂U ∂t
= κ∇^2 U.
In case U does not depend on t we call U the steady-state temperature, and in this case (2.1.14), for example, reduces to
∇^2 U = 0 or
∂x^2
∂y^2
∂z^2
which is Laplace’s equation (compare with (2.1.9) on page 9).
In the study of options in Financial Mathematics, a very important model is the so- called Black-Scholes equation (see, for example, [2] and [10])
∂V ∂t
σ^2 x^2 2
∂x^2
∂x
− rV = 0 (2.1.16)
Finite difference methods provide a powerful approach to solving differential equations and are widely used. Equations with variable coefficients and even nonlinear prob- lems can be treated with these techniques. Generally the error of an approximating solution can be made arbitrarily small. Rounding errors, which inevitably rise during the computational process, can be controlled by a preliminary analysis of the numer- ical stability of finite-difference schemes. In this chapter we introduce the most used finite-difference approximations. We use the Taylor series expansion to derive these approximations.
Let us consider a function u(x, t) depending on two variables x ∈ [0, L] and t ∈ [0, T ]. A discretisation of function u is obtained by considering only the values ui,j on a finite number of points (xi, tj )
ui,j = u(xi, tj ) = u(i∆x, j∆t), i = 1,... , m, j = 0,... , n, (3.1.1)
where ∆x = L/m, ∆t = T /n (see Figure 3.1).
NB. The notation uji is also used in place of ui,j.
A basic role to estimate the error involved in finite-difference approximations of the derivatives of a function is played by the well-known Taylor’s series expansion (see Section ??):
f (x + ∆x) =
X^ n
k=
f k(x)
(∆x)k k!
(∆x)(n+1) (n + 1)!
where ξ lies between x and x + ∆x and f (k)^ denotes the kth derivative of f
f 0 (x) = f (x)
. Noting that the last term is of order (∆x)n+1, (3.1.2) can be written as
f (x + ∆x) =
X^ n
k=
f k(x)
(∆x)k k!
(∆x)n+
(refer to Section ?? for the meaning of O (·))
∆x i∆x^
x
(i∆x, i∆t) i∆t
∆t
t
Figure 3.1: Space-time grid.
3.2 Finite-difference approximation of derivatives
Applying Taylor’s series expansion (3.1.3) to u(xi, tj + ∆t) gives
u(xi, tj + ∆t) = u(xi, tj ) + ut(xi, tj )∆t + O
(∆t)^2
which in view of notation (3.1.1) is written as
ui,j+1 = ui,j + (ut)i,j ∆t + O
(∆t)^2
ie
(ut)i,j =
ui,j+1 − ui,j ∆t
Thus we have the forward difference approximation of ut,
(ut)i,j ≈
ui,j+1 − ui,j ∆t
with a leading error of order ∆t. Similarly for ux we have that
(ux)i,j =
ui+1,j − ui,j ∆x
from which the forward difference approximation formula for ux follows:
(ux)i,j ≈
ui+1,j − ui,j ∆x
and has a leading error of order ∆x.
which imply
(ut)i,j =
4 ui,j+1 − 3 ui,j − ui,j+ 2∆t
(∆t)^2
We thus have a three-point forward approximation formula,
(ut)i,j ≈
4 ui,j+1 − 3 ui,j − ui,j+ 2∆t
for the approximation of ut with a leading error of order O ((∆t)^2 ). Similarly, the three-point forward approximation for ux can be obtained; it is
(ux)i,j ≈
4 ui+1,j − 3 ui,j − ui+2,j 2∆x
with a leading error of order O ((∆x)^2 ).
We start with the forward approximation of utt. Multiplying (3.2.15) by 2 and sub- tracting the result from (3.2.16) leads to
(utt)i,j =
ui,j+2 − 2 ui,j+1 + ui,j (∆t)^2
and hence we have the finite-difference approximation
(utt)i,j ≈
ui,j+2 − 2 ui,j+1 + ui,j (∆t)^2
with a leading error of order ∆t. A similar result holds for uxx:
(uxx)i,j ≈
ui+2,j − 2 ui+1,j + ui,j (∆x)^2
with a leading error of order ∆x. By analogous arguments the backward approxima- tions of the same derivatives can be inferred:
(utt)i,j ≈
ui,j− 2 − 2 ui,j− 1 + ui,j (∆t)^2
(uxx)i,j ≈
ui− 2 ,j − 2 ui− 1 ,j + ui,j (∆x)^2
with leading errors of ∆t and ∆x, respectively.
Summing up (3.2.9) and (3.2.10) and solving for (utt)i,j gives the central-difference approximation formula:
(utt)i,j ≈
ui,j+1 − 2 ui,j + ui,j− 1 (∆t)^2
with leading errors of (∆t)^2. A similar result holds for uxx:
(uxx)i,j ≈
ui+1,j − 2 ui,j + ui− 1 ,j (∆x)^2
with leading errors of (∆x)^2.
Remark 3.2.1 Formulas (3.2.25) and (3.2.26) are more accurate approximations of utt and uxx, respectively, than the forward or backward approximation formulas, (3.2.21), (3.2.22), (3.2.23) and (3.2.24).
We now consider approximations for mixed derivatives uxt. For the forward approxi- mations of uxt, if we apply (3.2.3) to the derivative ux we obtain
(uxt)i,j =
(ux)i,j+1 − (ux)i,j ∆t
which in view of (3.2.3) implies
(uxt)i,j =
ui+1,j+1 − ui,j+1 − ui+1,j + ui,j ∆x∆t
Therefore we have
(uxt)i,j ≈
ui+1,j+1 − ui,j+1 − ui+1,j + ui,j ∆x∆t
with a leading error of order O (∆x) + O (∆t).
Remark 3.2.2 Formula (3.2.29) can also be found by applying (3.2.5) to ut and using (3.2.3). (Verify this!)
Similarly we can derive the back-ward approximation
(uxt)i,j ≈
ui− 1 ,j− 1 − ui,j− 1 − ui− 1 ,j + ui,j ∆x∆t
and the central-difference approximation
(uxt)i,j ≈
ui+1,j+1 − ui− 1 ,j+1 − ui+1,j− 1 + ui− 1 ,j− 1 4∆x∆t
with leading errors of order O (∆x) + O (∆t) and O ((∆x)^2 ) + O ((∆t)^2 ), respectively.
NB. Higher order derivatives can be evaluated by analogous arguments.
3.3 Problem Set 3
(ut)i,j ≈
− 4 ui,j− 1 + 3ui,j + ui,j− 2 2∆t
(ux)i,j ≈
− 4 ui− 1 ,j + 3ui,j + ui− 2 ,j 2∆x
