Numerical Methods Book of Notes
Uploaded by
Shyam MahendraNumerical Methods Book of Notes
Uploaded by
Shyam MahendralOMoARcPSD|51243766
ENGR20005 Book of Notes
Numerical Methods for Engineering (University of Melbourne)
Scan to open on Studocu
Studocu is not sponsored or endorsed by any college or university
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
Numerical Methods in Engineering (ENGR20005)
Book
A. Ooi
a.ooi@unimelb.edu.au
July 24, 2020
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
Contents
1 Mathematical Preliminaries 5
2 Root Finding 11
2.1 Finding roots of equations . . . . . . . . . . . . . . . . . . . . . . . . 12
2.1.1 Graphical Method . . . . . . . . . . . . . . . . . . . . . . . . 13
2.2 Bracketing methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
2.2.1 The Bisection Method . . . . . . . . . . . . . . . . . . . . . . 15
2.2.2 Method of False Position . . . . . . . . . . . . . . . . . . . . . 17
2.3 Open methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
2.3.1 Fixed (One) Point Iteration . . . . . . . . . . . . . . . . . . . 21
2.3.2 Newton Raphson Method . . . . . . . . . . . . . . . . . . . . 24
2.3.3 Secant Method . . . . . . . . . . . . . . . . . . . . . . . . . . 29
2.4 System of nonlinear equations . . . . . . . . . . . . . . . . . . . . . . 30
3 Systems of linear algebraic equations 35
3.1 Direct Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
3.1.1 Gauss Elimination . . . . . . . . . . . . . . . . . . . . . . . . 37
3.1.2 LU Decomposition . . . . . . . . . . . . . . . . . . . . . . . . 42
3.1.3 Using LU decomposition to find the inverse of a matrix . . . . 51
3.2 Iterative Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
3.2.1 Point Jacobi . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
3.2.2 Gauss-Seidel . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
3.2.3 Convergence of iterative scheme . . . . . . . . . . . . . . . . . 59
3.3 Newton’s method for nonlinear equations (revisited) . . . . . . . . . . 62
4 Least Squares 67
4.1 Linear least squares approximation . . . . . . . . . . . . . . . . . . . 68
4.2 Polynomial least squares approximation . . . . . . . . . . . . . . . . . 69
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
4 CONTENTS
5 Interpolation 73
5.1 Newton Interpolation . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
5.1.1 Linear Newton Interpolation . . . . . . . . . . . . . . . . . . . 74
5.1.2 Quadratic Newton Interpolation . . . . . . . . . . . . . . . . . 75
5.1.3 General form of Newton’s Interpolating Polynomials . . . . . . 77
5.2 Lagrange Interpolating Polynomials . . . . . . . . . . . . . . . . . . . 79
5.3 Spline Interpolation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
5.3.1 Linear splines . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
5.3.2 Quadratic spline . . . . . . . . . . . . . . . . . . . . . . . . . 89
5.3.3 Cubic Splines . . . . . . . . . . . . . . . . . . . . . . . . . . . 94
6 Differentiation and Integration 105
6.1 Approximation of the 2nd Derivative . . . . . . . . . . . . . . . . . . 107
6.2 Integration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109
6.3 Simpson’s rule . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113
7 Boundary Value Problem 119
8 Ordinary Differential Equations 123
8.1 Euler’s method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123
8.2 Backward Euler Method . . . . . . . . . . . . . . . . . . . . . . . . . 126
8.3 Trapezoidal or Crank Nicolson method . . . . . . . . . . . . . . . . . 127
8.4 Runge Kutta methods . . . . . . . . . . . . . . . . . . . . . . . . . . 129
8.4.1 Second Order Runge Kutta Method . . . . . . . . . . . . . . . 130
8.4.2 4th Order Runge Kutta Scheme (RK-4) . . . . . . . . . . . . 133
8.5 Stability and error analysis . . . . . . . . . . . . . . . . . . . . . . . . 135
8.6 Systems of Ordinary Differential Equations . . . . . . . . . . . . . . . 138
8.6.1 Multi-dimensional Runge-Kutta methods . . . . . . . . . . . . 139
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
Chapter 1
Mathematical Preliminaries
In this chapter, some of the theorems that will be used for the rest of the course will
be reviewed. Only a brief descrition is given here. More details can be found in the
excellent text by [1].
Theorem 1 (Rolle’s Theorem) Suppose f is a continuous function in the domain
(a, b). If f (a) = f (b) = 0, then a number c exist in (a, b) such that f ′ (c) = 0. This
theorem is illustrated in Fig. 1.1.
Exercise 1.1: For the function f (x) = sin(x) in the domain (0, π),
(a) Sketch the function f (x) in the domain specified above.
(b) What is the value of f (0) and f (π).
(c) What is the value of c such that f (c) = 0 ? Your answer should confirm that
0 < c < π.
Theorem 2 (Mean Value Theorem) If f is a contiuous function and differen-
tiable between (a, b), then a number c exist such that
f (b) − f (a)
f ′ (c) = . (1.1)
b−a
This theorem is illustrated in Fig. 1.2.
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
6 CHAPTER 1. MATHEMATICAL PRELIMINARIES
f(x)
df(x)/dx=0
b
a c x
Figure 1.1: Figure illustrating Rolle’s theorem
Slope=(f(b)-f(a))/(b-a)
Slope=df(c)/dx
y=f(x)
a c b x
Figure 1.2: Figure illustrating Mean value theorem.
Theorem 3 (Extreme Value Theorem) If f (x) is a continuous function in [a, b],
then there exist two numbers, c1 and c2 in [a, b] such that f (c1 ) ≤ f (x) ≤ f (c2 ).
Furthermore, if f is differentiable on (a, b), then c1 and c2 occur either at endpoints
a and b or where f ′ (x) = 0. This theorem is illustrated in Fig. 1.3.
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
f(c2)
y=f(x)
f(c1)
a c2 c1 b x
Figure 1.3: Figure illustrating extreme value theorem.
Exercise 1.2: Find the maximum value for |f (x)| where
1
f (x) = (2 − ex + 2x) . (1.2)
3
Do this on the intervals [0,1] and [0,0.5].
Theorem 4 (Taylor’s Theorem) If f ∈ C n [a, b], that f (n+1) exist on [a, b], and
x0 ∈ [a, b]. For every x ∈ [a, b] there exists a number ξ(x) between x0 and x with
f (x) = Pn (x) + Rn (x),
where
n
X f (k) (x0 )
Pn (x) = (x − x0 )k
k=0
k!
f ′′ (x0 ) 2 f (n) (x0 )
′
= f (x0 ) + f (x0 )(x − x0 ) + (x − x0 ) + · · · + (x − x0 )n
2! n!
and
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
8 CHAPTER 1. MATHEMATICAL PRELIMINARIES
f (n+1) (ξ(x))
Rn (x) = (x − x0 )n+1
(n + 1)!
Pn (x) is the nth Taylor Polynomial for f about x0 and Rn (x) is the remainder
term which is also known at the truncation error.
Exercise 1.3:
Find the second Taylor polynomial P2 (x) for the function f (x) = ex cos(x)
about x0 = 0.
(a) For 0 ≤ x ≤ 0.5, find upper bound for the error between P2 (x) and f (x)
(|f (x) − P2 (x)|) using the error formula, and compare it to the actual error.
(b) For 0 ≤ x ≤ 1.0, find upper bound for the error between P2 (x) and f (x)
(|f (x) − P2 (x)|) using the error formula. Compare it to the actual error.
R1 R1
(c) Approximate 0 f (x)dx using 0 P2 (x)dx.
R1
(d) Find an upper bound for the error in (c) using 0 |R2 (x)|dx, and compare
it to the actual error.
Part solution to the above exercise is shown in Fig. 1.4.
Theorem 5 (Weighted Mean Value Theorem for Integrals) If f is a contin-
uous function and differentiable between (a, b) and if g(x) does not change sign in
(a, b), then there exist a number c such that
Z b Z b
f (x)g(x)dx = f (c) g(x)dx. (1.3)
a a
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
1.8 P2(x)=1+x
1.6
1.4 f(x)
1.2
0.8
0.6
0.4
0.2
0
0 0.2 0.4 0.6 0.8 1
x
Figure 1.4: Figure showing part solution to Exercise 1.3.
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
10 CHAPTER 1. MATHEMATICAL PRELIMINARIES
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
Chapter 2
Root Finding
The text by [2] has this excellent example where root finding becomes important.
Consider the forces acting on a mass thrown out from a plane. According to Newton’s
law:
dv
nett force on an object = m (2.1)
dt
where m is the mass of the turtle, v velocity of the turtle and t is time.
For a mass in free fall,
nett force = weight of mass + drag on turtle
= mg − cv
Therefore
dv c
=g− v (2.2)
dt m
Exercise 2.1: Solve Eq.(2.2) and show that
gm ct
v(t) = 1 − e− m (2.3)
c
if v = 0 at t = 0. This function is shown in Fig. 2.1.
Suppose that the velocity of the mass is governed by Eq. (2.3), find the value of
the drag coefficient, c, such that the mass of, m = 5, can attain a prescribed velocity,
v = 10, at a particular time, t = 9. Use g = 10. In the context of the problem
11
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
12 CHAPTER 2. ROOT FINDING
gm/c
v(t)
Figure 2.1: Figure showing the function v(t) defined by Eq. (2.3)
mentioned above, m, v, and t are constants (parameters) of the problem. The only
variable that we can adjust is c. To solve this problem, rewrite Eq. (2.3) as
gm ct
f (c) = 1 − e− m − v(t).
c
Substituting the numerical values of m, t, g and v gives
50 9c
f (c) = 1 − e− 5 − 10. (2.4)
c
So the task of solving the problem reduces to merely finding the value of c, such that
f (c) = 0, for the given values of v, m and t. The graph of f (c) (Eq. (2.4)) is shown in
Fig. 2.2. We need to find a value of c such that f (c) = 0. In numerical analysis this
problem is called finding the roots of equation. There are many methods of finding
the the roots of an equation. They are given in the following sections.
2.1 Finding roots of equations
The methods of solving roots of equations can be broadly classified as
• Graphical method
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
2.1. FINDING ROOTS OF EQUATIONS 13
10
2
f(c)
−2
−4
−6
−8
−10
0 1 2 3 4 5 6 7 8 9 10
c
Figure 2.2: Figure showing the function f (c) defined by Eq. (2.4)
• Bracketing methods
• Open methods
All these methods have their advantages and disadvantages.
2.1.1 Graphical Method
One crude way of finding roots for an equation is to just plot the graph and visually
inspect where the function evaluates to zero. For example, plot f (c) and find the
value of c where f (c) = 0. From Fig. 2.2, the root at f (c) occur when c is approx-
imately 5. The exact solution is c = 4.99382. As you can probably imagine, this
method is not very practical because the accuracy of the solution is very poor. You
are just relying on how accurate your eyes can read the graph ! Furthermore, it is
very inefficient if you ever need to solve f (c) = 0 for different values of m, g and t.
However, the graphical method have a few advantages that must be stated. By
drawing a graph of the function, you can see exactly how the function behaves. It
should also be clear how many roots the equation can have. The graphical method
is also often used to obtain the initial guess for the root of the equation for more
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
14 CHAPTER 2. ROOT FINDING
accurate methods such as the Bracketing and Open methods (see below). There are
various Bracketing and Open methods available and they are shown in Fig. 2.3.
Numerical methods for finding
roots of equations
Bracketing methods Open methods
Bisection Method of One point Secant Newton-
method false postion iteration method Raphson
method
Figure 2.3: Classification of root finding methods
2.2 Bracketing methods
There are two different Bracketing methods. They are
• Bisection Method
• Method of false position
These methods are known as bracketing methods because they rely on having two
initial guesses, xl , which is the lower bound and xu , which is the upper bound. xl
and xu must bracket (be either side of) the root. This requirement can give rise to
various situations which is illustrated in Fig. 2.4.
• If f (xu ) and f (xl ) are of different sign, then there must be at least one root,
xr , in between xu and xl . i.e. xl < xr < xu . This is shown in Fig. 2.4 (a).
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
2.2. BRACKETING METHODS 15
• If f (xu ) and f (xl ) are of different sign, there might be an odd number of roots
(see Fig. 2.4 (b)).
• Figure 2.4 (c) shows that if f (xu ) and f (xl ) are of the same sign, then there
might not be any root in between xu and xl .
There are, however, several exceptions to the rules above. When the function is
tangential to the x-axis, multiple roots occur (see Fig. 2.5 (a)). As is shown in Fig.
2.5 (b), in regions in the vicinity of discontinuities, there might also be excetions to
the rules above.
2.2.1 The Bisection Method
To find roots of an equation using the Bisection method, follow the steps below
1. Obtain two guesses for the root, xl and xu , that “brackets” the real root, xr .
These two guesses can be found by plotting a graph of f (x) and see where it
approximately crosses the x axis. In general, if
f (xu )f (xl ) < 0
then you can be reasonably sure that you have at least one root in between xl
and xu .
2. Estimate the root of f (x) to be
′ xl + xu
xr =
2
3. Make a choice:
′ ′ ′
• If f (xr ) = 0, then xr is the root i.e. xr = xr . Stop. There is nothing more
to do.
′
• If f (xl )f (xr ) < 0, then the root is in the lower subinterval. Therefore set
xu = xr and go back to STEP 2.
′
• If f (xu )f (xr ) < 0, then the root is in the upper subinterval. Therefore set
xl = xr and go back to STEP 2
You usually cannot get exactly f (xr ) = 0. In practice, if |f (xl )f (xr )| < ǫ, where ǫ is
a small value (the stopping criterion), you would take xr to be the root.
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
16 CHAPTER 2. ROOT FINDING
f(x) (a) f(x) (b)
xr
xu xu
xl xl xr xr
xr x x
f(x) (c) f(x) (d)
xr xr
xl xu xl xu
x x
Figure 2.4: Illustration of various situations that could occur when one brackets the
roots.
Exercise 2.2: Use the Bisection method to obtain the root of the following
equation
50 9c
f (c) = 1 − e− 5 − 10.
c
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
2.2. BRACKETING METHODS 17
f(x) (a) f(x) (b)
xr
xu xl
xl xr x xr xu x
Multipleroots
occurhere
Figure 2.5: Illustration of the exception to the rules of the bracketing methods.
The matlab file for solving this exercise is given in Matlab Program 2.1. You
might like to check your answer with the output from this Matlab Program.
2.2.2 Method of False Position
The method of false position is based on the following observation. If f (xl ) is closer
to zero than f (xu ), then the root is likely to be closer to xl than xu . See Fig. 2.6.
Using similar triangles, we can obtain the following expression
−f (xl ) f (xu )
= (2.5)
xr − xl xu − xr
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
18 CHAPTER 2. ROOT FINDING
Matlab Program 2.1 A matlab program illustrating the Bisection method.
%%%%%%%%%%%%%%%%%%%%%%%%%%
% This file can be found in CM/MATLAB/RootFinding/Bisection.m
% Written by Andrew Ooi 3/6/2003
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%clearing all variables
clear all;
EPS=1.0e-6;% EPS, small value which will determine the accuracy of your answer
maxiter=50; %maximum number of iterations
c_l=3.0; %Initial lower guess of root
c_u=9.0; %Initial upper guess of root
for i=1:maxiter,
c_r=(c_l+c_u)/2.; %Guess that the real root is in between c_l and c_u
f_l=(50/c_l)*(1-exp(-(9*c_l/5)))-10;
f_u=(50/c_u)*(1-exp(-(9*c_u/5)))-10;
f_r=(50/c_r)*(1-exp(-(9*c_r/5)))-10;
% Storing results in a matrix
Results(i,1)=i;Results(i,2)=c_r;Results(i,3)=f_r;
% testing to see if we have found the root
% If the root is found then exit the for loop
if abs(f_r) < EPS
break;
end
%Algorithm to find out whether the upper or lower guess
% is closer to the actual root
if (f_r*f_l)<0
c_u=c_r;
elseif (f_r*f_u)<0
c_l=c_r;
else
break;
end;
end;
% Displaying results
disp(’ iter c_r f(c_r)’);
disp(Results);
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
2.2. BRACKETING METHODS 19
f(x)
f(xu)
xl xr
xu x
f(xl)
Figure 2.6: Figure illustrating the method of false position.
Exercise 2.3: Show that Eq. (2.5) reduces to
xu f (xl ) − xl f (xu )
xr =
f (xl ) − f (xu )
or
f (xu )(xl − xu )
xr = xu − (2.6)
f (xl ) − f (xu )
So, to use the method of False Position to find roots, use Steps 1-3 of the Bisection
method but replace Step 2 with the formula given in Eq. (2.6).
Exercise 2.4: Use the method of false position to obtain the root of
50 9c
f (c) =1 − e− 5 − 10
c
You can check your answer by modifying Matlab Program 2.1 to use the
method of false position to find the roots of f (c). Compare with your hand
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
20 CHAPTER 2. ROOT FINDING
written solution.
Figure 2.7 shows the solution of solving Eq. (2.4) using both Bisection and
method of false position with initial lower and upper guess cl = 3 and cu = 9
respective. It can be seen that there is a monotonic convergence to the exact solution
using the method of False position. On the other hand the Bisection method oscillates
towards the exact answer. It is also clear that the method of False position converges
slightly faster to the exact answer.
9
Exact solution
Bisection method
Method of false position
4
1 2 3 4 5 6 7 8 9 10
iter
Figure 2.7: Figure illustrating the differences between the Bisection method and the
method of false position.
2.3 Open methods
The open methods methods which require either only one initial guess or two initial
guesses which do not have to bracket the actual root. It follows that it is easier to
obtain the initial guess for open methods (you can choose almost anything !). In
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
2.3. OPEN METHODS 21
general, if they open methods converges, they have a tendency to do so much more
rapidly than bracketing methods. However, unlike bracketing methods which always
converge (i.e. the root for a particular equation is found), the open methods of root
finding MAY NOT converge. The open methods that will be discussed in this course
are
• Fixed-point iteration
• Newton-Raphson method
• Secant method
2.3.1 Fixed (One) Point Iteration
To find the root of
f (x) = 0, (2.7)
arrange the function so that x is on the left hand side of the equation. i.e.
x = g(x).
We will rewrite the above equation as
xi+1 = g(xi ). (2.8)
Eq. (2.8) above provides us with a iterative formula to find the roots of Eq. (2.7).
Note that the functional form of Eq. (2.7) is not changed. Eq. (2.8) is just a
rearrangement of Eq. (2.7).
Suppose we want to find the root to
f (x) = e−x − x = 0. (2.9)
We will rearrange it to get
x = e−x
I.e. g(x) = e−x and h(x) = x. Hence, the formula for one point iteration scheme for
solving Eq. (2.9) is
xi+1 = e−xi
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
22 CHAPTER 2. ROOT FINDING
Exercise 2.5: Use one point iteration to obtain the root of
50 9c
f (c) = 1 − e− 5 − 10
c
The convergence of fixed point iteration
The iterative equation for fixed-point iteration is
xi+1 = g(xi ) (2.10)
Since f (xr ) = 0 and Eq. (2.10) is obtained by just rearraging f (x), xr must satisfy
xr = g(xr ). (2.11)
Subtracting Eq. (2.11) from (2.10) gives
xi+1 − xr = g(xi ) − g(xr ).
The above equation can be re-written as
g(xi ) − g(xr )
xi+1 − xr = (xi − xr ) (2.12)
(xi − xr )
The derivative mean value theorem states that
dg g(b) − g(a)
(x = ξ) =
dx b−a
where ξ is somewhere in between a and b. In plain English, derivative mean value
theorem say that there is a point (x = ξ) of a function (g) in between a and b such the
derivative of that function is given by g(b) − g(a)/(b − a). If we apply the derivative
mean value theorem to the bracketed term in Eq. (2.12), we get
dg g(xr ) − g(xi )
(ξ) = (2.13)
dx xr − xi
where ξ is somewhere in between xr and xi . Substituting Eq. (2.13) into Eq. (2.12),
we get
dg
xi+1 − xr = (ξ) (xi − xr )
dx
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
2.3. OPEN METHODS 23
The above equation says that the newer guess will be closer to the actual root than the
older guess providing that |dg(ξ)/dx| < 1 . So this is one condition for convergence.
Note that it does not mean that if you have |dg(ξ)/dx| > 1, the solution process will
blow up. Use this condition only as a good guide. Sometimes you can still get a
converged solution if you start have |dg(ξ)/dx| > 1. It really depends on the problem
you are trying to solve.
Exercise 2.6: Find the root for the following equation
f (x) = ex − 10x
by using the fixed point iteration. Try to find the solution with g(x) = ex /10
and x0 = 1.0 and x0 = 4.0 as your initial guesses. Which one of your solutions
blow up ? Why ? Can you get a solution with x0 = 3.5 ?
Exercise 2.7: Let
f (x) = x2 − 2x − 3 = 0
(a) Find the two roots for f (x) analytically.
(b) Show that f (x) can be rearranged into the following forms
√
x = g1 (x) = 2x + 3
3
x = g2 (x) = x−2
(x2 −3)
x = g3 (x) = 2
(c) With these three expressions for gi (x), use the fixed point iteration method
to find the approximate root for f (x). Start with the initial guess of x0 = 4.
Compare your solutions to the answer you get in part a).
(d) Sketch all the gi (x) and show what is happening graphically.
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
24 CHAPTER 2. ROOT FINDING
2.3.2 Newton Raphson Method
The Newton Raphson method is the most widely used root finding method. It
converges quickly and only need one initial guess.
From fig. 2.8, it is easy to express xi in terms of f (xi ), f ′ (xi ) and xi+1 .
f(x)
Slopedf(xi)/dx
xi+1 xi x
Figure 2.8: Figure illustrating the Newton Raphson method
f (xi )
f ′ (xi ) =
xi − xi+1
f (xi )
xi − xi+1 = ′
f (xi )
which gives the following iterative formula for the Newton Raphson method
f (xi )
xi+1 = xi − (2.14)
f ′ (xi )
So to find a root using Newton-Raphson method, do the following:
1. Let the initial guess be xi
2. Find xi+1 by using Eq. (2.14)
3. Let xi = xi+1 repeat steps 2 and 3 until you feel your answer is accurate enough.
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
2.3. OPEN METHODS 25
Exercise 2.8: Use Newton Raphson method to find the root of
50 9c
f (c) =1 − e− 5 − 10
c
You should be able to check your answer by modifying Matlab program 2.1 or
2.2.
Exercise 2.9:
Use Newton Raphson method to find the root of
f (x) = cos x − x
The matlab file showing how one can use Newton Raphson method for solving
Ex. 2.9 is given in Matlab Program 2.2. You can check your hand written
solution with the output from this matlab program.
Exercise 2.10: Let
f (x) = x2 − 2x − 3 = 0
(a) Find the two roots for f (x) analytically.
(b) Find the roots of f (x) using the Newton-Raphson method. Start with x0 =4
and x0 = −2
(c) Sketch f (x) and show the convergence graphically.
One of the reasons why Newton Raphsons method is so popular is that it con-
verges very fast. Proof:-
The Taylor series expansion of a function f (x) can be written as
f ′′ (ξ)
f (xi+1 ) = f (xi ) + f ′ (xi )(xi+1 − xi ) + (xi+1 − xi )2
2!
where ξ lies somewhere in the interval from xi to xi+1 . If we let xi+1 = xr , then we
will have f (xi+1 ) = f (xr ) = 0. Put this information into the above equation to get
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
26 CHAPTER 2. ROOT FINDING
Matlab Program 2.2 A matlab program illustrating the Newton Raphson method.
%%%%%%%%%%%%%%%%%%%%%%%%%%
% This file can be found in CM/MATLAB/RootFinding/NewtonRaphson.m
% Written by Andrew Ooi 3/6/2003
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;%clearing all variables
EPS=1.0e-8;% EPS, small value which determine the accuracy of your answer
maxiter=50; %maximum number of iterations
xi=0.2; %Initial lower guess of root
fi=cos(xi)-xi;
for i=1:maxiter,
% Storing results in a matrix
Results(i,1)=i;Results(i,2)=xi;Results(i,3)=fi;
fi=cos(xi)-xi;
dfi=-sin(xi)-1;
xi=xi-fi/dfi;
% testing to see if we have found the root
% If the root is found then exit the for loop
if abs(fi) < EPS
break;
end
end;
% Displaying results
disp(’ iter x_i f(x_i)’);
disp(Results);
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
2.3. OPEN METHODS 27
f ′′ (ξ)
0 = f (xi ) + f ′ (xi )(xr − xi ) + (xr − xi )2 (2.15)
2!
We can also rearrange Eq. (2.14) to get
0 = f (xi ) + f ′ (xi ) (xi+1 − xi ) . (2.16)
Subtracting Eq. (2.16) from Eq. (2.15) will give
f ′′ (ξ)
0 = f ′ (xi )(xr − xi+1 ) + (xr − xi )2 .
2!
One can then rearrange the above equation to get
f ′′ (ξ)
xr − xi+1 = − (xr − xi )2 (2.17)
2f ′ (xi )
thus
f ′′ (ξ)
|xr − xi+1 | = (xr − xi )2 (2.18)
2f ′ (xi )
According to Eq. (2.18), the magnitude of the error in Newton Raphson method is
roughly proportional to the square of the previous error. Therefore we have quadratic
convergence. Note that quadratic convergence of the Newton Raphson method can
only be achieved if f ′ (xr ) 6= 0. If f ′ (xr ) = 0, then another method has to be
introduced in order to achieve quadratic convergence.
Exercise 2.11:
The formula for Newton-Raphson method is given by
f (xi )
xi+1 = xi − (2.19)
f ′ (xi )
which can be used to find the root, xr , of a function, f (x). It has been been
shown in lectures that the Newton-Raphson method has quadratic convergence
providing that f ′ (xr ) 6= 0. In cases where f ′ (xr ) = 0 one will not be able to
achieve quadratic convergence with Eq. (2.19). In order to achieve quadratic
convergence for cases where f ′ (xr ) = 0, we have to define another function, µ(x),
where
f (x)
µ(x) = .
f ′ (x)
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
28 CHAPTER 2. ROOT FINDING
(a) Assume that f (x) could be written as
f (x) = (x − xr )2 q(x), (2.20)
where q(xr ) 6= 0
(i) Show that f (xr ) = f ′ (xr ) = 0
(ii) Show that µ(xr ) = 0 but µ′ (xr ) 6= 0.
(iii) The results above show that finding the root of f (x) is equivalent
to finding the root of µ(x). Hence, show that a modified version of
the Newton-Raphson method that can be used in order to achieve
quadratic convergence is
f (xi )f ′ (xi )
xi+1 = xi − . (2.21)
[f ′ (xi )]2 − f (xi )f ′′ (xi )
(iv) Use Eqs. (2.19) and (2.21) to obtain the root of
f (x) = e(x−1) − x. (2.22)
You might like to plot the function to see how it looks like. Start with
an initial guess of x0 = 5. Which of the two methods give you faster
convergence ? The exact answer is xr = 1.
(b) Solve
x3 + 4x − 9 = 0 (2.23)
using both Eq. (2.19) and Eq. (2.21). Start with an initial guess of x0 = 5.
(c) Compare the convergence and discuss the advantages and disadvantages of
each schemes.
Several pitfalls of Newton Raphson method
• Requires the evaluation of the derivative of a function. Sometimes, this is
difficult or inconvenient to do.
• Have lots of problems in the vicinity of local maxima or minima because
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
2.3. OPEN METHODS 29
f ′ (xi) = 0 if xi is at a local maxima or minima. This is illustrated in Fig.
2.9.
f(x)
xi xi+1 x
Figure 2.9: Figure illustrating the pitfalls of Newton Raphson method.
2.3.3 Secant Method
One of the main reasons in introducing the secant method is that in the Newton
Raphson method, you have to analytically calculate the derivative f ′ (xi ). Sometimes,
this is very inconvenient. In the secant method, the derivative of a function is
approximated by a backward finite divided difference
f (xi ) − f (xi−1 )
f ′ (xi ) ≈
xi − xi−1
Put the above equation into Eq. (2.14) gives the following formula for the Secant
method
(xi − xi−1 ) f (xi )
xi+1 = xi − (2.24)
(f (xi ) − f (xi−1 ))
Notice that Eq. (2.24 requires two initial estimates of x, i.e. xi and xi−1 . However,
because the function f(x) is not require to change signs between the estimates, it is
not classified as a bracketing method.
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
30 CHAPTER 2. ROOT FINDING
False position Secant
f(x) f(x)
xl xi-1
xr ' xu x xr ' xi x
f(x) f(x)
xr'
xl xr'
xu x xi xi-1 x
Figure 2.10: Figure illustrating the difference between the Secant method (which is
classified as an open method) and the method of false position (which is classified as
bracketing method).
2.4 System of nonlinear equations
So far, we have only talked about methods of finding roots for one equation. What if
we need to simultaneously find roots for two or more equations ? For example, what
if you to find a set of x and y values so that they satisfy the following two equations.
y = x2 + 1 (2.25)
y = 3 cos(x) (2.26)
These two equations are illustrated in Fig. 2.11.
In general, what we would like to do is to find all the xi ’s so that we satisfy all
of the following set of equations.
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
2.4. SYSTEM OF NONLINEAR EQUATIONS 31
4 y=x2+1
2
Solution to the system of equations
0 y=3cos(x)
−1
−2
0 0.5 1 1.5 2 2.5 3
x
Figure 2.11: Figure showing Eqs. (2.25) and (2.26). The point of intersection is the
solution we are after.
f1 (x1 , x2 , ......., xn ) = 0
f2 (x1 , x2 , ......., xn ) = 0
f3 (x1 , x2 , ......., xn ) = 0
.
.
.
fn (x1 , x2 , ......., xn ) = 0
For the purpose of this course, let’s just confine ourselves to n = 2. So we wish
to find x1 = x and x2 = y so that the following two equations are satisfied
u(x, y) = 0
v(x, y) = 0
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
32 CHAPTER 2. ROOT FINDING
The multivariable Taylor series for u and v can be written as
∂u ∂u
ui+1 = ui + (xi+1 − xi ) (xi , yi ) + (yi+1 − yi ) (xi , yi )
∂x ∂y
∂v ∂v
vi+1 = vi + (xi+1 − xi ) (xi , yi ) + (yi+1 − yi ) (xi , yi )
∂x ∂y
For simplicity, the following notation will be used
∂u ∂u
(xi , yi ) ≡
∂x ∂x i
So the above two equations can be re-written as
∂u ∂u
ui+1 = ui + (xi+1 − xi ) + (yi+1 − yi )
∂x i ∂y i
∂v ∂v
vi+1 = vi + (xi+1 − xi ) + (yi+1 − yi )
∂x i ∂y i
which can be simplified to be
∂u ∂u ∂u ∂u
xi+1 + yi+1 = −ui + xi + yi (2.27)
∂x i ∂y i ∂x i ∂y i
∂v ∂v ∂v ∂v
xi+1 + yi+1 = −vi + xi + yi (2.28)
∂x i ∂y i ∂x i ∂y i
if we set ui+1 = vi+1 = 0. The above two equations are a set of two linear equations.
They can easily be solved.
Exercise 2.12: Solve Eqs. (2.27) and (2.28) for xi+1 and yi+1 and show that
vi ∂u
∂y
− u i ∂y
∂v
i i
xi+1 = xi + ∂v
∂u
∂x i ∂y
− ∂u ∂y
∂v
∂x i
i i
∂v
∂u
ui ∂x i − vi ∂x i
yi+1 = yi + ∂v
∂u ∂u ∂v
∂x i ∂y
− ∂y ∂x i
i i
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
2.4. SYSTEM OF NONLINEAR EQUATIONS 33
Exercise 2.13: Provide a simultaneous solution to the two equations out-
lined below using the two equation Newton-Raphson method
y = x2 + 1
y = 3 cos(x)
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
34 CHAPTER 2. ROOT FINDING
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
Chapter 3
Systems of linear algebraic
equations
In many problems of engineering interest, you will need to be able to solve a set of
linear algebraic equations. The methods can be broadly classified as shown in Fig.
3.1.
Solution of linear
algebraic equations
Direct methods Iterative methods
Gauss-Elimination Gauss-Seidel
Gauss-Jordan Gauss-Jacobi
LU Decomposition
Figure 3.1: Classification of different methods for solving a system of algebraic equa-
tions
Consider the following set of equations
35
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
36 CHAPTER 3. SYSTEMS OF LINEAR ALGEBRAIC EQUATIONS
a11 x1 + a12 x2 + a13 x3 + ...... + a1n xn = c1
a21 x1 + a22 x2 + a23 x3 + ...... + a2n xn = c2
.
.
.
.
an1 x1 + an2 x2 + an3 x3 + ...... + ann xn = cn
where all the a’s and c’s are constants. We need to find all the xi ’s such that all
the above equations are satisfied. Sometimes, the matrix notation will be used
[A] {X} = {C} (3.1)
where
a11
a12 . . a1n
a21
a22 . . a2n
[A] =
a31
a32 . . a3n ,
.
. . . .
an1
an2 . . ann
x1
x2
{X} = x3
.
xn
and
c1
c2
{C} = c3 .
.
cn
Equation (3.5) can be solved by just inverting the matrix [A] and then multiply
[C] by [A]−1 .
{X} = [A]−1 {C}
But this might be difficult and time consuming, sometimes impossible to do if A is
large. Can we numerically try to solve it ?
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
3.1. DIRECT METHODS 37
3.1 Direct Methods
In these methods, no iterations are involved.
3.1.1 Gauss Elimination
If you want to use the Gauss elimination procedure to solve the following system of
linear equations,
a11 x1 + a12 x2 + a13 x3 + ...... + a1n xn = c1 = a1,n+1
a21 x1 + a22 x2 + a23 x3 + ...... + a2n xn = c2 = a2,n+1
.
.
.
.
an1 x1 + an2 x2 + an3 x3 + ...... + ann xn = cn = an,n+1
you have to:
• first augment [A] with {C}. We have to next do forward elimination to get the
matrix into a form that is easily solvable.
(1) (1) (1) (1) (1)
a11 a12 a13 . . . . a1n | a1,n+1
(1) (1) (1) (1) (1)
a21 a22 a23 . . . . a2n | a2,n+1
(1) (1) (1) (1) (1)
a31 a32 a33 . . . . a3n | a3,n+1
. . . . . . . . | .
. . . . . . . . | .
. . . . . . . . | .
. . . . . . . . | .
. . . . . . . . | .
. . . . . . . . | .
(1) (1) (1) (1) (1)
an1 an2 an3 . . . . ann | an,n+1
• Eliminate x1 from the second and subsequent equations. For example, to elim-
(1) (1)
inate x1 from the ith row, multiply row 1 by ai1 /a11 and subtract from row
i.
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
38 CHAPTER 3. SYSTEMS OF LINEAR ALGEBRAIC EQUATIONS
(1) (1) (1) (1) (1)
a11 a12 a13 . . . . a1n | a1,n+1
(1) (1) (1) (1) (1)
a21 a22 a23 . . . . a2n | a2,n+1
(1) (1) (1) (1) (1)
a31 a32 a33 . . . . a3n | a3,n+1
. . . . . . . . | .
. . . . . . . . | .
. . . . . . . . | .
. . . . . . . . | .
. . . . . . . . | .
. . . . . . . . | .
(1) (1) (1) (1) (1)
an1 an2 an3 . . . . ann | an,n+1
↓
(2) (2) (2) (2) (2)
a11 a12 a13 . . . . a1n | a1,n+1
(2) (2) (2) (2)
0 a22 a23 . . . . a2n | a2,n+1
(2) (2) (2) (2)
0 a32 a33 . . . . a3n | a3,n+1
. . . . . . . . | .
. . . . . . . . | .
. . . . . . . . | .
. . . . . . . . | .
. . . . . . . . | .
. . . . . . . . | .
(2) (2) (2) (2)
0 an2 an3 . . . . ann | an,n+1
where
(2) (1)
aij = aij for i= 1 and j= 1, 2, 3.., ,n+1
(2)
aij = 0 for i=2, 3,...,n and j= 1
(1)
(2) (1) a (1)
aij = aij − i1 (1) a1j for i= 2, 3,...,n and j= 2, 3,....n+1
a11
• Next, eliminate x2 from the third and subsequent equations. For example, to
(2) (2)
eliminate x2 from the ith row, multiply row 2 by ai2 /a22 and subtract from
row i. Hence,
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
3.1. DIRECT METHODS 39
(2) (2) (2) (2) (2)
a11 a12 a13 . . . . a1n | a1,n+1
(2) (2) (2) (2)
0 a22 a23 . . . . a2n | a2,n+1
(2) (2) (2) (2)
0 a32 a33 . . . . a3n | a3,n+1
. . . . . . . . | .
. . . . . . . . | .
. . . . . . . . | .
. . . . . . . . | .
. . . . . . . . | .
. . . . . . . . | .
(2) (2) (2) (2)
0 an2 an3 . . . . ann | an,n+1
↓
(3) (3) (3) (3) (3)
a11 a12 a13 . . . . a1n | a1,n+1
(3) (3) (3) (3)
0 a22 a23 . . . . a2n | a2,n+1
(3) (3) (3)
0 0 a33 . . . . a3n | a3,n+1
. . . . . . . . | .
. . . . . . . . | .
. . . . . . . . | .
. . . . . . . . | .
. . . . . . . . | .
. . . . . . . . | .
(2) (2) (2)
0 0 an3 . . . . ann | an,n+1
where
(3) (2)
aij = aij
(3)
aij = 0
!
(2)
(3) (2) ai2 (2)
aij = aij − (2)
a2j
a22
• If you keep repeating the steps above, you will end up with a matrix that look
something like
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
40 CHAPTER 3. SYSTEMS OF LINEAR ALGEBRAIC EQUATIONS
(n) (n) (n) (n) (n) (n)
a11 a12 a13 a14 . . . a1n | a1,n+1
(n) (n) (n) (n) (n)
0 a22 a23 a24 . . . a2n | a2,n+1
(n) (n) (n) (n)
0 0 a33 a34 . . . a3n | a3,n+1
(n) (n) (n)
0 0 0 a44 . . . a4n | a4,n+1
. . . . . . . . | .
(3.2)
. . . . . . . . | .
. . . . . . . . | .
. . . . . . . . | .
. . . . . . . . | .
(n) (n)
0 0 0 . . . . ann | an,n+1
where
(k) (k−1)
aij = aij
(k)
aij = 0
(k−1)
!
(k) (k−1) ai,k−1 (k−1)
aij = aij − (k−1)
ak−1,j
ak−1,k−1
We can solve the system of equations given by Eq. (3.2) by using back substitution.
First look at the last row of Eq. (3.2)
(n)
a(n)
nn xn = an,n+1
(n)
an,n+1
xn = (n)
(3.3)
ann
So we can get a numerical value for xn . Next, look at the 2nd last row of Eq. (3.2)
(n) (n) (n)
an−1,n−1 xn−1 + an−1,n xn = an−1,n+1
(n) (n)
an−1,n+1 − an−1,n−1 xn
xn−1 = (n)
(3.4)
an−1,n−1
Since we already know xn from Eq. (3.3), we can also get a numerical value for xn−1
from Eq. (3.4). Similarly, looking at the 3rd last row of Eq. (3.2)
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
3.1. DIRECT METHODS 41
(n) (n) (n) (n)
an−2,n−2 xn−2 + an−2,n−1 xn−1 + an−2,n xn = an−2,n+1
(n) (n) (n)
an−2,n+1 − an−2,n xn + an−2,n−1 xn−1
xn−2 = (n)
an−2,n−2
In general, the formula for getting the numerical value for xi is
(n)
an,n+1
xn = (n)
ann
and
n
P
(n) (n)
ai,n+1 − aij xj
j=i+1
xi = (n)
aii
for i = n − 1, n − 2, n − 3, .., 2, 1.
The pitfalls of the Gauss elimination method for solving a linear system of equa-
tions are:
• It is very expensive, compared to other methods that are currently available.
• The formula for the forward elimination stage requires that we divide by the
diagonal components of the matrix A. So, we have to always rearrange our
equations so to ensure that we are not dividing by zero. For example the
following system
2x2 + 3x3 = 1
3x1 + 5x2 + 6x3 = 2
9x1 + 2x2 + 3x3 = 3
For this case, a11 = 0, so you will run into problems in the first forward
elimination stage. If you re-write the above equation as
3x1 + 5x2 + 6x3 = 2
2x2 + 3x3 = 1
9x1 + 2x2 + 3x3 = 3
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
42 CHAPTER 3. SYSTEMS OF LINEAR ALGEBRAIC EQUATIONS
will eliminate the problem. The process of interchanging rows so that you
are not dividing by zero is called partial pivoting. In general, if the diagonal
elements are very small, you would like to pivot the rows so that you are not
dividing by a small number.
Exercise 3.1: Solve the following system of 4 equations using Gauss Elimi-
nation
x1 + 2x2 + 3x3 + 4x4 = 3
3x1 + 4x2 + 8x3 + 9x4 = 4
10x1 + 12x2 + 4x3 + 3x4 = 8
5x1 + 6x2 + 7x3 + 8x4 = 10
3.1.2 LU Decomposition
Suppose we want to solve the following linear equation
[A] {X} = {C} . (3.5)
If we can somehow re-express [A] such that
[A] = [L][U ] (3.6)
where [L] and [U ] are lower- and upper-triangular matrices. The structure of [L] and
[U ] are shown below
l11 0 0 0 . . . . . 0
l21 l22 0 0 . . . . . 0
l31 l32 l33 0 . . . . . 0
. . . . . . . . . .
. . . . . . . . . .
[L] =
. . . . . . . . . .
. . . . . . . . . .
. . . . . . . . . .
. . . . . . . . . .
ln1 ln2 ln3 ln4 . . . . . lnn
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
3.1. DIRECT METHODS 43
1 u12 u12 u12 . . . . . u1n
0 1 u23 u24 . . . . . u2n
0 0 1 u34 . . . . . u3n
. . . . . . . . . .
. . . . . . . . . .
[U ] =
. . . . . . . . . .
. . . . . . . . . .
. . . . . . . . . .
. . . . . . . . . .
0 0 0 0 . . . . . 1
If we now substitute Eq. (3.6) into Eq. (3.5), then we will have
[L] [U ] {X} = {C} (3.7)
Let
[U ] {X} = {R} (3.8)
and substitute Eq. (3.8) into Eq. (3.7) to give
[L] {R} = {C} (3.9)
Thus we have replaced the original problem, Eq. (3.5), by two problems given by
Eqs. (3.8) and (3.9). However, Eqs. (3.8) and (3.9) are much easier to solve because
[L] and [U ] are lower- and upper-triangular matrices respectively. This means that to
solve Eqs. (3.8) and (3.9), all you have to do is do backward and forward substitution
respectively.
A strategy of solving Eq. (3.5) using Eqs. (3.8) and (3.9) is the following
• Decompose [A] into [L] and [U]
• Solve Eq. (3.9) and obtain {R} by forward subtitution.
• Use the {R} obtain in the step above and substitute into Eq. (3.8) so that you
can solve for {X} by backward substitution.
The flow of information in the LU decomposition method is shown in Fig. 3.2.
The problem now is how do we find [L] and [R] from [A] ? To illustrate the
procedure of finding [L] and [R], we will look at how it can be done for a small 4 × 4
[A] matrix.
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
44 CHAPTER 3. SYSTEMS OF LINEAR ALGEBRAIC EQUATIONS
[A]{X}={C}
1 decomposition
[L] [U]
[L]{R}={C}
} {R}
2 forward substitution
[U]{X}={R}
}
3 backward substitution
{X}
Figure 3.2: The flow of information in the LU decomposition method.
a11 a12 a13 a14 l11 0 0 0 1 u12 u13 u14
a21 a22 a23 a24 l21 l22 0 0
0 1 u23 u24
a31 =
a32 a33 a34 l31 l32 l33 0 0 0 1 u34
a41 a42 a43 a44 l41 l42 l43 l44 0 0 0 1
If we multiply the rows of [L] by the first column of [U ] and compare it with [A], we
get
l11 = a11 , l21 = a21 , l31 = a31 , l41 = a41 .
So the first column of [L] is the first column of [A]
Multiply the first row of [L] by the columns of [U ] to get
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
3.1. DIRECT METHODS 45
l11 u12 = a12 , l11 u13 = a13 , l11 u14 = a14 .
Hence, the first row of [U ] can be determined to be
a12 a13 a14
u12 = , u13 = , u14 = (3.10)
l11 l11 l11
We keep alternating between getting a column of [L] and a row of [U ]. So we
next get the equation for the second column of [L] by multiplying the rows of [L] by
the second column of [U ].
l21 u12 + l22 = a22
l31 u12 + l32 = a32
l41 u12 + l42 = a42
This gives the 2nd column of [L] to be
l22 = a22 − l21 u12
l32 = a32 − l31 u12 (3.11)
l42 = a42 − l41 u12
Since we already know u12 , l21 , l31 and l41 from Eqs. (3.9) and (3.13). Multiply the
second row of [L] by the columns of [U ] to get
l21 u13 + l22 u23 = a23
l21 u14 + l22 u24 = a24
This allows us to obtain the 2nd row of [U ] to be
a23 −l21 u13
u23 = l22
a23 −l21 u14
u24 = l22
Since Eqs. (3.9), (3.13) and (3.14) will give you l21 , u13 , u14 and l22 . If you keep
going you can get
l33 = a33 − l31 u13 − l32 u23 , l43 = a43 − l41 u13 − l42 u23
u34 = a34 −l31 ul33
14 −l32 u24
l44 = a44 − l41 u14 − l42 u24 − l43 u34
In general, the formula for getting the elements of [L] and [U ] from the elements an
n × n [A] matrix can be written as
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
46 CHAPTER 3. SYSTEMS OF LINEAR ALGEBRAIC EQUATIONS
li1 = ai1 for i = 1, 2, . . . , n
a1j
u1j = l11 for j = 2, 3, . . . , n
Begin for loop j = 2, 3, . . . , n
j−1
P
lij = aij − lik ukj for i = j, j + 1, . . . , n
k=1 (3.12)
j−1
P
ajk − lji uik
ujk = i=1
ljj
for k = j + 1, j + 2, . . . , n
End for loop
uii = 1.0 for i = 1, 2, . . . , n
The LU decomposition algorithm outlined in Eq. (3.12) is also known as the Crout’s
LU decomposition algorithm. Matlab implementation of Crout’s method is shown in
Matlab Program 3.1
After obtaining the [L] and the [U ] matrix, we need to solve for {R}
[L]{R} = {C}
l11 0 0 0 r1
c 1
l21 l22 0 0 r2
c2
l31 l32 l33 0 r3 = (3.13)
c3
l41 l42 l43 l44 r4 c4
by forward substitution. From the first row,
l11 r1 = c1
c1
r1 = (3.14)
l11
From the 2nd row,
l21 r1 + l22 r2 = c2
c2 − l21 r1
r2 =
l22
we can get a numerical value for r2 because we know r1 (from Eq. (3.14)). From the
third and fourth rows of Eq. (3.13) we can get
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
3.1. DIRECT METHODS 47
Matlab Program 3.1 A matlab program to perform LU decomposition using
Crout’s method.
clear all;
A=[1 2 3 4 8;
3 4 8 9 9;
10 12 4 3 1 ;
5 6 7 8 2;
1 2 5 7 9];
[m,n]=size(A);
L=zeros(m,n);
U=zeros(m,n);
for i=1:n
L(i,1)=A(i,1);
end;
for j=2:n
U(1,j)=A(1,j)/L(1,1);
end;
for j=2:n
for i=j:n
sum=0.0;
for k=1:j-1
sum=sum+L(i,k)*U(k,j);
end;
L(i,j)=A(i,j)-sum;
end;
for k=j+1:n
sum=0.0;
for i=1:j-1
sum=sum+L(j,i)*U(i,k);
end;
U(j,k)=(A(j,k)-sum)/L(j,j);
end;
end;
for i=1:n
U(i,i)=1.0;
end;
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
48 CHAPTER 3. SYSTEMS OF LINEAR ALGEBRAIC EQUATIONS
c3 −l31 r1 −l32 r2
r3 = l33
c4 −l41 r1 −l42 r2 −l43 r3
r4 = l44
In general, the formula for ri for a n × n system is
r1 = c1 /l11 !
i−1
ci −
P
lij rj (3.15)
j=1
ri = lii
for i = 2, 3, 4,.....,n
Lastly, we need to solve
[U ]{X} = {R}
1 u12 u13 u14 x1
r1
0 1
u23 u24 x2
r2
0 0 =
1 u34 x3 r3
0 0 0 1 x4 r4
Use the method of back substitution which we was introduced in the method of
Gauss elimination. Start with the last (fourth row)
x 4 = r4 (3.16)
Next, use the 2nd last (third row) to give
x3 + u34 x4 = r3
x3 = r3 − u34 x4 (3.17)
You can get a value for x3 because we know x4 from Eq. (3.16). The second and
first row will give
x2 = r2 − u23 x3 − u24 x4
x1 = r1 − u12 x2 − u13 x3 − u14 x4
which can be use to solve for x2 and x1 (in that sequence). In general, the formula
for a n × n system is
x n = rn
n
xi = ri −
P
uij xj for i = n − 1, n − 2, n − 3, ....., 2, 1 (3.18)
j=i+1
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
3.1. DIRECT METHODS 49
Exercise 3.2: Solve the following system of 4 equations using LU Decom-
position
x1 + 2x2 + 3x3 + 4x4 = 3
3x1 + 4x2 + 8x3 + 9x4 = 4
10x1 + 12x2 + 4x3 + 3x4 = 8
5x1 + 6x2 + 7x3 + 8x4 = 10
In the example above, the LU decomposition was conducted using an [U ] matrix
with ones on the diagonal. This is called the Crout’s method. One can also factorize
a matrix by assuming ones along the diagonal of the [L] matrix instead of the [U ]
matrix. This is called the Doolittle’s method, which for a 4×4 system, would look
something like
a11 a12 a13 a14 1 0 0 0 u11 u12 u13 u14
a21 a22 a23 a24 l21 1 0
0 0 u22 u23
u24
a31 =
a32 a33 a34 l31 l32 1 0 0 0 u33 u34
a41 a42 a43 a44 l41 l42 l43 1 0 0 0 u44
The algorithm to perform Doolittle’s method is very similar to that to do Crout’s
method. The Matlab program to perform Doolittle’s method is shown in Matlab
Program 3.2
Exercise 3.3: Use the Doolittle’s method and find the [L] and [U ] matrices
for the following 3×3 matrix
1 1 −3
3 1 4
2 0 1
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
50 CHAPTER 3. SYSTEMS OF LINEAR ALGEBRAIC EQUATIONS
Matlab Program 3.2 Matlab program to perform LU decomposition using Doolit-
tle’s method.
clear all;
A=[1 2 3 4; 3 4 8 9 ;10 12 4 3 ;5 6 7 8 ];
[m,n]=size(A);
L=zeros(m,n);
U=zeros(m,n);
for j=1:n
U(1,j)=A(1,j);
L(j,j)=1.0;
end;
for i=2:n
L(i,1)=A(i,1)/U(1,1);
end;
for i=2:n
for j=i:n
sum=0.0;
for k=1:i-1
sum=sum+L(i,k)*U(k,j);
end;
U(i,j)=A(i,j)-sum;
end;
for m=i+1:n
sum=0.0;
for nn=1:i-1
sum=sum+L(m,nn)*U(nn,i);
end;
L(m,i)=(A(m,i)-sum)/U(i,i);
end;
end;
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
3.1. DIRECT METHODS 51
3.1.3 Using LU decomposition to find the inverse of a matrix
Sometimes, it is useful to find the inverse of a matrix. LU decomposition can be
used to find inverse of a matrix. The inverse of a nonsingular matrix [A], [A]−1 is
defined as
[A] [A]−1 = [I] (3.19)
where [I] is the identity matrix. I will demonstrate the method for a 4×4 system.
The extension to a n × n system is straightforward. For a 4×4 system, Eq. (3.19)
can be written as
a11 a12 a13 a14 b11 b12 b13 b14 1 0 0 0
a21 a22 a23 a24 b21 b22 b23 b24 0 1 0 0
a31 a32 a33 a34 b31 b32 b33 b34 = 0 0 1 0
a41 a42 a43 a44 b41 b42 b43 b44 0 0 0 1
Observe that one can break the above equations down into the following 4 equations
a11 a12 a13 a14 b 11
1
a21 a22 a23 a24 b21 0
a31 a32 a33 a34 b31 = 0 (3.20)
a41 a42 a43 a44 b41 0
a11 a12 a13 a14 b12 0
a21 a22 a23 a24 b22 1
a31 a32 a33 a34 b32 = 0 (3.21)
a41 a42 a43 a44 b42 0
a11 a12 a13 a14 b 13
0
a21 a22 a23 a24 b23 0
a31 a32 a33 a34 b33 = 1 (3.22)
a41 a42 a43 a44 b43 0
a11 a12 a13 a14 b 14
0
a21 a22 a23 a24 b24 0
a31 a32 a33 a34 b34 = 0 (3.23)
a41 a42 a43 a44 b44 1
Thus solve Eqs. (3.20) - (3.23) will give you all the bij s which are the elements of the
[A]−1 matrix. Note that this method is very efficient because the [A] matrix remain
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
52 CHAPTER 3. SYSTEMS OF LINEAR ALGEBRAIC EQUATIONS
the same. Hence one can use the same [L] and [U ] matrices when solving Eqs. (3.20)
- (3.23).
Note that you do not need to use to use LU Decomposition to solve Eqs. (3.20)
- (3.23). You can use any method you like. As long as you can solve Eqs. (3.20) -
(3.23), you should be able to write down [A]−1 .
Example 3.1: Let’s try to use LU Decomposition to find the inverse of
1 2 3 4
3 4 8 9
[A] =
10 12 4 3
5 6 7 8
Following the discussion above, we will need to define the following 4 equations
1 2 3 4
b11
1
3 4 8 9 b21
0
10 12 4 3 b31 = (3.24)
0
b41 0
5 6 7 8
1 2 3 4
b12
0
3 4 8 9 b22
1
10 12 4 3 b32 = (3.25)
0
b42 0
5 6 7 8
1 2 3 4
b13
0
3 4 8 9 b23
0
10 12 4 3 b33 = (3.26)
1
b43 0
5 6 7 8
1 2 3 4
b14
0
3 4 8 9 b24
0
10 12 4 3 b34 = (3.27)
0
b44 1
5 6 7 8
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
3.1. DIRECT METHODS 53
Using LU decomposition to solve Eq. (3.24)- (3.27) will give you
− 10 9
b11
37
b21 36
=
b31
− 13
18
b41
5
9
7
−9
b12
4
b22 9
=
b32
13
9
b42
− 10
9
1
−3
b13
1
b23 3
=
b33
1
3
b43
− 13
14
9
b14
− 41
b24 36
=
b34
− 25
18
b44
11
9
Thus the inverse of the original matrix can be written down as
10
− 9 − 97 − 13 14 9
37 4 1
−1
36 9 3
− 41
36
[A] = 13 13
1
25
− 18 9 3
− 18
5
9
− 10
9
− 13 11
9
Exercise 3.4: Given that the values of a, b, c and d that satisfies
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
54 CHAPTER 3. SYSTEMS OF LINEAR ALGEBRAIC EQUATIONS
1 1 0 2 a p−1
2 5 2 0 b q−2
3 8 3 1 c = r − 2
5 0 1 3 d s−1
are
a p − 3 + 3q − 2r
b
2 q + p − 5/4 r − 9/4 − 1/4 s
β= =
41 r + 61 − 7/2 p − 15/2 q + 5/8 s
c
8 8
d 1/8 s + 178
− 1/2 p − 5/2 q + 13
8
r
Use the information above to find
−1
1 1 0 2
2 5 2 0
3 8 3 1
5 0 1 3
3.2 Iterative Methods
As the name implies, these methods require iterative procedure in order to obtain
the solutions. Suppose you would like to solve a system that looks like
[A] {x} = {C} (3.28)
where the square brackets denote a matrix and the curly brackets denote a vector.
Suppose that the matrix [A] is rewritten as
[A] = [M ] − [N ]
Substituting into Eq. (3.28) gives
[M ] {x} = [N ] {x} + {C} (3.29)
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
3.2. ITERATIVE METHODS 55
An iterative method of solving for the vector {x} would be
[M ] {x}(k+1) = [N ] {x}(k) + {C} (3.30)
{x}(k) is the current (k’th) guess of the true solution, {x}. Naturally, one would like
to choose the matrix [M ] such that Eq. (3.30) is easier to solve than Eq. (3.28). In
this course, we will cover two different methods,
1. If [M ] consist of only the diagonal elements of [A], then we have the point
Jacobi method.
2. If [M ] consist of the lower triangular elements of [A], then we have the Gauss-
Seidel method.
So one can take the following steps when using iterative techniques to solve linear
systems
1. Define the [M ] and [N ] matrices
2. Obtain an initial guess, {x}(0)
3. Solve Eq. (3.30) to obtain a new guess.
4. Define and obtain the maximum absolute value of the residual vector,
{r} = [A] {x}(k) − {C}
If the maximum absolute value of {r} is greater than a small specified value,
then repeat step 3. Else, end.
3.2.1 Point Jacobi
For this method, the matrix [M ] consist of only the diagonal elements of [A]. The
negative of all other elements is placed in [N ]. This method will be illustrated with
the 3 × 3 system below
a11 a12 a13 x1 c1
a21 a22 a23 x2 = c2
a31 a32 a33 x3 c3
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
56 CHAPTER 3. SYSTEMS OF LINEAR ALGEBRAIC EQUATIONS
[M ] will contain the diagonal elements of [A]
a11 0 0
[M ] = 0 a22 0
0 0 a33
and the -ve of all other elements are placed in the matrix [N ]
0 −a12 −a13
[N ] = −a21 0 −a23
−a31 −a32 0
For this case, the general formula given in Eq. (3.30) will actually look like
(k+1) (k)
a11 0 0 x1 0 −a12 −a13 x1 c1
0 a22 0 x2 = −a21 0 −a23 x2 + c2
0 0 a33 x3 −a31 −a32 0 x3 c3
Solving this system of equation will give you
1
(k+1) (k) (k)
x1 = −a12 x2 − a13 x3 + c1
a11
1
(k+1) (k) (k)
x2 = −a21 x1 − a23 x3 + c2
a22
1
(k+1) (k) (k)
x3 = −a31 x1 − a32 x2 + c3
a33
It is straightforward to see that if [A] is an N × N matrix, then the generic formula
(k+1)
for xi is
N
(k+1) 1 X (k)
xi = −aij x j + c i
aii j=1
j6=i
3.2.2 Gauss-Seidel
This method is the most commonly used iterative method. For conciseness, I will
illustrate the method with the 3 × 3 system below
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
3.2. ITERATIVE METHODS 57
a11 x1 + a12 x2 + a13 x3 = c1 (3.31)
a21 x1 + a22 x2 + a23 x3 = c2 (3.32)
a31 x1 + a32 x2 + a33 x3 = c3 (3.33)
For this method, the matrix [M ] consist of only the lower triangular elements of [A].
[M ] will contain the diagonal elements of [A]
a11 0 0
[M ] = a21 a22 0
a31 a32 a33
The negative of all other elements is placed in [N ].
0 −a12 −a13
[N ] = 0 0 −a23
0 0 0
For this case, the general formula given in Eq. (3.30) will actually look like
(k+1) (k)
a11 0 0 x1 0 −a12 −a13 x1 c1
a21 a22 0 x2 = 0 0 −a23 x2 + c2
a31 a32 a33 x3 0 0 0 x3 c3
Solving this system in sequence from the top row to the bottom row will give you
(k) (k)
(k+1) c1 − a12 x2 − a13 x3
x1 = (3.34)
a11
(k+1) (k)
(k+1) c2 − a21 x1 − a23 x3
x2 = (3.35)
a22
(k+1) (k+1)
(k+1) c3 − a31 x1 − a32 x2
x3 = (3.36)
a33
We can now start the iterative process. WARNING, you have to be careful that none
of your aii ’s are zero !!
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
58 CHAPTER 3. SYSTEMS OF LINEAR ALGEBRAIC EQUATIONS
(0) (0) (0)
1. Assume initial guess values for x1 , x2 , and x3 : x1 , x2 and x3 .
(0) (0)
2. Substitute the guess values for x2 , and x3 into Eq. (3.34) to obtain a better
(1)
estimate for x1 . Call this value x1 .
(1) (0) (1)
3. Substitute x1 and x3 into Eq. (3.35) to obtain a new estimate for x2 = x2 .
(1) (1) (1)
4. Use x1 and x2 in Eq. (3.36) to obtain a better estimate for x3 = x3 .
The flow of information for the Gauss Seidel method is shown in Fig. 3.3.
c1 - a12 x2 - a13 x3
x1 =
a11
c2 - a21x1 - a23 x3
x2 = First iteration
a22
c3 - a31x1 - a32 x2
x3 =
a33
c1 - a12 x2 - a13 x3
x1 =
a11
c2 - a21x1 - a23 x3
x2 = Second iteration
a22
c3 - a31x1 - a32 x2
x3 =
a33
Figure 3.3: The flow of information in the Gauss Seidel method.
Note that one can extend the point Jacobi and Gauss-Seidel methods to solve non-
linear systems too !
Exercise 3.5: Use point Jacobi and the Gauss-seidel iterative method to
solve the following linear system of equations
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
3.2. ITERATIVE METHODS 59
9x1 + 4x2 +x3 = −17
x1 − 2x2 −6x3 = 14
x1 + 6x2 =4
Which method converges faster ?
It is straightforward to see that if [A] is an N × N matrix, then the generic formula
(k)
for xi using the Gauss-Seidel method
j<i N
!
(k+1) 1 X (k+1)
X (k)
xi = − aij xj − aij xj + ci
aii j=1 j>i
3.2.3 Convergence of iterative scheme
Let’s now look at under what conditions will the approximated solution, {x}(k) con-
verge to the exact solution {x}. Define the error matrix at the k’th iteration to
be
{ǫ}(k) = {x}(k) − {x}
Substracting Eq. (3.29) from Eq. (3.30) gives
{ǫ}(k+1) = [M ]−1 [N ] {ǫ}(k)
{ǫ}(k+1) = [P ] {ǫ}(k)
where
[P ] = [M ]−1 [N ]
One should be able to deduce from this equation that the error at the k + 1’th
iteration is related to the initial error
{ǫ}(k+1) = [P ]k+1 {ǫ}(0)
If we assume that the all eigenvalues of the matrix [P ] are linearly independent, then
we can write
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
60 CHAPTER 3. SYSTEMS OF LINEAR ALGEBRAIC EQUATIONS
Matlab Program 3.3 A matlab program containing the Gauss-Seidel method to
solve the equations in Exercise 3.5
clear all
maxiter=1000;
A=[9 4 1;
1 6 0;
1 -2 -6];
C=[-17;
4;
14];
%%%%%%%%%%%%%%%%%%%%%%
%Initial Guess
%%%%%%%%%%%%%%%%%%%%%%%%
x=zeros(size(C));
[rows,columns]=size(A);
for k=1:maxiter
for i=1:rows
sum=0.0;
for j=1:columns
if j~=i
sum=sum-A(i,j)*x(j);
end
end
x(i)=(sum+C(i))/A(i,i);
end
residual=A*x-C
if max(abs(residual)) < 1.0e-8
fprintf(’Solution obtained !’);
x
break;
end;
end;
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
3.2. ITERATIVE METHODS 61
[P ] = [S] [Λ] [S]−1
where the columns of [S] are the eigenvectors of [P ] and [Λ] is a diagonal matrix
whose elements are the eigenvalues of [P ]. Thus we can write the error at the k’th
iteration as
{ǫ}(k+1) = [S] [Λ]k+1 [S]−1 {ǫ}(0)
So in order for the error to go to zero as fast as possible, we would need that the
magnitude of all the eigenvalues of [P ] less than 1.
Exercise 3.6: Consider the system of equations in Exercise 3.5. Define the
[M ] and [N ] matrices. Calculate the eigenvalues of [M ]−1 [N ] and show that they
are consistent with your results in Exercise 3.5.
Exercise 3.7: The equation
[A] {x} = {C} (3.37)
can be solved by iterative schemes which can be written as
[M ] {x}(k+1) = [N ] {x}(k) + {C} (3.38)
{x}(k) is the approximate value of the vector {x} at the k’th iteration. [A] =
[M ] − [N ] are matrices {C} is a constant vector.
1. Consider a system of two equations where
a b
[A] = [M ] − [N ] =
c d
What is the condition among the elements of [A] for the iterative scheme
to converge if
a 0
[M ] = (3.39)
c d
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
62 CHAPTER 3. SYSTEMS OF LINEAR ALGEBRAIC EQUATIONS
2. Solve
1 7 x1 4
= (3.40)
2 3 x2 5
with the iterative scheme and the structure of matrix [M ] given by Eq.
(3.39). Start with
(0)
x1 2
=
x2 0
and write down the the behaviour of your approximated solution for 4
iterations. Can you get a converged solution ?
3. Swap the rows of Eq. (3.40) and try to find the solution again. Can you
get a converged solution ? Explain why.
3.3 Newton’s method for nonlinear equations (re-
visited)
Earlier in the course, we talked about finding roots of two simultaneous nonlin-
ear equations. What if we want to find roots of 3 or more simultaneous nonlinear
equations ? The methods that we have mentioned so far can only cope with linear
equations. What can we do for nonlinear equations ?
For simplicity, the theory for solving a system of 3 or more simultaneous nonlinear
equations will be developed for 3 × 3 system. Extension to larger systems is straight
forward. Consider the following system of 3 simultaneous nonlinear equations
u(x, y, z) = 0
v(x, y, z) = 0
w(x, y, z) = 0
Do a Taylor series expansion on the left hand side and truncate at the first order
terms to get
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
3.3. NEWTON’S METHOD FOR NONLINEAR EQUATIONS (REVISITED) 63
u(xi+1 , yi+1 , zi+1 ) ≈ u(xi , yi , zi )
+ ux (xi , yi , zi )(xi+1 − xi )
+ uy (xi , yi , zi )(yi+1 − yi )
+ uz (xi , yi , zi )(zi+1 − zi )
= 0
v(xi+1 , yi+1 , zi+1 ) ≈ v(xi , yi , zi )
+ vx (xi , yi , zi )(xi+1 − xi )
+ vy (xi , yi , zi )(yi+1 − yi )
+ vz (xi , yi , zi )(zi+1 − zi )
= 0
w(xi+1 , yi+1 , zi+1 ) ≈ w(xi , yi , zi )
+ wx (xi , yi , zi )(xi+1 − xi )
+ wy (xi , yi , zi )(yi+1 − yi )
+ wz (xi , yi , zi )(zi+1 − zi )
= 0
where
∂u ∂u ∂u
ux = , uy = , uz =
∂x ∂y ∂x
∂v ∂v ∂v
vx = , vy = , vz =
∂x ∂y ∂z
∂w ∂w ∂w
wx = , wy = , wz =
∂x ∂y ∂z
We can put the above system of equations into matrix form
ux uy uz xi+1 − xi u
vx vy vz yi+1 − yi =− v (3.41)
wx wy wz zi+1 − zi w
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
64 CHAPTER 3. SYSTEMS OF LINEAR ALGEBRAIC EQUATIONS
Equation (3.41) is a linear equation so we should be able to use any of the methods
discussed earlier to solve it !
Hence, to solve the original system of equations, do the following
1. Calculate all the required derivatives.
2. Guess the (initial) values for xi , yi and zi .
3. Solve Eq. (3.41) to obtain (xi+1 − xi ), (yi+1 − yi ) and (zi+1 − zi ). You can use
Gauss elimination, LU decomposition, Gauss-Seidel etc.
4. Calculate xi+1 , yi+1 , zi+1 since you know xi , yi and zi from step 2 above.
5. Repeat steps 2, 3 and 4 until a converged solution is obtained.
This method is a more general version of the Newton-Raphson method mentioned
earlier.
Example 3.2:
If it is desired to solve the following set of nonlinear algebraic equations
1
3x + cos(yz) =
2
2 2
x − 81(y + 0.1) + sin(z) = −1.06
10π − 3
e−xy + 20z + = 0
3
using the Newton-Raphson method, one would have to define the functions u(x, y, z),
v(x, y, z), w(x, y, z) and calculate the corresponding partial derivatives. Firstly, let
1
u(x, y, z) = 3x + cos(yz) −
2
v(x, y, z) = x − 81(y + 0.1)2 + sin(z) + 1.06
2
10π − 3
w(x, y, z) = e−xy + 20z + .
3
The corresponding partial derivatives are
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
3.3. NEWTON’S METHOD FOR NONLINEAR EQUATIONS (REVISITED) 65
∂u ∂u ∂u
=3 ∂y
= −z sin(yz) = −y sin(yz)
∂x ∂z
∂v ∂v ∂v
= 2x ∂y
= −16.2 − 162y = cos(z)
∂x ∂z
∂w ∂w ∂w
= −ye−xy ∂y
= −xe−xy = 20
∂x ∂z
In order to use Eq. (3.41), we would need to guess initial values of x, y and z.
Starting with (x0 , y0 , z0 ) = (0, 0, 0) and substituting these numbers into Eq. (3.41)
gives
3.0 0.0 0.0 x 1 − x 0
−0.5000000000
0.0 −16.2 1.0 y − y = −0.25
1 0
0.0 0.0 20.0 z −z
1 0 −10.47197551
3.0 0.0046301 0.00014934 x2 − x1
0.000039
−0.33334 −13.464 0.86602 y2 − y1 = −0.02827
z −z
0.016841 0.16619 20.0 2 1 0.0028
3.0 0.0040504 0.00011437 x3 − x2
0.00000061
−0.33331 −13.805 0.86608 y3 − y2 = 0.000359
z −z
0.014744 0.166246 20.0 3 2 −0.0000
and so on. Note that with every iteration, the right hand side becomes smaller and
smaller. Eventually, the right hand side becomes a zero vector and xn+1 = xn =
−0.16667, yn+1 = yn = −0.01481, zn+1 = zn = −0.523476 and we get a converged
solution when n is large.
Exercise 3.8: Solve the following set of nonlinear equations
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
66 CHAPTER 3. SYSTEMS OF LINEAR ALGEBRAIC EQUATIONS
xyz − x2 + y 2 = 1.34
xy − z 2 = 0.09
ex − ey + z = 0.41
using the multi-dimensional Newton-Raphson method.
One can easily extend the example above to a system of N equations. Suppose
it is desired to solve
f1 (x1 , x2 , x3 , x4 , ......., xN ) = 0
f2 (x1 , x2 , x3 , x4 , ......., xN ) = 0
f3 (x1 , x2 , x3 , x4 , ......., xN ) = 0
,
.
.
fN (x1 , x2 , x3 , x4 , ......., xN ) = 0
Then Eq. (3.41) extended to this system is
(k−1) (k−1) (k−1)
∂f1 ∂f1
··· ∂f1 (k) (k−1)
f 1 x , x , . . . , x
∂x1 ∂x2 ∂xN
x1 − x1
1 2 N
∂f2 ... ..
x(k) (k−1) f2 x(k−1) , x(k−1) , . . . , x(k−1)
2 − x2
∂x1
. 1 2 N
.. .. ..
.. =− .
. . . .
..
∂fN ∂fN (k) (k−1)
∂x1
··· ··· ∂xN
xN − xN
fN x1(k−1) , x2(k−1) , . . . , xN
(k−1)
(k−1) (k)
xi is the previous guess of xi and xi is your next (hopefully better !) guess of
xi . The matrix of partial derivatives on the left hand side, ∂fi /∂xj , is usually known
(k−1)
as the Jacobian. It is evaluated at xi .
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
Chapter 4
Least Squares
Suppose from an experiment you get a set of x, y values that look something like in
Fig. 4.1. The data looks like it can be well approximated by a straight line. How
can you draw a line a best fit through the data ? I.e. which of lines p, q or r above
best fit the data ? A method of systematically calculating the line of best fit is called
Least Squares approximation.
y(x)
p
q
Figure 4.1: Fitting a straight lines through a set of data points.
67
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
68 CHAPTER 4. LEAST SQUARES
4.1 Linear least squares approximation
Let’s represent any one of the lines in Fig. 4.1 as
Y = ax + b (4.1)
So our task now is to find a and b such that the predicted by Eq. (4.1) above is
as close to the actual value of yi at a particular value of xi . The error between the
predicted y value and yi is denoted by
ei = yi − (axi + b)
|{z} | {z }
actual value predicted value
Because we have N data points, let’s square the error at every point and sum
the all up. That way we have the total error between the line and the actual data
points.
S = e21 + e22 + e23 + e24 + ....... + e2N
X N
= e2i
i=1
N
X
S= (yi − axi − b)2
i=1
We would like to find a and b, such that S is minimized. So we differentiate S with
respect to a and b and the derivatives to zero. Remember that xi and yi are constants
(data points) when we do the differentiation.
N
P
∂S
∂a
= 2(yi − axi − b)(−xi ) = 0
i=1
N (4.2)
∂S
P
∂b
= 2(yi − axi − b)(−1) = 0
i=1
Exercise 4.1: Show that Eq. (4.2) can be re-written as
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
4.2. POLYNOMIAL LEAST SQUARES APPROXIMATION 69
N
X N
X N
X
a x2i +b xi = xi yi
i=1 i=1 i=1
N
X N
X
a xi + bN = yi
i=1 i=1
The above two equations are just a system of 2 × 2 linear equations. They can
be solved to get a and b.
N N
N
P 2 P P
xi xi xi yi
i=1 i=1 a i=1
P N
b = N
P
xi N
yi
i=1 i=1
Exercise 4.2: Example: Suppose that you have done an experiment and
collected the following data on how the resistance of a wire varies with tempera-
ture
T (o C) R (ohms)
20.5 765
32.7 826
51.0 873
73.2 942
95.7 1032
Obtain a linear least squares fit of the data.
4.2 Polynomial least squares approximation
Suppose we have data as shown in Fig. 4.2. Clearly, it is not appropriate to fit a
straight line through the data. For this data, it is more appropriate to fit a polynomial
through it. Assume a functional relationship of the form
Y = a0 + a1 x + a2 x2 + ...... + an xn
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
70 CHAPTER 4. LEAST SQUARES
1.005
1.004
1.003
1.002
1.001
1
y
0.999
0.998
0.997
0.996
0.995
0 0.2 0.4 0.6 0.8 1
Figure 4.2: Data points requiring nonlinear fit of data.
The error between the actual data and and the polynomial above is defined as
e i = y i − Yi
= yi − a0 + a1 xi + a2 x2i + ...... + an xni
= yi − a0 − a1 xi − a2 x2i − ...... − an xni
Let’s square the error and sum it
N
X N
X 2
S= e2i = yi − a0 − a1 xi − a2 x2i − ...... − an xni (4.3)
i=1 i=1
Equation (4.3) represent the total error. We would like to minimize the. To do this,
calculate the partial derivatives with respect to ai ’s and set them to zero.
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
4.2. POLYNOMIAL LEAST SQUARES APPROXIMATION 71
N
P
∂S
∂a0
= 2 (yi − a0 − a1 xi − a2 x2i − ...... − an xni ) (−1) = 0,
i=1
PN
∂S
∂a1
= 2 (yi − a0 − a1 xi − a2 x2i − ...... − an xni ) (−xi ) = 0,
i=1
PN
∂S
= 2 (yi − a0 − a1 xi − a2 x2i − ...... − an xni ) (−x2i ) = 0,
∂a2
i=1 (4.4)
.
.
.
N
P
∂S
∂an
= 2 (yi − a0 − a1 xi − a2 x2i − ...... − an xni ) (−xni ) = 0,
i=1
Equation (4.4) can be written in the following matrix form
P P 2 P 3 P n P
PN P 2 x i P 3 x i P x4i ··· P xn+1 i
a0
P yi
P x2i P x3i P x4i P x5i ··· x a1 xi yi
P n+2i
P x2 yi
P x3i P x4i P x5i P x6i ··· x a2
P n+3i
= P 3i
xi xi xi xi ··· xi
a3 xi yi
.. .. .. .. .. .. ..
..
P. n P . n+1 P . n+2 P . n+3 . P . 2n . P .n
xi xi xi xi ··· xi an xi yi
The above matrix system can then be solved using techniques such as Gauss elimi-
nation or LU decomposition to give the ai ’s.
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
72 CHAPTER 4. LEAST SQUARES
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
Chapter 5
Interpolation
Sometimes, you may need to estimate the value between data points. One of the
more common method is polynomial interpolation. We fit the data with the function
f (x) = a0 + a1 x + a2 x2 + · · · + an xn (5.1)
and require that the function passes through all the data points. See examples shown
in Fig. 5.1.
f(x) f(x) f(x)
x x x
Figure 5.1: Polynomial interpolation through a set of data points.
73
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
74 CHAPTER 5. INTERPOLATION
Shape of actual function f(x)
d
f(x1) b
f(x)
c
Linear interpolated function f1(x)
f(x0)
a
x0 x x1
Figure 5.2: Linear Newton Interpolation.
5.1 Newton Interpolation
5.1.1 Linear Newton Interpolation
For linear interpolation, we simply draw a straight line between the two given data
points a and b (see Fig. 5.2). If we want to estimate the value at point x, linear
interpolation will give us the value at c whereas the true value is d. Thus there is an
error associated with the linear interpolation.
From similar triangles,
f1 (x) − f (x0 ) f (x1 ) − f (x0 )
= (5.2)
x − x0 x1 − x0
Exercise 5.1: Rearrange Eq. (5.2) to give
f (x1 ) − f (x0 )
f1 (x) = f (x0 ) + (x − x0 ) (5.3)
x1 − x0
and
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
5.1. NEWTON INTERPOLATION 75
x1 − x x − x0
f1 (x) = f (x0 ) + f (x1 ) (5.4)
x1 − x0 x1 − x0
A few things worth noting about the linear interpolation
• Equation (5.3) is called the linear-interpolation formula
• The term (f (x1 ) − f (x0 ))/(x1 − x0 ) in Eq. (5.3) is called the finite divided
difference approximation of the first derivative (see later).
• In general, the smaller the interval between x0 and x1 , the better the approxi-
mation at point x.
• The subscript 1 in f1 (x) denote a first order interpolating polynomial.
Exercise 5.2: Use linear interpolation to estimate the value of the function
f (x) = 1 − e−x
at x = 1.0. Use the interval x0 = 0 and x1 = 5.0. Repeat the calculation with
x1 = 4.0, x1 = 3.0 and x1 = 2.0. Illustrate on a graph how you are approaching
the exact value of f (1) = 0.6321.
5.1.2 Quadratic Newton Interpolation
For quadratic interpolation, let’s assume that the polynomial is of the following form
f2 (x) = b0 + b1 (x − x0 ) + b2 (x − x0 )(x − x1 ) (5.5)
The interpolated function is required to have the the same values as the actual
function at the data points, hence
• f2 (x0 ) = f (x0 )
• f2 (x1 ) = f (x1 )
• f2 (x2 ) = f (x2 )
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
76 CHAPTER 5. INTERPOLATION
Shapeofactualfunctionf(x)
Quadraticallyinterpolated
functionf2(x)
f(x2)
f(x1)
f(x0)
x0 x1 x2
Figure 5.3: Quadratic Newton Interpolation.
Note that the subscript 2 in f2 (x) denote that this is a second order interpolating
polynomial. The following exercise shows that Eq. (5.5) is just another form of the
polynomial function defined in Eq. (5.1).
Exercise 5.3: Show that Eq. (5.5) can be re-written as
f2 (x) = a0 + a1 x + a2 x2
where
a0 = b0 − b1 x0 + b2 x0 x1
a1 = b1 − b2 x0 − b2 x1
a2 = b2
Hence, in order to do quadratic interpolation, all we need to do is to find all the
b’s in Eq. (5.5).
(1) Set x = x0 in Eq. (5.5). Therefore
b0 = f (x0 ) (5.6)
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
5.1. NEWTON INTERPOLATION 77
(2) If we now let x = x1 in Eq. (5.5) and use the result from step (1) above, to get
f (x1 ) − f (x0 )
b1 = (5.7)
x1 − x0
(3) Since f2 (x2 ) = f (x2 ), we can use Eqs. (5.6) and (5.7) together with Eq. (5.5)
to obtain
f (x2 )−f (x1 ) f (x1 )−f (x0 )
x2 −x1
− x1 −x0
b2 = (5.8)
x2 − x0
Equations (5.7) and (5.8) can be simplified by introducing the following notation
f (x1 ) − f (x0 )
b1 = f [x1 , x0 ] = (5.9)
x1 − x0
f [x2 , x1 ] − f [x1 , x0 ]
b2 = f [x2 , x1 , x0 ] = (5.10)
x2 − x0
The advantage of using the square brackets will be clear when we look at the general
form of Newton’s interpolating polynomials in the next section.
5.1.3 General form of Newton’s Interpolating Polynomials
In general, if you have n data points, you will fit a polynomial of order n − 1 through
all the n data points.
fn (x) = b0 + b1 (x − x0 ) + b2 (x − x0 )(x − x1 ) + · · · + bn (x − x0 )(x − x1 ) · · · (x − xn−1 )
where all the b’s are defined as
b0 = f (x0 )
b1 = f [x1 , x0 ]
b2 = f [x2 , x1 , x0 ]
.. .. ..
. . .
bn = f [xn , xn−1 , · · · , x2 , x1 , x0 ]
where
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
78 CHAPTER 5. INTERPOLATION
f (xi ) − f (xj )
f [xi , xj ] =
xi − xj
is called the first divided difference,
f [xi , xj ] − f [xj , xk ]
f [xi , xj , xk ] =
xi − xk
is called the second divided difference
f [xi , xj , xk ] − f [xj , xk , xl ]
f [xi , xj , xk , xl ] =
xi − xl
is called the third divided difference and
f [xn , xn−1 , . . . , x1 ] − f [xn−1 , xn−2 , . . . , x0 ]
f [xn , xn−1 , . . . , x1 , x0 ] =
xn − x0
is called the n’th divided difference. Notice that the second divided difference is
calculated from the first divided difference. The third divided difference is calculated
from the second divided difference and so on. This is illustrated in Fig. 5.4.
i xi f ( xi ) First Second Third
0 x0 f ( x0 ) f [ x1 , x0 ] f [ x2 , x1 , x0 ] f [ x3, x2, x1, x0 ]
1 x1 f ( x1 ) f [ x2 , x1 ] f [ x3 , x2 , x1 ]
]
2 x2 f ( x2 ) f [ x3 , x2 ]
3 x3 f ( x3 )
Figure 5.4: Calculation of Divided Difference.
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
5.2. LAGRANGE INTERPOLATING POLYNOMIALS 79
f(x1)
f(x0)
x0 x1
Figure 5.5: Two data points
Exercise 5.4: Use the Newton’s interpolating polynomial to fit a function
through
f (x) = 1 − e−x .
Assume that you have data at points x =0, 1, 3 and 8.
5.2 Lagrange Interpolating Polynomials
Suppose we want to interpolate between two data points, (x0 , f (x0 )), (x1 , f (x1 ))
shown in Fig. 5.5.
Assume that we have two functions, L0 (x) and L1 (x). They are defined such that
1 if x = x0
L0 (x) =
0 if x = x1
0 if x = x0
L1 (x) =
1 if x = x1
It should be easy to convince yourself that L1 (x) and L2 (x) can be written as
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
80 CHAPTER 5. INTERPOLATION
L0(x) L1(x)
x0 x1
Figure 5.6: Illustrative sketch of L0 (x) and L1 (x).
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
5.2. LAGRANGE INTERPOLATING POLYNOMIALS 81
f1(x)=f(x0)L0(x)+f(x1)L1(x)
f(x1)
f(x0)L0(x)
f(x1)L1(x)
f(x0)
x0 x1
Figure 5.7: Illustrative sketch of Lagrange interpolation polynomial.
(x − x1 )
L0 (x) =
(x0 − x1 )
(x − x0 )
L1 (x) =
(x1 − x0 )
These two functions are sketched in Fig. 5.6. The first order Lagrange interpo-
lating polynomial is obtained by a linear combination of L0 (x) and L1 (x).
f1 (x) = f (x0 )L0 (x) + f (x1 )L1 (x) (5.11)
A sketch of f1 (x) is shown in Fig. 5.7.
Suppose we three data points, (x0 , f (x0 )), (x1 , f (x1 )), (x2 , f (x2 )) and we want to
fit a polynomial through it (see Fig. 5.8). Assume that we have three functions,
L0 (x), L1 (x) and L2 (x). They are defined such that
1 if x = x0
L0 (x) = 0 if x = x1 (5.12)
0 if x = x2
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
82 CHAPTER 5. INTERPOLATION
f(x2)
f(x1)
f(x0)
x0 x1 x2
Figure 5.8: Three points for polynomial Lagrange interpolation.
0 if x = x0
L1 (x) = 1 if x = x1 (5.13)
0 if x = x2
0 if x = x0
L2 (x) = 0 if x = x1 (5.14)
1 if x = x2
It should be easy to convince yourself that the following expressions for L0 (x), L1 (x)
and L2 (x) satisfies the conditions listed in Eqs. (5.12) - (5.14).
(x − x1 )(x − x2 )
L0 (x) = (5.15)
(x0 − x1 )(x0 − x2 )
(x − x2 )(x − x0 )
L1 (x) = (5.16)
(x1 − x2 )(x1 − x0 )
(x − x1 )(x − x0 )
L2 (x) = (5.17)
(x2 − x1 )(x2 − x0 )
For illustrative purposes, Eqs. (5.15) - (5.17) are shown in Fig. 5.9.
The second order Lagrange interpolating polynomial is obtained by a linear com-
bination of L0 (x), L1 (x) and L2 (x).
f2 (x) = f (x0 )L0 (x) + f (x1 )L1 (x) + f (x2 )L2 (x) (5.18)
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
5.2. LAGRANGE INTERPOLATING POLYNOMIALS 83
1.5
L1(x)
11
L0(x)
0.5 L2(x)
00
-0.5
0
x10 2 3
x1
4 5
x62 7 8
Figure 5.9: Graphical illustration of Eqs. (5.15) - (5.17).
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
84 CHAPTER 5. INTERPOLATION
In general, the nth order Lagrange polynomial (i.e. the Lagrange polynomial
that can be fit through n + 1 data points) can be represented concisely as
n
X
fn (x) = Li (x)f (xi ) (5.19)
i=0
where
n
Y x − xj
Li (x) = (5.20)
x − xj
j=0 i
j6=i
Exercise 5.5: Given the following data
x0 = 1 f (x0 ) = 2
x1 = 4 f (x1 ) = 3.386294
x2 = 6 f (x2 ) = 3.791760
Estimate the value of f (x = 2) using both first and second order Lagrange
Interpolating polynomials. The matlab program solution for this exercise is given
in Matlab Program 5.1
Note that
• The Lagrange interpolating polynomial is just a different form of the Newton
interpolating polynomials.
• The main advantage is that it is slightly easier to program.
• It is slightly slower to compute than the Newton Interpolating polynomials.
Exercise 5.6: Show that both the Lagrange and Newton interpolating poly-
nomials are identical by reducing the first and second order Newton polynomials
to first and second order Lagrange polynomials respectively.
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
5.3. SPLINE INTERPOLATION 85
Matlab Program 5.1 Matlab program for Exercise 5.5.
clear all;
xi=[1 4 6];
fxi=[2 3.386294 3.791760];
x=2
N=length(xi);
fx=0;
for i=1:N
Lix=1.0;
for j=1:N
if j ~= i
Lix=Lix*(x-xi(j))/(xi(i)-xi(j));
end
end
fx=fx+fxi(i)*Lix;
end
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
86 CHAPTER 5. INTERPOLATION
Cubic spline interpolation
20
15
f(x)
10 Lagrange interpolation
5 Exact function
Data points
0
0 2 4 6 8
x
Figure 5.10: Illustration of the advantage of spline interpolation.
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
5.3. SPLINE INTERPOLATION 87
(a) (b)
S(x)
fi(x)isasingle
polynomialfunction
S0(x) S1(x) S (x)
2
x x
Figure 5.11: Newton or Lagrange interpolation polynomial (a) and Spline function
(b).
5.3 Spline Interpolation
Newton or Lagrange interpolation can lead to erroneous results when there is an
abrupt change in data (see example in Fig. 5.10)
Splines are made up of piecewise polynomials connecting only two data points. This
is different to Newton or Lagrange polynomials where you have only one polynomial
connecting all data points. See Fig. 5.11.
• The spline function,S(x), is made up of several polynomial functions, Si (x).
• S(x) = Si (x) for xi < x < xi+1 .
• Si (x) is ONLY defined between xi+1 and xi . Outside of the interval, Si (x) is
not defined and has got NO meaning.
• If Si (x) are linear functions, then S(x) is called a linear spline.
Si (x) quadratic → S(x) quadratic spline.
Si (x) cubic → S(x) cubic spline
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
88 CHAPTER 5. INTERPOLATION
S2(x)
S1(x)
S0(x)
f(x3)
f(x1)
f(x2)
f(x0)
x0 x1 x2 x3
Figure 5.12: Linear spline.
5.3.1 Linear splines
• Si (x) = ai + bi (x − xi )
• S(x) = Si (x) for xi < x < xi+1 .
• Note that S(x) is continuous but the derivative of S(x), S ′ (x), is not continuous
at x = xi .
• If there are n + 1 data points, there are only n intervals and hence, n number
of defined Si (x).
Need to find all the ai ’s and bi ’s in order to find S(x). For this case, there are
2n number of unknowns. In order to find all the unknowns, we need to impose the
following requirements
1. Require S(x) to have the values of f (xi ) at x = xi , hence,
Si (xi ) = ai = f (xi ) (5.21)
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
5.3. SPLINE INTERPOLATION 89
2. Require that the function values at adjacent polynomials must be equal at the
interior knots.
Si (xi+1 ) = Si+1 (xi+1 )
ai + bi (xi+1 − xi ) = ai+1
Using Eq. (5.21) we conclude that
f (xi+1 ) − f (xi )
bi = (5.22)
xi+1 − xi
5.3.2 Quadratic spline
Spline made up of quadratic elements is shown in Fig. 5.14. It is important to note
that
• Si (x) = ai + bi (x − xi ) + ci (x − xi )2
• S(x) = Si (x) for xi < x < xi+1 .
• S(x) is continuous and the derivative of S(x), S ′ (x), is also continuous at
x = xi .
• If there are n + 1 data points, there are only n number of defined Si (x).
We need to find all the ai ’s, bi ’s and ci ’s in order to completely define S(x). To
do this we will require that
1. S(x) to have the values of f (xi ) at x = xi , hence,
Si (xi ) = ai = f (xi ) (5.23)
2. the function values at adjacent polynomials must be equal at the interior knots.
Si (xi+1 ) = Si+1 (xi+1 )
ai + bi (xi+1 − xi ) + ci (xi+1 − xi )2 = ai+1
For the sake of conciseness, let hi = xi+1 − xi , so
ai + bi hi + ci h2i = ai+1 (5.24)
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
90 CHAPTER 5. INTERPOLATION
S2(x)
S0(x) S1(x)
f(x3)
f(x1)
f(x2)
f(x0)
x0 x1 x2 x3
Figure 5.13: Quadratic spline.
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
5.3. SPLINE INTERPOLATION 91
′
3. the derivative of Si (x), Si (x), to be continuous at the interior nodes
′ ′
Si (xi+1 ) = Si+1 (xi+1 ) (5.25)
Equation (5.25) leads to
bi + 2ci hi = bi+1 (5.26)
From Eq. (5.24) we can get
ai+1 − ai
bi = − ci h i (5.27)
hi
Exercise 5.7: Substitute Eq. (5.27) into (5.26) and show that
1 aj+1 1 1 aj−1
cj = − aj + + − cj−1 hj−1 (5.28)
hj hj hj−1 hj hj−1
Hence, to construct quadratic splines, do the following
1. Make use of Eq. (5.23) and set ai = f (xi ).
2. Assume c0 = 0. This is equivalent as saying that the second derivative at x0 is
zero !
3. Use Eq. (5.28) to obtain all other values of ci ’s.
4. Use Eq. (5.27) to obtain the bi ’s.
After steps 1-4, we can evaluate the function values at the point you are interested
from the following formula Si (x) = ai + bi (x − xi ) + ci (x − xi )2 . The only tricky thing
left is to figure out which interval your x value belongs to.
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
92 CHAPTER 5. INTERPOLATION
S2(x) S3(x)
S1(x)
S0(x)
f(x3)
f(x1)
f(x2)
f(x0)
x0 x1 x2 x3
Figure 5.14: Quadratic/Cubic spline.
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
5.3. SPLINE INTERPOLATION 93
Exercise 5.8:
You are given the four data points shown in Fig. 5.14 and asked to fit a quadratic
spline through the four data points. The quadratic spline can be thought as
being made up of three ’real’ quadratic polynomial S0 (x), S1 (x), S2 (x) and one
’imaginary’ polynomial S3 (x). The three polynomials can be written as
S0 (x) = a0 + b0 (x − x0 ) + c0 (x − x0 )2 (5.29)
S1 (x) = a1 + b1 (x − x1 ) + c1 (x − x1 )2 (5.30)
S2 (x) = a2 + b2 (x − x2 ) + c2 (x − x2 )2 (5.31)
S3 (x) = a3 + b3 (x − x3 ) + c3 (x − x3 )2 (5.32)
(a) Show that
ai = f (xi ) for i =0, 1, 2, 3 (5.33)
(b) By requiring that adjacent Si (x) to have common values at the interior
points, show that
a1 − a0
b0 = − c0 h 0 (5.34)
h0
a2 − a1
b1 = − c1 h 1 (5.35)
h1
a3 − a2
b2 = − c2 h 2 (5.36)
h2
where hi = xi+1 − xi .
(c) By requiring that the derivative of the adjacent Si (x) to have common values
at the interior points, show that
b0 + 2c0 h0 = b1 (5.37)
b1 + 2c1 h1 = b2 (5.38)
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
94 CHAPTER 5. INTERPOLATION
(d) Use Eqs. (5.34) to (5.38) together with the assumption that c0 = 0 and show
that the ci ’s could be obtained by solving the following matrix equation
0
1 0 0 c0
a2 a0
− a1 h10 + 1
h 0 h 1 0 c1 = h1 h1
+ h0 (5.39)
0 h1 h2 c2 a3
− a2 h11 + 1
+ a1
h2 h2 h1
(e) Explain how you can obtain the coefficients for S0 (x), S1 (x), S2 (x) from the
results in parts (a)-(d).
In order to illustrate the implementation of the numerical algorithm, consider the
set of data points
xi f (xi )
0.0 5.0
1.0 2.0
4.0 4.0
6.0 9.0
8.0 2.0
10.0 0.0
The sample matlab program to implement quadratic spline with the data set above
is shown in Matlab program 5.2. The data points are shown in circles and the solid
line represents the quadratic spline which was obtain by evaluating S(x) at 100 data
points between 0 and 10. Note that the first spline is linear because we have set
c0 = 0. Note also that the spline must pass through all data points.
5.3.3 Cubic Splines
Figure 5.16 shows a cubic spline. At each interval, the splines now consist of a cubic
polynomial. Cubic splines are usually smoother than quadratic splines. Note that
• Si (x) = ai + bi (x − xi ) + ci (x − xi )2 + di (x − xi )3
• S(x) = Si (x) for xi < x < xi+1 .
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
5.3. SPLINE INTERPOLATION 95
Matlab Program 5.2 A matlab program implementing the algorithm for quadratic
spline.
clear all
x_data=[0 1 4 6 8 10 ]
f_data=[5 2 4 9 2 0]
N=length(x_data)-1; %number of intervals
for i=1:N+1
a(i)=f_data(i);
end
c(1)=0.0;
h(1)=x_data(2)-x_data(1);
for j=2:N
h(j)=x_data(j+1)-x_data(j);
c(j)=(1/h(j))*(a(j+1)/h(j)-a(j)*((1/h(j-1))+(1/h(j)))+a(j-1)/h(j-1)-c(j-1)*h(j-1));
end
for i=1:N
b(i)=(a(i+1)-a(i))/h(i)-c(i)*h(i);
end
x=linspace(x_data(1),x_data(N+1),100);
for i=1:length(x)
for j=1:N
if x_data(j)>x(i)
j=j-1;
break
end
end
f(i)=a(j)+b(j)*(x(i)-x_data(j))+c(j)*(x(i)-x_data(j))*(x(i)-x_data(j));
end
plot(x_data,f_data,’ko’,x,f);
xlabel(’x’);ylabel(’f(x)’);
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
96 CHAPTER 5. INTERPOLATION
10
4
f(x)
−2
−4
0 1 2 3 4 5 6 7 8 9 10
x
Figure 5.15: The output from Matlab program 5.2. quadratic spline, ◦ data
points.
• Note that S(x) is continuous, the derivative of S(x), S ′ (x), is continuous and
the second derivative of S(x), S ′′ (x), is also continuous at x = xi .
• If there are n + 1 data points, there are only n number of defined Si (x).
We need to find all the ai ’s, bi ’s, ci ’s and di ’s (i = 0..n − 1) in order to completely
determine S(x). For reasons that will become apparent later, we need to ”introduce”
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
5.3. SPLINE INTERPOLATION 97
S2(x)
S0(x) S1(x)
f(x3)
f(x1)
f(x2)
f(x0)
x0 x1 x2 x3
Figure 5.16: Cubic spline.
an and cn in order to completely determine S(x). Hence, there are n + 1 number of
ai ’s and ci ’s (i = 0..n) but only n number of bi and di (i = 0..n − 1) . The unknowns
that need to be solve for in order to completely determine cubic splines are
a0 , a1 , ...., an
b0 , b1 , ...., bn−1
c0 , c1 , ...., cn
d0 , d1 , ...., dn−1
Note that there are n + 1 ai ’s and ci ’s but only n number of bi ’s and di ’s. Thus the
total number of unknowns are (n + 1) + n + (n + 1) + n = 4n + 2. In order to find
all the unknowns, we need 4n + 2 equations. To find these equations,
1. S(x) have the values of f at x = xi , hence,
Si (xi ) = ai = f (xi ) (5.40)
where i = 0..n. Equation (5.40) will give you n + 1 equations.
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
98 CHAPTER 5. INTERPOLATION
2. the function values at adjacent polynomials must be equal at the interior knots.
Si (xi+1 ) = Si+1 (xi+1 )
ai + bi (xi+1 − xi ) + ci (xi+1 − xi )2 + di (xi+1 − xi )2 = ai+1
For the sake of conciseness, let hi = xi+1 − xi , so
ai + bi hi + ci h2i + di h3i = ai+1 (5.41)
where i = 0..n − 1. This will give us another n number of equations.
3. We also want the derivative of Si (x) to be continuous at the interior nodes
′ ′
Si (xi+1 ) = Si+1 (xi+1 ) (5.42)
which leads to
bi + 2ci hi + 3di h2i = bi+1 (5.43)
where i = 0..n−2. Note that for most cubic splines, bn does not exist. Equation
(5.43) will give us another n − 1 equations.
4. the double derivative of Si (x) to be continuous at the interior nodes.
′′ ′′
Si (xi+1 ) = Si+1 (xi+1 ) (5.44)
which gives
ci + 3di hi = ci+1 (5.45)
where i = 0..n − 1. This will give us another n equations.
The four conditions above will give us (n + 1) + n + (n − 1) + n = 4n number
of equations (Eqs. (5.40), (5.41), (5.43), (5.45)). We are still in need of two more
equations because there is a total of 4n + 2 number of unknowns. If you know more
information about the function you are trying to approximate (e.g. its derivative
at the end points etc.) you can use this information to generate the two additional
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
5.3. SPLINE INTERPOLATION 99
information. In the absence of any information, it is common to require that the
double derivative of the cubic spline is zero at the two end points. This will lead to
the following two equations
c0 = 0 (5.46)
cn = 0 (5.47)
These two additional equations will give us a total of 4n + 2 equations. So it is now
possible to solve for all the ai ’s, bi ’s, ci ’s and di ’s.
In order to solve for all the unknowns, do the following. Re-arrange (5.45) to get
ci+1 − ci
di = (5.48)
3hi
Put Eq. (5.48) into Eq. (5.43) to get
bi+1 = bi + (ci + ci+1 ) hi (5.49)
Put Eq. (5.48) into Eq. (5.41) to get
1 hi
bi = (ai+1 − ai ) − (2ci + ci+1 ) (5.50)
hi 3
Exercise 5.9: Substitute Eq. (5.50) into Eq. (5.49) and show that you will
get the following relationship between ai ’s and ci ’s
3 3
hj−1 cj−1 + 2 (hj + hj−1 ) cj + hj cj+1 = (aj+1 − aj ) + (aj−1 − aj ) (5.51)
hj hj−1
where j = 1..n − 1.
Equation (5.51) is a tridiagonal system and can be solve if we assume that c0 =
cn = 0. The linear system of equation equation that we need to solve looks like the
following
[A]{X} = {C} (5.52)
where
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
100 CHAPTER 5. INTERPOLATION
1 0 0 0 ··· 0
h0 2 (h0 + h1 ) h1 0 ··· 0
..
0 h1 2 (h1 + h2 ) h2 . 0
[A] = .. .. .. ..
0
0 . . . .
0 0 ··· hn−2 2 (hn−2 + hn−1 ) hn−1
0 0 0 ··· 0 1
c0
c1
..
.
{X} = ..
.
..
.
cn
0
3 3
h1
(a2 − a1 ) + h0
(a0 − a1 )
3 3
h2
(a3 − a2 ) + h1
(a1 − a2 )
{C} = ..
.
3 3
hn−1
(an − an−1 ) + hn−2 (an−2 − an−1 )
0
Hence, to construct cubic splines, do the following
1. Make use of Eq. (5.40) and set ai = f (xi ).
2. Solve Eq. (5.51) to obtain all other values of ci ’s. Note that by solving Eq.
(5.51) we are assuming that c0 = 0 and cn = 0. This is equivalent to saying
that the second derivative at x0 and xn are zero !
3. Use Eq. (5.50) to obtain the bi ’s.
4. Use Eq. (5.48) to obtain the di ’s.
After steps 1-4, we can evaluate the function values at the point you are interested
from the following formula Si (x) = ai + bi (x − xi ) + ci (x − xi )2 + di (x − xi )3 . The
only thing tricky is to find the interval that corresponds to your point.
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
5.3. SPLINE INTERPOLATION 101
Exercise 5.10: Given the set of data points below, use cubic splines to find
the value of f at x = 5.
xi f (xi )
3.0 2.5
4.5 1.0
7.0 2.5
9.0 0.5
Exercise 5.11:
You are given the four data points shown in Fig. 5.14 and asked to fit a cubic
spline through the four data points. The cubic spline can be thought as be-
ing made up of three ’real’ quadratic polynomial S0 (x), S1 (x), S2 (x) and one
’imaginary’ polynomial S3 (x). The three polynomials can be written as
S0 (x) = a0 + b0 (x − x0 ) + c0 (x − x0 )2 + d0 (x − x0 )3 (5.53)
S1 (x) = a1 + b1 (x − x1 ) + c1 (x − x1 )2 + d1 (x − x1 )3 (5.54)
S2 (x) = a2 + b2 (x − x2 ) + c2 (x − x2 )2 + d2 (x − x2 )3 (5.55)
S3 (x) = a3 + b3 (x − x3 ) + c3 (x − x3 )2 + d3 (x − x3 )3 (5.56)
(a) Show that
ai = f (xi ) for i =0, 1, 2, 3 (5.57)
(b) By requiring that adjacent Si (x) to have common values at the interior
points, show that
a0 + b0 h0 + c0 h20 + d0 h30 = a1 (5.58)
a1 + b1 h1 + c1 h21 + d1 h31 = a2 (5.59)
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
102 CHAPTER 5. INTERPOLATION
a2 + b2 h2 + c2 h22 + d2 h32 = a3 (5.60)
where hi = xi+1 − xi .
(c) By requiring that the derivative of the adjacent Si (x) to have common values
at the interior points, show that
b0 + 2c0 h0 + 3d0 h20 = b1 (5.61)
b1 + 2c1 h1 + 3d1 h21 = b2 (5.62)
(d) By requiring that the double derivative of the adjacent Si (x) to have common
values at the interior points, show that
c1 − c0
d0 = (5.63)
3h0
c2 − c1
d1 = (5.64)
3h1
c3 − c2
d2 = (5.65)
3h2
(e) Substitute Eqs. (5.63) - (5.65) into Eqs. (5.61) - (5.62) and show that
b1 = b0 + (c0 + c1 ) h0 (5.66)
b2 = b1 + (c1 + c2 ) h1 (5.67)
(f ) Substitute Eqs. (5.63) - (5.65) into Eqs. (5.58) - (5.60) and show that
a1 − a0 h0
b0 = − (2c0 + c1 ) (5.68)
h0 3
a2 − a1 h1
b1 = − (2c1 + c2 ) (5.69)
h1 3
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
5.3. SPLINE INTERPOLATION 103
a3 − a2 h2
b2 = − (2c2 + c3 ) (5.70)
h2 3
(g) Use Eqs. (5.66) - (5.70) together with the assumption that c0 = 0 and c3 = 0
and show that the ci ’s could be obtained by solving the following matrix
equation
1 0 0 0 c0
h0 2(h0 + h1 ) h1 0 c1
=
0 h1 2(h1 + h2 ) h2 c2
0 0 0 1 c3
0
3 (a2 − a1 ) + 3 (a0 − a1 )
h31 h0 (5.71)
(a3 − a2 ) + 3 (a1 − a2 )
h2 h1
0
(f ) Explain how you can obtain the coefficients for S0 (x), S1 (x), S2 (x) from the
results in parts (a)-(g).
Exercise 5.12:
Find the cubic spline that interpolates the following data points
xi f (xi )
1.0 2.0
3.0 4.0
4.0 9.0
7.0 2.0
Exercise 5.13: A natural cubic spline, S(x), which interpolates the data
given below
xi f (xi )
1 1
2 1
3 0
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
104 CHAPTER 5. INTERPOLATION
is defined to be
S0 (x) = 1 + P (x − 1) − Q(x − 1)3 if α ≤ x < β
S(x) =
S1 (x) = 1 + p(x − 2) − 43 (x − 2)2 + q(x − 2)3 if γ ≤ x < λ
Find the values of α, β, γ, λ, P , Q, p and q.
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
Chapter 6
Differentiation and Integration
In order to find a formula to differentiate a function, consider the Taylor series
1 ′′ 1 ′′′
f (xi ) (xi+1 − xi )2 + f (ξ1 ) (xi+1 − xi )3
′
f (xi+1 ) = f (xi ) + f (xi ) (xi+1 − xi ) +
2! 3!
(6.1)
where
xi+1 = xi + ∆
′df
f (x) =
dx
and ξ1 is somewhere in between xi and xi+1 . Equation (6.1) can be re-arranged to
obtain
f (xi+1 ) − f (xi ) 1 ′′ 1 ′′′
− f (xi ) (xi+1 − xi ) + f (ξ1 ) (xi+1 − xi )2
′
f (xi ) = (6.2)
(xi+1 − xi ) 2! 3!
′
Hence f (x) can be approximated as
′ f (xi+1 ) − f (xi )
f (xi ) = (6.3)
(xi+1 − xi )
Equation (6.3) is called the Forward Difference Scheme (FDS) The other terms
that have been neglected in Eq. (6.2) gives an approximation of the error in approx-
′
imating f (xi ) with Eq. (6.3). the leading term in the truncation error for the FDS
approximation is
105
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
106 CHAPTER 6. DIFFERENTIATION AND INTEGRATION
1 ′′
f (ξ1 ) (xi+1 − xi )
EF DS =
2!
The Taylor series expansion Eq. (6.1)could also be used to obtain an expression
for f (xi−1 ).
′
f (xi−1 ) = f (xi ) + f (xi ) (xi−1 − xi )
1 ′′ 1 ′′′
+ f (xi ) (xi−1 − xi )2 + f (ξ2 ) (xi−1 − xi )3
2! 3!
′
= f (xi ) − f (xi ) (xi − xi−1 )
1 ′′ 1 ′′′
+ f (xi ) (xi − xi−1 )2 − f (ξ2 ) (xi − xi−1 )3 (6.4)
2! 3!
′
Equation (6.4) can be re-arranged to obtain an another approximation for f (xi )
′ f (xi ) − f (xi−1 )
f (xi ) = (6.5)
(xi − xi−1 )
Equation (6.5) is called the Backward Difference Scheme (BDS) approximation
of the first derivative. If we assume that xi−1 is very close to xi , then the leading
term in the truncation error for the FDS approximation is
1 ′′
EBDS = f (ξ2 ) (xi − xi−1 )
2!
Exercise 6.1: Subtract Eq. (6.4) from Eq. (6.1) and show that, after some
algebra, that the first derivative can be written as
′ f (xi+1 )−f (xi−1 ) ′′ (xi+1 −xi )2 −(xi −xi−1 )2
f (xi ) = xi+1 −xi−1
− f (x i ) 2!(xi+1 −xi−1 )
′′′ ′′′
f (ξ1 ) (xi+1 −xi )3 f (ξ2 ) (xi −xi−1 )3
− 3! (xi+1 −xi−1 )
− 3! (xi+1 −xi−1 )
′
Following Ex. (6.1), yet another formula for computing f (xi ) is
′ f (xi+1 ) − f (xi−1 )
f (xi ) = (6.6)
xi+1 − xi−1
This formula is called the Central Difference Scheme (CDS) and its leading order
error is given by
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
6.1. APPROXIMATION OF THE 2ND DERIVATIVE 107
′′ (xi+1 − xi )2 − (xi − xi−1 )2
ECDS = f (xi )
2! (xi+1 − xi−1 )
Exercise 6.2:
For the case where all the xi ’s are equally spaced, i.e. xi+1 − xi = ∆ = const,
show that Eqs. (6.3), (6.5), and (6.6) simplifies to the following expressions for
the FDS, BDS and CDS
′ f (xi+1 ) − f (xi )
f (xi ) =
∆
′ f (xi ) − f (xi−1 )
f (xi ) =
∆
′ f (xi+1 ) − f (xi−1 )
f (xi ) =
2∆
Also show that the leading error term in the FDS and BDS is O(∆) and the
leading error term in the CDS is O(∆2 ).
Exercise 6.3: An alternative way of obtaining approximations to the deriva-
tive of a function is to fit a Lagrange polynomial through the data points and
then differentiating the resulting curve.
1. Fit a Lagrange polynomial through two data points and show that you can
obtain the forward and backward differencing scheme.
2. Fit a Lagrange polynomial through three data points and show that you
can obtain the central differencing scheme.
6.1 Approximation of the 2nd Derivative
′
A formula for the 2nd derivative can be obtained by evaluating f (x) at points
halfway between xi+1 and xi and between xi and xi−1 (see Fig. 6.1).
df df
d2 f dx
(xi+1/2 ) − dx (xi−1/2 )
(x i ) = (6.7)
dx 2 1
xi + 2 (xi+1 − xi ) − xi−1 + 21 (xi − xi−1 )
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
108 CHAPTER 6. DIFFERENTIATION AND INTEGRATION
l l m m
xi-1 xi xi+1
xi-1/2 xi+1/2
Figure 6.1: Figure showing the data points used to calculate the second derivative
Using Eq. (6.6) to approximate the first derivative, Eq. (6.7) can be written as
f (xi+1 )−f (xi )
d2 f xi+1 −xi
− f (xxi )−f (xi−1 )
i −xi−1
(xi ) = 1 (6.8)
dx2 2
(xi+1 − xi−1 )
Exercise 6.4: Show that Eq. (6.8) simplifies to
d2 f f (xi+1 ) (xi − xi−1 ) − f (xi ) (xi+1 − xi−1 ) + f (xi−1 ) (xi+1 − xi )
2
(xi ) = 1
dx 2
(xi+1 − xi−1 ) (xi+1 − xi ) (xi − xi−1 )
In addition, if one assumes that xi+1 − xi = xi − xi−1 = ∆, then
d2 f f (xi+1 ) − 2f (xi ) + f (xi−1 )
2
(xi ) = (6.9)
dx ∆2
Exercise 6.5: Apply the forward difference scheme (FDS) twice and show
that an approximation for the 2nd derivative can be written as
d2 f f (xi+2 ) − 2f (xi+1 ) + f (xi )
2
(xi ) = (6.10)
dx ∆2
Similarly, apply the backward difference scheme (BDS) twice and show that an
approximation for the 2nd derivative can be written as
d2 f f (xi ) − 2f (xi−1 ) + f (xi−2 )
2
(xi ) = (6.11)
dx ∆2
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
6.2. INTEGRATION 109
f(x1)
f(x0)
x0 x1
Figure 6.2: Area under a curve
f1(x)
f(x1) f(x1)
f(x0) f(x0)
x0 x0 x1
x1
Area A Area A'
Figure 6.3: The trapezoidal approximation of area under a curve.
6.2 Integration
Integration means finding the area under a curve between two points. Let’s call these
two points x0 and x1 (see Fig. 6.2). If we only have information at these two points,
then the best we can do is to fit a straight line through these two points and find
the area under it. We are estimating A with A’ (see Fig. 6.3). We know that the
Lagrangian interpolating polynomial for two points can be written as
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
110 CHAPTER 6. DIFFERENTIATION AND INTEGRATION
(x − x1 ) (x − x0 ) 1
f1 (x) = f (x0 ) + f (x1 ) + f (2) (ξ(x)) (x − x0 )(x − x1 )
(x0 − x1 ) (x1 − x0 ) 2
So to obtain the area A’ all we have to do is to integrate f1 (x) from x0 to x1 .
Z x1 Z x1
f (x0 )
f1 (x)dx = (x − x1 ) dx
x0 (x0 − x1 ) x0
Z x1
f (x1 )
+ (x − x0 )dx
(x1 − x0 ) x0
Z x1
1 (2)
+ f (ξ1 ) (x − x0 )(x − x1 )dx (6.12)
2 x0
where ξ1 is some number in (x0 , x1 ). Note that the Weighted Mean Value Theorem
for Integrals (see page 8) have been used in evaluating the integral involving the error
term. This can be done because (x − x0 )(x − x1 ) does not change sign if x0 ≤ x ≤ x1 .
Exercise 6.6: Integrate Eq. (6.12) and show that
Z x1
x1 − x0 ∆3 (2)
f1 (x)dx = (f (x0 ) + f (x1 )) − f (ξ) (6.13)
x0 2 12
where ∆ = x1 − x0 . You might like to use the following result
A+m∆
Z
3 m3 m2
(x − (A + k∆)) (x − (A + l∆)) dx = ∆ − (k + l) + klm
3 2
A
The Trapezoidal rule is derived from Eq. (6.13) and it approximates the integral
of f (x) as
Z x1
∆
f (x)dx ≈ (f (x0 ) + f (x1 )) (6.14)
x0 2
In general, at any two points, xi and xi+1 , the trapezoidal rule can be written as
Z xi+1
∆
f (x)dx ≈ (f (xi ) + f (xi+1 )) (6.15)
xi 2
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
6.2. INTEGRATION 111
The error of this approximation is
∆3 (2)
f (ξ)
12
a b a b
Area estimated by
Actual area under a curve the trapezoidal rule
Figure 6.4: Integration over larger interval.
So if the size of the interval, h, is large, the error would be large as well. What
happens if we want to minimize the error as we integrate over a large interval ?
Consider the situation shown in Fig. 6.4. Clearly, in this situation, one cannot
expect to get an accurate answer if one were to use the trapezoidal rule to integrate
over the interval a ≤ x ≤ b. In order to get a more accurate answer, it is common
to divide the domain into smaller intervals, usually of equal length, and apply the
Trapezoidal rule in each of the subintervals (see Fig. 6.5).
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
112 CHAPTER 6. DIFFERENTIATION AND INTEGRATION
∆1 ∆1 ∆2 ∆2 ∆2 ∆2
a=x0 x1 b=x2 a=x0 x1 x2 x3 b=x4
Figure 6.5: Decomposing a large interval into smaller intervals.
Intuitively, it should be clear that as you divide the large interval up into smaller in-
tervals, you should get more accurate answers than simply just applying Trapezoidal
rule once over the entire large interval. What is the general formula for Trapezoidal
rule if we apply it to smaller subintervals ? Consider the situation shown in Fig. 6.6.
If we apply the Trapezoidal rule to each of the subintervals,
∆
A0 = (f (x0 ) + f (x1 ))
2
∆
A1 = (f (x1 ) + f (x2 ))
2
∆
A2 = (f (x2 ) + f (x3 ))
2
∆
A3 = (f (x3 ) + f (x4 ))
2
The trapezoidal rule estimate of the area under the curve would be just the sum of
all the Ai ’s
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
6.3. SIMPSON’S RULE 113
A = A0 + A1 + A2 + A3
∆ ∆ ∆ ∆
= (f (x0 ) + f (x1 )) + (f (x1 ) + f (x2 )) + (f (x2 ) + f (x3 )) + (f (x3 ) + f (x4 ))
2 2 2 2
∆
= (f (x0 ) + 2f (x1 ) + 2f (x2 ) + 2f (x3 ) + f (x4 ))
2 !
3
∆ X
= f (x0 ) + 2 f (xi ) + f (x4 )
2 i=1
In general, the formula for Trapezoidal rule applied to a large interval a < x < b
which is divided into N subintervals of equal size h is
Zb N −1
!
∆ X
f (x)dx ≈ f (a) + 2 f (xi ) + f (b) (6.16)
2 i=1
a
Exercise 6.7: Do the following integral
Z10
1
dx
x
2
using the Trapezoidal rule with increasing number of integrals. Compare with
the exact value of 1.6094379.
6.3 Simpson’s rule
The Trapezoidal rule was derived by fitting the Lagrange polynomial through two
points. A more accurate method of calculating integrals would be try and fit a
Lagrange polynomial through three equally spaced points. Then find the area under
the Lagrange polynomial (see Fig. 6.7).
What we are trying to do here is to estimate the area A under the actual function
with A′ which is the area under the Lagrange polynomial (see Fig. 6.8). For three
data points, the Lagrange interpolating Polynomial can be written as
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
114 CHAPTER 6. DIFFERENTIATION AND INTEGRATION
f(x)
A0 A1 A2 A3
∆ ∆ ∆ ∆
a=x0 x1 x2 x3 b=x4
Figure 6.6: Decomposing a large interval into smaller intervals. The total area under
the curve can be quite accurately represented by A0 + A1 + A2 + A3 .
(x − x1 ) (x − x2 )
f2 (x) = f (x0 )
(x0 − x1 ) (x0 − x2 )
(x − x0 ) (x − x2 )
+ f (x1 )
(x1 − x0 ) (x1 − x2 )
(x − x0 ) (x − x1 )
+ f (x2 )
(x2 − x0 ) (x2 − x1 )
1
+ (x − x0 )(x − x1 )(x − x2 )f (3) (ξ(x))dx
6
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
6.3. SIMPSON’S RULE 115
Actual function
f(x2)
Approximation of
the actual function
f(x1) with Lagrange
polynomial
f(x0) through 3 data
points.
∆ ∆
x0 x1 x2
Figure 6.7: Fitting a Lagrange polynomial through three data points.
where the last term is the error term associated with the Lagrange polynomial.
Exercise 6.8: Assume that x0 , x1 and x2 are equally spaced. Integrate
between x0 and x2 and show that
Zx2
∆
f2 (x) = (f (x0 ) + 4f (x1 ) + f (x2 )) (6.17)
3
x0
where
x1 − x0 = ∆
x2 − x1 = ∆
Ignore the error term. You might like to use the following result
A+m∆
Z
3 m3 m2
(x − (A + k∆)) (x − (A + l∆)) dx = ∆ − (k + l) + klm
3 2
A
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
116 CHAPTER 6. DIFFERENTIATION AND INTEGRATION
f(x2) f(x2)
f(x1) f(x1)
f(x0) f(x0)
∆ ∆ ∆ ∆
x0 x1 x2 x0 x1 x2
Actual area A Calculated area A’
Figure 6.8: Physical illustration of Simpsons rule.
In general the simpson’s rule can be written as
xi+2
Z
∆
f2 (x) = (f (xi ) + 4f (xi+1 ) + f (xi+2 )) (6.18)
3
xi
It is more accurate than the Trapezoidal rule. If you want to integrate a function
over a large interval, it is better to split the function into smaller interval and use
the Simpsons rule over smaller sub-intervals. See the example in the section on
Trapezoidal rule.
For example, if we divide up the domain of a function into 10 intervals, i.e. if we
have 11 data points, (x0 , f (x0 )), (x1 , f (x1 )), (x2 , f (x2 )),....,(x10 , f (x10 )),
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
6.3. SIMPSON’S RULE 117
Z x10
∆
f (x)dx ≈ (f (x0 ) + 4f (x1 ) + f (x2 ))
x0 3
∆
+ (f (x2 ) + 4f (x3 ) + f (x4 ))
3
∆
+ (f (x4 ) + 4f (x5 ) + f (x6 ))
3
∆
+ (f (x6 ) + 4f (x7 ) + f (x8 ))
3
∆
+ (f (x8 ) + 4f (x9 ) + f (x10 ))
3
∆
= [f (x0 )
3
+ 4 (f (x1 ) + f (x3 ) + f (x5 ) + f (x7 ) + f (x9 ))
+ 2 (f (x2 ) + f (x4 ) + f (x6 ) + f (x8 ))
+ f (x10 )]
In general, the simpson’s rule can be written as
N −1 N −2
∆ X X
f (x0 ) + 4 f (xi ) + 2 f (xi ) + f (xN ) (6.19)
3 i=1 i=1
iodd ieven
where N is the number of intervals. Note that Simpson’s rule only works of N is
even i.e. the total number of data points (N + 1) is odd.
Exercise 6.9: A large interval, x0 < x < xN , can be divided into N equally
spaced subintervals of size ∆ with N + 1 data points (see Figure 6.9). Assuming
that N is even, use Eq. (6.17) on two subintervals of size ∆ and show that the
integral of f (x) over the large interval can be approximated as
Z xN N −1
∆ 4∆ X
f (x)dx ≈ f (x0 ) + f (xn )
x0 3 3
n odd
N −2
2∆ X ∆
+ f (xn ) + f (xN ) (6.20)
3 n even 3
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
118 CHAPTER 6. DIFFERENTIATION AND INTEGRATION
f(xN)
f(xN-1)
f(x2) f(x)
f(x1)
f(x0)
∆ ∆ ∆
x0 x1 x2 xN-1 xN
Figure 6.9: Equally spaced subintervals of size ∆ which are inside the large interval,
x0 ≤ x ≤ xN .
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
Chapter 7
Boundary Value Problem
In many engineering problems, you will have to solve ordinary differential equations
that look something like
d2 y dy
2
+ p(x) + q(x)y = r(x) (7.1)
dx dx
where
a≤x≤b
We are also usually given the conditions at the boundaries
y(a) = α
y(b) = β
Exercise 7.1: Using CDS approximation for all the derivatives, show that
Eq. (7.1) can be approximated as
− xi+1p−x
i
i−1
+ 2
(xi+1 −xi−1 )(xi −xi−1 )
yi−1
2
+ − (xi+1 −xi )(x i −xi−1 )
+ qi y i (7.2)
+ xi+1p−x i
i−1
+ (xi+1 −xi−12)(xi+1 −xi ) yi+1 = ri
where
yi = y(xi )
pi = p(xi )
qi = q(xi )
ri = r(xi )
119
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
120 CHAPTER 7. BOUNDARY VALUE PROBLEM
From the boundary conditions, we know that
y0 = α
(7.3)
yN = β
Eq. (7.2) together with Eq. (7.3) represent a linear system that can be solved.
They can be represented by
[A]{X} = {C}
where
1 0 0 0 ··· 0
χ 1 ε 1 η1 0 ··· 0
...
0 χ 2 ε2 η2 0
[A] = ..
0 0 ... ... ...
.
0 0 · · · χN −2 εN −2 ηN −1
0 0 0 0 ··· 1
y0
y1
..
.
{X} = ..
.
..
.
yN
α
r1
r2
{C} = ..
.
rN −1
β
and
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
121
pi 2
χi = − +
x − xi−1 (xi+1 − xi−1 ) (xi − xi−1 )
i+1
2
εi = − + qi
(xi+1 − xi ) (xi − xi−1 )
pi 2
ηi = +
xi+1 − xi−1 (xi+1 − xi−1 ) (xi+1 − xi )
Example 7.1: Suppose you are asked to solve
d2 y dy
2
+ 3 + 2y = 0 (7.4)
dx dx
in the domain 0 ≤ x ≤ 1. You are also given the boundary condition y(0) = 1.0 and
y(1) = e−2 . From the previous chapter, you know that the derivative and double
derivative can be approximated as
dy yi+1 − yi−1
(xi ) = (7.5)
dx 2∆
d2 y yi+1 − 2yi + yi−1
2
(xi ) = (7.6)
dx ∆2
If we now substitute Eqs. (7.5) and Eqs. (7.6) into Eq. (7.4), you will end up with
1 3 2 1 3
− yi−1 + − 2 + 2 yi + + yi+1 = 0 (7.7)
∆2 2∆ ∆ ∆2 2∆
Equation (7.7) is the discretised form of Eq. (7.4). Suppose we want to implement
Eq. (7.7) on a grid as shown in Fig. (7.1). This grid is equally spaced with the
spacing between grid points, x∆ = 0.25. The corresponding values of y on the grid
shown in Fig. (7.1) are y0 , y1 , y2 , y3 , y4 respectively. The yi ’s are the unknowns that
you need to solve for. Thus for this grid, you have 5 unknowns. So implementing
Eq. (7.7) on x = x1 , x = x2 , x = x3 , will give
10y0 − 30y1 + 22y2 = 0 (7.8)
10y1 − 30y2 + 22y3 = 0 (7.9)
10y2 − 30y3 + 22y4 = 0 (7.10)
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
122 CHAPTER 7. BOUNDARY VALUE PROBLEM
We also know from the given boundary conditions that
y0 = 1.0 (7.11)
y4 = e−2 (7.12)
You can put Eqs. (7.8)-(7.12) in matrix form as
1 0 0 0 0 y0
1
10 −30 22 0 0 y1
0
0 10 −30 22 0 y2 = 0 (7.13)
0 0 10 −30 22 y3 0
−2
0 0 0 0 1 y4 e
One can now solve Eq. (7.13) to get
y0
1
y
1
0.5972
y2 = 0.3598 (7.14)
y
0.2192
3
−2
y4 e
which is the approximate solution to Eq. (7.4) at x = 0, 0.25, 0.5, 0.75, 1.00.
∆=0.25
x0=0 x1=0.25 x2=0.5 x3=0.75 x4=1.0
Figure 7.1: Figure showing the grid used in Example 7.1
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
Chapter 8
Ordinary Differential Equations
In many engineering problems, you will need to solve differential equations that look
something like
dx
= f (t, x) (8.1)
dt
in the domain
a≤t≤b
with the initial condition
x(t = a) = α
8.1 Euler’s method
Euler’s method is probably the simplest method used to solve Eq. (8.1). Consider
Taylor’s theorem
dx (tn+1 − tn )2 d2 x
x(tn+1 ) = x(tn ) + (tn+1 − tn ) (tn ) + (ξn ) (8.2)
dt 2 dt2
where ξn is somewhere in between tn and tn+1 . If we let tn+1 − tn = h and substitute
Eq. (8.1) into Eq. (8.2), we get
h2 d2 x
x(tn+1 ) = x(tn ) + hf (tn , xn ) + (ξn ).
2 dt2
If we assume that h is small, then we can neglect the second order term in the
equation above. Thus, we get the formula for Euler’s method
123
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
124 CHAPTER 8. ORDINARY DIFFERENTIAL EQUATIONS
xn+1 = xn + hf (tn , xn ) (8.3)
where xn is the numerical approximation of the exact solution x(tn ). Euler’s method
is the simplest method but it is not very accurate.
Matlab program 8.1 shows how Eulers method can be used to solve
dx
= −2x
dt
with the initial condition x(0) = 1.0 for 0 ≤ t ≤ 3. The matlab program compares
the numerical solution with the analytical solution, which is x(t) = exp(−2t). The
output from this program for h = 0.3 is shown in Fig. 8.1
0.9
0.8
0.7
0.6
0.5
x
0.4
0.3
0.2
0.1
0
0 0.5 1 1.5 2 2.5 3
t
Figure 8.1: The output from Matlab program 8.1
Exercise 8.1: Using Euler’s method, solve
dx
= e−t
dt
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
8.1. EULER’S METHOD 125
Matlab Program 8.1 A matlab program illustrating the Eulers method.
clear all
% Sample program to do Euler integration
% Parameters for simulation
N=10;
h=3/N;
% Initial condition
x(1)=1;
t(1)=0;
%Loop for calculating solution with Eulers method
for n=1:N
dxdt=-2*x(n);
x(n+1)=x(n)+dxdt*h;
t(n+1)=t(n)+h;
end;
%Specifying exact solution
t_plot=linspace(0,10);
x_exact=exp(-2*t_plot);
%Plotting numerical solution with analytical solution
plot(t_plot,x_exact,’k-’,t,x,’ko’);
xlabel(’t’);
ylabel(’x’);
axis([0 3 0 1]);
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
126 CHAPTER 8. ORDINARY DIFFERENTIAL EQUATIONS
for 0 < t < 10. Use x(t = 0) = 0 and
(a) h = 2.0
(b) h = 1.0
(c) h = 0.5
(d) h = 0.1
The solution to Exercise 8.1 is shown in Figure 8.2. Note that, as expected, the
numerical solution gets closer to the exact solution for smaller values of h.
Exact solution x(t)=1-e-t
Approximate solution using Euler’s method
3
3
2.5
2.5
h=1
x(t) 2
h=2 2
1.5
1.5
1
1
0.5
0.5
0
0 2 4 6 8 10 0
0 2 4 6 8 10
t
3 3
2.5
h=0.5 2.5
h=0.1
2 2
1.5 1.5
1 1
0.5 0.5
0 0
0 2 4 6 8 10 0 2 4 6 8 10
Figure 8.2: Figure showing the solution to Exercise 8.1.
8.2 Backward Euler Method
In many applications (see later) the Euler method described in section 8.1 is very
unstable. A more stable method is the backward Euler scheme given by the following
formula:
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
8.3. TRAPEZOIDAL OR CRANK NICOLSON METHOD 127
xn+1 = xn + hf (tn+1 , xn+1 ) (8.4)
Note that the Backward Euler method is an implicit method. xn+1 exist on both
sides of Eq. (8.4). If f is a simple function, then it may be possible to obtain an
explicit expression for xn+1 from Eq. (8.4). If f is a complicated function, then
one may need to use root finding methods such as the Newton Raphson method in
Section 2.3.2 to solve for xn+1 .
8.3 Trapezoidal or Crank Nicolson method
The solution to Eq. (8.1) can also be obtain by integration
Z tn+1
x(t) = x(tn ) + f (t, x(t))
tn
The integration can be approximated using the trapezoidal rule (see 2nd year, Com-
putational mechanics lecture notes)
Z tn+1
h
f (t, x(t)) = (f (tn , xn ) + f (tn+1 , xn+1 )) + O(h3 )
tn 2
where h = tn+1 − tn = h. Hence,
h
xn+1 = xn + (f (tn , xn ) + f (tn+1 , xn+1 )) (8.5)
2
Equation (8.5) is the Crank-Nicolson method. It is sometimes also called the Trape-
zoidal method. This formula is an implicit formula for xn+1 . If f is a simple function,
then we can usually find an explicit expression for xn+1 . If f is a more complex func-
tion, then xn+1 can only be found by using some kind of numerical scheme e.g. the
Newton-Raphson method.
As will be shown later, Crank-Nicolson method and the Backward Euler scheme
are usually more difficult to implement than the explicit Euler scheme. However, as
will be shown later, these two schemes have the added advantage of being a lot more
stable than the explicit Euler scheme.
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
128 CHAPTER 8. ORDINARY DIFFERENTIAL EQUATIONS
Exercise 8.2: Solve
dx
+ 2x = 0
dt
for 0 < t < 3 with x(t = 0) = 1 using
(a) Euler method
(b) Backward Euler method
(c) Crank-Nicolson method
Example 8.1: Consider the nonlinear first order ordinary differential equation
dx
+ ex = 0 (8.6)
dt
with the initial condition x(0) = 1.0. The exact solution to Eq. (8.6) is
1
x(t) = − log t + (8.7)
e
We will now show how the Crank Nicolson method can be used to obtain an
approximate solution to Eq. (8.6). Applying Eq. (8.5) to Eq. (8.6) gives
h h
xn+1 + exn+1 = xn − exn
2 2
One can rearrange the above equation to obtain
h h
xn+1 + exn+1 − xn + exn = 0. (8.8)
2 2
xn is known from the previous time step. One then needs to find xn+1 . It is not
easy to find an explicit expression for xn+1 . The root finding methods from Chapter
2 must be used to find xn+1 numerically. If we used the Newton-Raphson method,
the iterative scheme could be written as
(k) (k)
(k+1) (k) xn+1 + h2 exn+1 − xn + h2 exn
xn+1 = xn+1 − (k)
1 + h2 exn+1
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
8.4. RUNGE KUTTA METHODS 129
(k)
where xn+1 is the k’th guess for the value of xn+1 that satisfies Eq. (8.8). While one
(k) (0)
may use any value for xn+1 , it would be logical to use xn+1 = xn .
Figure 8.3 shows the numerical solution to Eq. (8.6) obtained using both the
Crank-Nicolson and the explicit Euler’s method. The numerical data was obtained
using h = 0.2. As is evident in the figure, the Crank Nicolson method produces a
more accurate solution than the Eulers method.
0.5
0
x
−0.5
−1
0 0.2 0.4 0.6 0.8 1
t
Figure 8.3: Solution to Example 8.1. exact solution (Eq. (8.7)), ◦ Crank-
Nicholson solution, solution using Eulers method.
8.4 Runge Kutta methods
These methods are probably the most popular methods in solving initial value prob-
lems. However, many variations of the Runge-Kutta methods exist. We will derive
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
130 CHAPTER 8. ORDINARY DIFFERENTIAL EQUATIONS
one set below and see why the various formulae are not unique. Let’s start with the
problem we would like to solve
dx
= f (t, x)
dt
In general, the Runge-Kutta schemes can be written as
xn+1 = xn + φ(xn , tn , h)h (8.9)
where h is the interval size, i.e. h = ∆t = ti+1 − ti . φ is known as the incremental
function and it can be interpreted as the slope which is used to predict the new value
of x. In general, φ can be written as
φ = a1 k1 + a2 k2 + a3 k3 + a4 k4 + · · · aN kN (8.10)
where
k1 = f (tn , xn )
k2 = f (tn + p1 h, xn + q11 k1 h)
k3 = f (tn + p2 h, xn + q21 k1 h + q22 k2 h)
k4 = f (tn + p3 h, xn + q31 k1 h + q32 k2 h + q33 k3 h)
.. .. ..
. . .
.. .. ..
. . .
.. .. ..
. . .
kN = f (tn + pN −1 h, xi + qN −1,1 k1 h + qN −1,2 k2 h + · · · + qN −1,N −1 kN −1 h)
For N = 1, we get the first order Runge-Kutta scheme. This is just the same as the
Euler integration scheme presented earlier.
8.4.1 Second Order Runge Kutta Method
If we put N = 2 into Eq. (8.10), we get
φ = a1 k1 + a2 k2
Substituting the above into Eq. (8.9) gives
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
8.4. RUNGE KUTTA METHODS 131
xn+1 = xn + (a1 k1 + a2 k2 ) h
= xn + a1 f (tn , xn )h + a2 f (tn + p1 h, xn + q11 k1 h)h (8.11)
The last term in the above equation is a 2 variable function and the Taylor series
expansion (to the lineariazed approximation) for a two variable function is given by
∂f ∂f
f (t + h, x + ∆) = f (t, x) + h +∆
∂t ∂x
Using this relationship,
∂f ∂f
f (tn + p1 h, xn + q11 k1 h) = f (tn , xn ) + p1 h + q11 k1 h (8.12)
∂t ∂x
We know that
k1 = f (tn , xn )
Hence Eq. (8.12) can be written as
∂f ∂f
f (tn + p1 h, xn + q11 k1 h) = f (tn , xn ) + p1 h + q11 f (tn , xn )h (8.13)
∂t ∂x
Substituting Eq. (8.13) into Eq. (8.11) gives
∂f 2 ∂f
xn+1 = xn + (a1 + a2 )f (tn , xn )h + a2 p1 h + a2 q11 f (tn , xn ) h2 (8.14)
∂t ∂x
We can also write a Taylor series expansion for x in terms of t as
dx d2 x h2
x(tn+1 ) = x(tn ) +
(tn ) h + 2 (tn ) + HOT
dt dt 2!
where HOT stands for higher order terms. Let’s assume that they are small. So the
above equation becomes
dx d 2 x h2
xn+1 = xn + h + 2 (8.15)
dt dt 2!
The problem we are trying to solve has
dx
= f (t, x)
dt
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
132 CHAPTER 8. ORDINARY DIFFERENTIAL EQUATIONS
Substituting the above into Eq. (8.15) gives
df (tn , xn ) h2
xn+1 = xn + f (tn , xn )h +
dt 2
∂f ∂f
df = dt + dx
∂t ∂x
df ∂f ∂f dx
= +
dt ∂t ∂x dt
(8.16)
Hence
1 ∂f 2 1 ∂f
xn+1 = xn + f (tn , xn )h + h + f (tn , xn )h2 (8.17)
2 ∂t 2 ∂x
Comparing Eqs. (8.17) with Eq. (8.14) will give you the following three equations
a1 + a2 = 1
1
a2 p1 =
2
1
a2 q11 =
2
We have four unknowns (a1 , a2 , p1 and q11 ) but only the above three equations. So
there cannot be a unique solution. You have an infinite number of solutions. One
possible solution is to set
1
a2 =
2
then
1
a1 =
2
p1 = 1
q11 = 1
Hence, one possible second order Runge-Kutta (RK-2) time stepping scheme is
1 1
xn+1 = xn + k1 + k2 h (8.18)
2 2
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
8.4. RUNGE KUTTA METHODS 133
where
k1 = f (tn , xn )
k2 = f (tn + h, xn + hk1 )
8.4.2 4th Order Runge Kutta Scheme (RK-4)
The steps described in the previous section could be extended to derive a method
that is 4th order accurate. The 4th order Runge Kutta scheme is by far, the most
popular numerical method for solving ODE. The formula for this scheme can be
written as
1 1 1
xn+1 = xn + k1 + (k2 + k3 ) + k4 h (8.19)
6 3 6
where
k1 = f (tn , xn )
h h
k 2 = f t n + , xn + k 1
2 2
h h
k 3 = f t n + , xn + k 2
2 2
k4 = f (tn + h, xn + hk3 )
A Matlab implementation using Eq. (8.19) to solve the sample problem
dx
= −2x
dt
with initial conditions x(0) = 1 is shown in Matlab program 8.2. The output from
this program with h = 0.3 is shown in Fig. 8.4. It is evident that there is very good
agreement with the analytical solution. By comparing this figure with Fig. 8.1, it
is clear that the 4th order Runge-Kutta method is much more accurate than the
Euler’s method.
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
134 CHAPTER 8. ORDINARY DIFFERENTIAL EQUATIONS
Matlab Program 8.2 A matlab program illustrating the 4th order Runge-Kutta
method.
clear all
% Sample program to do Rungge-Kutta
% integration
% Parameters for simulation
N=10;
h=3/N;
% Initial condition
x(1)=1;
t(1)=0;
%Loop for calculating solution with RK-4 method
for n=1:N
k1=-2*x(n);
k2=-2*(x(n)+(h/2)*k1);
k3=-2*(x(n)+(h/2)*k2);
k4=-2*(x(n)+h*k3);
x(n+1)=x(n)+((1/6)*k1+(1/3)*(k2+k3)+(1/6)*k4)*h;
t(n+1)=t(n)+h;
end;
%Specifying exact solution
t_plot=linspace(0,10);
x_exact=exp(-2*t_plot);
%Plotting numerical solution with analytical solution
plot(t_plot,x_exact,’k-’,t,x,’ko’);
xlabel(’t’);
ylabel(’x’);
axis([0 3 0 1]);
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
8.5. STABILITY AND ERROR ANALYSIS 135
0.9
0.8
0.7
0.6
0.5
x
0.4
0.3
0.2
0.1
0
0 0.5 1 1.5 2 2.5 3
t
Figure 8.4: The output from Matlab program 8.2
8.5 Stability and error analysis
When discussing stability of the numerical methods, one usually considers the model
problem
dx
= λx (8.20)
dt
λ is a constant which can be a complex number. In most engineering problems, the
real part of λ is usually negative. This means that the solution that you are after
will typically decay with time. Applying the Euler method (with timestep ∆t = h)
to Eq. (8.20) gives
xn+1 = xn + λhxn (8.21)
where xn is the numerical solution of x(tn ) where tn is the time tn = nh.
To determine the region of stability, simplify Eq. (8.21) to give
xn+1 = (1 + λh) xn
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
136 CHAPTER 8. ORDINARY DIFFERENTIAL EQUATIONS
Thus the error at any time step n can be written as
xn = x0 (1 + λh)n
= x0 (1 + (λR + iλI )h)n
= x0 σ n
(8.22)
where λR and λI are the real and imaginary parts of λ = λR + iλI .
σ = (1 + hλR + ihλI )
is usually called the amplification factor. If |σ| < 1 then the error will decay with
time. The opposite is true if |σ| > 1. Hence, in order to ensure the stability of the
numerical method, |σ| < 1.
|σ|2 = (1 + hλR )2 + (hλI )2 < 1
This is just a circle of radius 1 centred on (-1,0). This plot is called the stability
diagram and is shown in Fig. 8.5
λ Ih
λ Rh
2
Figure 8.5: The stability diagram of Euler method.
If λ is real and negative, then in order for the numerical method to be stable,
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
8.5. STABILITY AND ERROR ANALYSIS 137
2
h≤ (8.23)
|λ|
Exercise 8.3: Solve
dx
+ 2x = 0 (8.24)
dt
with
x(t = 0) = 1
Use h = 1.2, 0.8 and 0.3 in the domain 0 ≤ t ≤ 30
Exercise 8.4: Perform a stability analysis on the Crank-Nicholson method
and show that it is unconditionally stable, i.e. the method is stable for all values
of h.
Exercise 8.5: Perform stability analysis for the second order Runge-Kutta
method on the model equation
dx
= λx
dt
and show that the region of stability can be obtained by solving the following
equation
λ2 h 2
1 + λh + − eiθ = 0 (8.25)
2
Show also that for the 4th order Runge-Kutta method, the stability region is
obtained by solving
λ2 h 2 λ 3 h 3 λ 4 h 4
λh + + + + 1 − eiθ = 0 (8.26)
2 6 24
The solutions to Eqs. (8.25) and (8.26) are shown in Fig. 8.6.
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
138 CHAPTER 8. ORDINARY DIFFERENTIAL EQUATIONS
8.6 Systems of Ordinary Differential Equations
In many engineering problems, it is essential to solve, not just one, but a set of
ordinary differential equations. In general, the set of equations that you might be
required to solve can be expressed as
dx1
= f1 (x1 , x2 , ...., xN )
dt
dx2
= f2 (x1 , x2 , ...., xN )
dt
dx3
= f3 (x1 , x2 , ...., xN )
dt
dx4
= f4 (x1 , x2 , ...., xN )
dt
. = .
. = .
. = .
dxN
= fN (x1 , x2 , ...., xN )
dt
If you have a linear system, then the above set of equations can be expressed in
matrix-vector notation as
d
{x} = [A] {x} (8.27)
dt
where the curly brackets ({}) denotes a vector and the square brackets [] represents
a matrix. {x} is the vector of dependent variables and [A] is usually a matrix of
constant coefficients. Equation (8.27) can be solved using the methods described
above. Applying the explicit Euler method
{x}n+1 − {x}n
= [A] {x}n
h
{x}n+1 = {x}n + h [A] {x}n
At every time step, one would need to perform a matrix-vector multiplication in
order to obtain the values of the vector {x} at the next time step.
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
8.6. SYSTEMS OF ORDINARY DIFFERENTIAL EQUATIONS 139
If an implicit scheme is required to solve the problem, then one would need to solve
an system of linear algebraic equations at every time step. For example, applying
the trapezoidal/Crank Nicholson scheme to Eq. (8.27) gives
{x}n+1 − {x}n 1 1
= [A] {x}n + [A] {x}n+1
h 2 2
h h
{x}n+1 − [A] {x}n+1 = {x}n + [A] {x}n
2 2
h n+1 h
I − [A] {x} = I + [A] {x}n (8.28)
2 2
Equation (8.28) can now be solved using methods for solving a system of differ-
ential (e.g. LU decomposition, Gauss-Seidel etc.)
Exercise 8.6: Solve the following second order ordinary differential equation
d2 x
+ ω2x = 0 (8.29)
dt2
with the initial conditions
x (t = 0) = p
dx
(t = 0) = 0.
dt
Outline the numerical algorithm for the Eulers and Trapezoidal method.
8.6.1 Multi-dimensional Runge-Kutta methods
In many situations, you might be required to solve a multi-dimensional system using
RK-4 and RK-2. In the examples below, it will be assumed that the system you are
interested in is autonomous, i.e. the right hand side is not a function of time. The
multi-dimensional form or RK-2 to solve the system of 2 equations
dx1
= f1 (t, x1 , x2 )
dt
dx2
= f2 (t, x1 , x2 )
dt
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
140 CHAPTER 8. ORDINARY DIFFERENTIAL EQUATIONS
(i.e. extending Eq. (8.18) to a set of 2 equations) can be written as
1 1
xn+1,1 = xn,1 + k11 + k21 h (8.30)
2 2
1 1
xn+1,2 = xn,2 + k12 + k22 h (8.31)
2 2
k11 , k12 , k21 and k22 can be obtained from
k11 = f1 (tn , xn,1 , xn,2 )
k12 = f2 (tn , xn,1 , xn,2 )
k21 = f1 (tn + h, xn,1 + hk11 , xn,2 + hk12 )
k22 = f2 (tn + h, xn,1 + hk11 , xn,2 + hk12 )
Further extending the RK-2 method to solve the system of 4 equations
dx1
= f1 (t, x1 , x2 , x3 , x4 )
dt
dx2
= f2 (t, x1 , x2 , x3 , x4 )
dt
dx3
= f3 (t, x1 , x2 , x3 , x4 )
dt
dx4
= f4 (t, x1 , x2 , x3 , x4 )
dt
then the solution can be written as
1 1
xn+1,1 = xn,1 + k11 + k21 h (8.32)
2 2
1 1
xn+1,2 = xn,2 + k12 + k22 h (8.33)
2 2
1 1
xn+1,3 = xn,3 + k13 + k23 h (8.34)
2 2
1 1
xn+1,4 = xn,4 + k14 + k24 h (8.35)
2 2
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
8.6. SYSTEMS OF ORDINARY DIFFERENTIAL EQUATIONS 141
where k11 , k12 , k13 , k14 and k21 , k22 , k23 and k24 can be obtained from
k11 = f1 (tn , xn,1 , xn,2 , xn,3 , xn,4 )
k12 = f2 (tn , xn,1 , xn,2 , xn,3 , xn,4 )
k13 = f3 (tn , xn,1 , xn,2 , xn,3 , xn,4 )
k14 = f4 (tn , xn,1 , xn,2 , xn,3 , xn,4 )
k21 = f1 (tn + h, xn,1 + hk11 , xn,2 + hk12 , xn,3 + hk13 , xn,4 + hk14 )
k22 = f2 (tn + h, xn,1 + hk11 , xn,2 + hk12 , xn,3 + hk13 , xn,4 + hk14 )
k23 = f3 (tn + h, xn,1 + hk11 , xn,2 + hk12 , xn,3 + hk13 , xn,4 + hk14 )
k24 = f4 (tn + h, xn,1 + hk11 , xn,2 + hk12 , xn,3 + hk13 , xn,4 + hk14 )
The RK-4 method can also be easily extended to a multi-dimensional system.
Suppose you are asked to solve a system of 4 ordinary differential equations that
looks like
dx1
= f1 (t, x1 , x2 , x3 , x4 )
dt
dx2
= f2 (t, x1 , x2 , x3 , x4 )
dt
dx3
= f3 (t, x1 , x2 , x3 , x4 )
dt
dx4
= f4 (t, x1 , x2 , x3 , x4 )
dt
Equation (8.19) can be written in multi-dimensional form as
1 1 1
xn+1,1 = xn,1 + k11 + (k21 + k31 ) + k41 h (8.36)
6 3 6
1 1 1
xn+1,2 = xn,2 + k12 + (k22 + k32 ) + k42 h (8.37)
6 3 6
1 1 1
xn+1,3 = xn,3 + k13 + (k23 + k33 ) + k43 h (8.38)
6 3 6
1 1 1
xn+1,4 = xn,4 + k14 + (k24 + k34 ) + k44 h (8.39)
6 3 6
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
142 CHAPTER 8. ORDINARY DIFFERENTIAL EQUATIONS
where the approximated solution to xi at time n + 1 is denoted as xn+1,i . The first
subscript refers to the time level and the second subscript of the variable x correspond
to the number of dependent variables. In order to approximate xi at various times,
kij need to be calculated as follows
k11 = f1 (tn , xn,1 , xn,2 , xn,3 , xn,4 )
k12 = f2 (tn , xn,1 , xn,2 , xn,3 , xn,4 )
k13 = f3 (tn , xn,1 , xn,2 , xn,3 , xn,4 )
k14 = f4 (tn , xn,1 , xn,2 , xn,3 , xn,4 )
h h h h h
k21 = f1 tn + , xn,1 + k11 , xn,2 + k12 , xn,3 + k13 , xn,4 + k14
2 2 2 2 2
h h h h h
k22 = f2 tn + , xn,1 + k11 , xn,2 + k12 , xn,3 + k13 , xn,4 + k14
2 2 2 2 2
h h h h h
k23 = f3 tn + , xn,1 + k11 , xn,2 + k12 , xn,3 + k13 , xn,4 + k14
2 2 2 2 2
h h h h h
k24 = f4 tn + , xn,1 + k11 , xn,2 + k12 , xn,3 + k13 , xn,4 + k14
2 2 2 2 2
h h h h h
k31 = f1 tn + , xn,1 + k21 , xn,2 + k22 , xn,3 + k23 , xn,4 + k24
2 2 2 2 2
h h h h h
k32 = f2 tn + , xn,1 + k21 , xn,2 + k22 , xn,3 + k23 , xn,4 + k24
2 2 2 2 2
h h h h h
k33 = f3 tn + , xn,1 + k21 , xn,2 + k22 , xn,3 + k23 , xn,4 + k24
2 2 2 2 2
h h h h h
k34 = f4 tn + , xn,1 + k21 , xn,2 + k22 , xn,3 + k23 , xn,4 + k24
2 2 2 2 2
and
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
8.6. SYSTEMS OF ORDINARY DIFFERENTIAL EQUATIONS 143
h h h h h
k41 = f1 tn + , xn,1 + k31 , xn,2 + k32 , xn,3 + k33 , xn,4 + k34
2 2 2 2 2
h h h h h
k42 = f2 tn + , xn,1 + k31 , xn,2 + k32 , xn,3 + k33 , xn,4 + k34
2 2 2 2 2
h h h h h
k43 = f3 tn + , xn,1 + k31 , xn,2 + k32 , xn,3 + k33 , xn,4 + k34
2 2 2 2 2
h h h h h
k44 = f4 tn + , xn,1 + k31 , xn,2 + k32 , xn,3 + k33 , xn,4 + k34
2 2 2 2 2
Generalization of RK-4 to a system of N dependent variables should be straight
forward.
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
144 CHAPTER 8. ORDINARY DIFFERENTIAL EQUATIONS
Runge−Kutta
−1
−2
−3
−5 −4 −3 −2 −1 0 1 2
Figure 8.6: Figure showing the stability region of the 2nd and 4th order Runge-Kutta
methods.
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
Solutions to selected exercises
145
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
146 Solution to Exercise 2.7
10.0
7.5
5.0
2.5
0.0
-3 -2 -1 0 1 2 3 4
x
-2.5
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
Solution to Exercise 2.7 147
4.0
3.5
3.0
2.5
2.0
2.0 2.5 3.0 3.5 4.0
x
g_1(x)
x
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
148 Solution to Exercise 2.7
0
-7 -6 -5 -4 -3 -2 -1 0 1 2 3 4
x
-2
-4
-6
g_2(x)
x
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
Solution to Exercise 2.7 149
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
150 Solution to Exercise 2.7
40
30
20
10
0
2 4 6 8 10
x
g_3(x)
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
Solution to Exercise 2.10 151
Maple worksheet to illustrate fixed point iteration scheme
Define the function f (x )
2
f := x/x K 2 x K 3 (NewtonRaphsonExample01.1)
Plot f(x) to see how it looks like
10.0
7.5
5.0
2.5
0.0
-3 -2 -1 0 1 2 3 4
x
-2.5
From the graph above, it is easy to see that the roots of the function f (x ) is at x=-1 and x=3.
Let's now see one can use the Newton Raphson. First we need to calculate the derivative
of f (x ).
dfdx := x/2 x K 2 (NewtonRaphsonExample01.2)
The Newton Raphson method is given by the following formula
f 0 xi 1
xi C 1 = xi K
df
dx i
0x 1
Using this formula to iterate with the initial value
x := 4 (NewtonRaphsonExample01.3)
x := 3.166666667
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
152 Solution to Exercise 2.10
x := 3.006410256
x := 3.000010240
x := 3.000000000
x := 3.000000000
x := 3.000000000 (NewtonRaphsonExample01.4)
As you can clearly see, starting from x=4 will converge to the closest root at x=3.
Plot graph close to the root to see if we can find out what is going on.
10.0
7.5
5.0
2.5
0.0
2.0 2.5 3.0 3.5 4.0 4.5 5.0
x
-2.5
Try to find root close to x=-2
x := K2 (NewtonRaphsonExample01.5)
x := K1.166666667
x := K1.006410256
x := K1.000010240
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
Solution to Exercise 2.10 153
x := K1.000000000
x := K1.000000000
x := K1.000000000 (NewtonRaphsonExample01.6)
As you can clearly see, starting from x=-2 will converge to the closest root at x=-1.
Plot graph close to the root to see if we can find out what is going on.
10.0
7.5
5.0
2.5
0.0
-3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0
x
-2.5
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
154 Solution to Exercise 3.3
I1 1 K3 L
J M
A := J 3 1 4M (1)
J M
K2 0 1N
Doolittle's method for LU decomposition requires that the [L] and [U] matrices be written as below for a
3x3 system
I1 0 0 L
J M
J 0 M
L := J l21 1 M
J M
l
K 31 32l 1 N
Iu u u L
J 11 12 13
M
J M
U := J 0 u22 u23 M (2)
J M
K 0 0 u33 N
Multiply the L and U matrices will give you
I u u12 u13 L
J 11 M
J M
temp := J l21 u11 l21 u12 C u22 l21 u13 C u23 (3)
M
J M
K l31 u11 l31 u12 C l32 u22 l31 u13 C l32 u23 C u33 N
The product of [L] and [U] must be equal to [A]. So equate them and you will get
I u u12 u13 L I L
J 11 M J 1 1 K3 M
J M J M
l u l21 u12 C u22 l21 u13 C u23
J 21 11 M =J 3 1 4 M
J M J M
K l31 u11 l31 u12 C l32 u22 l31 u13 C l32 u23 C u33 N K 2 0 1 N
By comparing both sides of the above equations, you will have 9 equations with 9 unknowns. In order
to solve for the 9 unknowns, you
will need to follow the following sequence.
Going across the first row will and compare with the [A] matrix will give you u11, u12 and u13
u11 := 1
u12 := 1
u13 := K3 (4)
Going down the first column will give you l21 and l31
l21 := 3
l31 := 2 (5)
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
Solution to Exercise 3.3 155
Going across the 2nd row will give you u22 and u
23
u22 := K2
u23 := 13 (6)
Going down the 2nd column will give you l32
l32 := 1 (7)
Finally, comparing elements of the 3rd row and 3rd column of both sides will give ou u33
u33 := K6 (8)
We have now found the [L] and [U] matrices using Doolittle's method
I1 0 0L
J M
L = J3 1 0M
J M
K2 1 1N
I1 1 K3 L
J M
U = J 0 K2 13 M (9)
J M
K0 0 K6 N
Check to see if everything is OK by computing L*U-A. If LU=A, then LU-A=0.
I0 0 0L
J M
LU K A = J 0 0 0M (10)
J M
K0 0 0N
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
156 Solution to Exercise 3.4
I pK3C3 qK2 r L
J M
J 5 9 1 M
J 2 q C p K rK K s M
J 4 4 4 M
J M
b := J 41 61 7 15 5 M (1)
J 8 r C K p K q C s
8 2 2 8 M
J M
J 1 M
J s C 17 K 1 p K 5 q C 13 r M
K 8 8 2 2 8 N
The first column of the inverse of A can be found by solving
I1 1 0 2 L Ia L I1 L
J MJ M J M
J2 5 2 0 M Jb M J0 M
J MJ M = J M
J3 8 3 1 M J c M J0 M
J MJ M J M
K5 0 1 3 N Kd N K0 N
In order to make the original equation look like the above, one would need to set p=2, q=2, r = 2 and
s = 1. Substituting these values into b will give you
I 1L
J M
J 1M
J M
J K7 M (2)
J 2 M
J M
J K1 M
J M
K 2 N
The second column of the inverse of A can be found by solving
I1 1 0 2 L Ia L I0 L
J MJ M J M
J2 5 2 0 M Jb M J1 M
J MJ M = J M
J3 8 3 1 M J c M J0 M
J MJ M J M
K5 0 1 3 N Kd N K0 N
In order to make the original equation look like the above, one would need to set p=1, q=3, r = 2 and
s = 1. Substituting these values into b will give you
I 3L
J M
J 2M
J M
J K15 M (3)
J 2 M
J M
J K5 M
J M
K 2 N
The third column of the inverse of A can be found by solving
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
Solution to Exercise 3.4 157
I1 1 0 2 L Ia L I0 L
J MJ M J M
J2 5 2 0 M Jb M J0 M
J MJ M = J M
J3 8 3 1 M J c M J1 M
J MJ M J M
K5 0 1 3 N Kd N K0 N
In order to make the original equation look like the above, one would need to set p=1, q=2, r = 3 and
s = 1. Substituting these values into b will give you
I K2 L
J M
J K5 M
J M
J 4 M
J M
J 41 M (4)
J 8 M
J M
J 13 M
J M
K 8 N
The fourth column of the inverse of A can be found by solving
I1 1 0 2 L Ia L I0 L
J MJ M J M
J2 5 2 0 M Jb M J0 M
J MJ M = J M
J3 8 3 1 M J c M J0 M
J MJ M J M
K5 0 1 3 N Kd N K1 N
In order to make the original equation look like the above, one would need to set p=1, q=2, r = 2 and
s = 2. Substituting these values into b will give you
I 0L
J M
J K1 M
J M
J 4 M
J M
J 5M (5)
J 8M
J M
J 1M
J M
K 8N
Putting all the columns (Eqs. (2) - (5)) will give the inverse of A to be
I 1 3 K2 0L
J M
J K5 K1 M
J 1 2 M
J 4 4 M
J M
Ainv := J K7 K15 41 5M (6)
J 2 2 8 8M
J M
J K1 K5 13 1 M
J M
K 2 2 8 8N
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
158 Solution to Exercise 3.4
Check to see if everything is OK by multiplying A*Ainv. This should give you the identity matrix.
I1 0 0 0L
J M
J0 1 0 0M
J M (7)
J0 0 1 0M
J M
K0 0 0 1N
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
Solution to Exercise 3.8 159
This worksheet outlines the steps one would take to solve the following system of equation
2 2
xyz K x C y K 1.34 = 0
2
xy K z K 0.09 = 0
x y
e K e C z K 0.41 = 0
Warning, the protected names norm and trace have been redefined and
unprotected
Warning, the name GramSchmidt has been rebound
Define the functions
2 2
u (x, y, z ) = x y z K x C y K 1.34
2
v (x, y, z ) = x y K z K 0.09
x y
w (x, y, z ) = e K e C z K 0.41 (1)
The Jacobian, [J] is defined to be
I vu vu vu L
J M
J vx vy vz M
J M
vv vv vv M
[ J ]= J
J vx vy vz M
J M
J vw vw vw M
J M
K vx vy vz N
For this system [ J ] can be calculated to be
Iy zK 2 x x zC 2 y x y L
J M
J M
[J ] = J y x K2 z M (2)
J M
x y
K e Ke 1 N
Ix L
I1 L
J 0M
J M
J M
The initial guess is J y0 M = J 1 M
J M
J M
z K1 N
K 0N
Set up matrix and iterate. The system to solve at every iteration is of the form
I vu vu vu L
J 0 xi, yi, zi 1 vy 0 xi, yi, zi 1 vz 0 xi, yi, zi 1 M I x K x L I u0x, y, z 1 L
J vx M iC 1 i i i i
J MJ M J M
J vv x , y , z vv vv M J M J M
J vx 0 i i i 1 vy 0 i i i 1 vz 0 i i i 1 M J i C 1
x , y , z x , y , z y K yiM = K J
v 0 x , y ,
i i i M
z 1
J MJ M J M
J vw vw vw M K zi C 1 K zi N K w 0 xi, yi, zi 1 N
J M
K vx 0 i i i 1 vy 0 i i i 1 vz 0 i i i 1 N
x, y, z x, y, z x, y, z
Solving
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
160 Solution to Exercise 3.8
Ix K x L
I K1. 3. 1. L J 1 0 I 0.34 L
M
J M J M
J M
J 1. M y
1. K2. , J 1 K y0 M, "=",
J 0.09 M
J M J M
J M
K 2.7183 K2.7183 1. N K z K z N K K.59 N
1 0
gives
Ix K x L
I K.103738540015165082 L
J 1 0
M
J M
J M
y K y0 , "=", J 0.0951802085692621258 M
J 1 M
J M
J M
z K z K K0.0492791657229514693 N
K 1 0N
Thus
Ix L
J 1M I 0.896261459984834862 L
J M J M
y
J 1 M, "=", J 1.09518020856926212 M
J M
J M
K 0.950720834277048565 N
K z1 N
Solving
Ix K x L
I K.7513 3.0425 0.98158 L J 2 1
M I 0.0106 L
J M J
M
J M
J 1.0952 0.89626 K1.9014 M, J y2 K y1 M, "=", J 0.01229 M
J M M
J M
K 2.4504 K2.9898 J
1. N K z K z N K K0.00132 N
2 1
gives
Ix K x L
J 2 1
M I 0.00598659254954604301 L
J M J M
y
J 2 K y1 M, "=",
J 0.00515167675034208187 M
J M J M
K z2 K z1 N K K0.000587063235234869426 N
Thus
Ix L
J 2M I 0.902248052534380895 L
J M J M
y
J 2 M, "=", J 1.10033188531960424 M
J M J M
K z2 N K 0.950133771041813735 N
Solving
Ix K x L
I K.7591 3.0578 0.99275 L J 3 2 I 0.0001 L
J M M
J M
J 1.1003 0.90225 K1.9003 M, y3 K y2 M, "=",
J
J J K0. M
M
J M J M
J M
K 2.4651 K3.0051 1. N K z K z N K K0.00013 N
3 2
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
Solution to Exercise 3.8 161
gives
Ix K x L
I K0.0000200216157804530296 L
J 3 2
M
J M
J M
y3 K y2 , "=", J 0.0000272899433490323116 M
J M
J M
J M
K 0.00000136429381857177163 N
K z3 K z2 N
Thus
Ix L
I 0.902228030918600488 L
J 3M
J M
J M
y J 1.10035917526295335 M
J 3 M, "=",
J M
J M
K 0.950135135335632319 N
K z3 N
Solving
Ix K x L
I K.7590 3.0580 0.99281 L J 4 3 I K0.0002 L
M
J M J M
J M
J 1.1004 0.90223 K1.9003 M, y4 K y3 , "=", J K0.00004 M
J M
J M J M
J M
K 2.4651 K3.0054 1. N K z K z N K 0.00016 N
4 3
gives
Ix K x L I K0.00000667436618378504944 L
J 4 3
M J M
Jy K y M J K0.0000629366550639356396 M
J 4 3 M, "=",
J M
J M K K0.0000126968430495036220 N
K z4 K z3 N
Thus
Ix L
J 4M I 0.902221356552416753 L
J M J M
y
J 4 M, "=", J 1.10029623860788939 M
J M J M
K z4 N K 0.950122438492582821 N
Solving
I L
I K.7590 3.0578 0.99271 L J x5 K x4 M I 0.0001 L
J M J M
J 1.1003 0.90222 K1.9002 M, J y5 K y4 M, "=", J 0.00002 M
J J M
M M J M
K 2.4651 K3.0051 1. N J z
K 5 K z K K0.00012 N
4N
gives
Ix K x L
J 5 4 I K0.0000104001665921580714 L
M
J J M
M
y5 K y4 , "=", J 0.0000307533753822318448 M
J M
J J M
M
K z5 K z4 N K K0.00000194558097252621844 N
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
162 Solution to Exercise 3.8
Thus
Ix L
I 0.902210956385824602 L
J 5M
J M
J M
y J 1.10032699198327166 M
J 5 M, "=",
J M
J M
K 0.950120492911610270 N
K z5 N
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
Solution to Exercise 5.12 163
Specifying the data points
Data := [ [1, 2 ], [3, 4 ], [4, 9 ], [7, 2 ] ] (1)
First define all the hi
h0 := 2
h1 := 1
h2 := 3 (2)
Then get all your ai
a0 := 2
a1 := 4
a2 := 9
a3 := 2 (3)
Need to solve the tridiagonal system
Ic L
I1 0 0 0L J 0M I 0L
J M J M J M
J2 6 1 0M J c1 M J 12 M
J M, J M, "=", J M (4)
J0 1 8 3M J c2 M J K22 M
J M J M J M
K0 0 0 1N J M K 0N
K c3 N
This gives
I 0L
Ic L J M
J 0M J 118 M
J M J
c
J 1M 47 M
J M
J M, "=", J K144 M
(5)
J c2 M J M
J M J 47 M
J M
K c3 N J M
K 0N
Substituting into the formula for bi gives
K95
b0 :=
141
613
b1 :=
141
535
b2 := (6)
141
Using for formula for di gives
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
164 Solution to Exercise 5.12
59
d0 :=
141
K262
d1 :=
141
16
d2 := (7)
47
The cubic spline will consist of your piecewise
polynomials
O 377 95 59
P K xC (x K 1 )
3
1%x and x%3
P 141 141 141
P
P 425 613 118 262
S(x) = Q K C xC (x K 3 )2 K (x K 3 )3 3%x and x%4 (8)
P 47 141 47 141
P
P 871 535 144 2 16 3
P K C xK (x K 4 ) C (x K 4 ) 4%x and x%7
R 141 141 47 47
Plotting your spline and data points to see that everything is OK.
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
Solution to Exercise 5.12 165
11
10
y 6
1
1 2 3 4 5 6 7
x
Data
Cubic Spline
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
166 Solution to Exercise 5.13
Specifying the data points
Data := [ [1, 1 ], [2, 1 ], [3, 0 ] ] (1)
Defining your piecewise polynomials
3
"So(x)=", 1 C P (x K 1 ) C Q (x K 1 )
3 2 3
"S1(x)=", 1 C p (x K 2 ) C (x K 2 ) C q (x K 2 ) (2)
4
The cubic spline will consist of your picewise polynomials. Need to find P, Q, p and q to completely
determine the cubic spline, S (x ).
O
P 1 C P (K1 C x ) C Q (K1 C x )
3
1%x and x%2
P
S(x) = Q (3)
P 1 C p (x K 2 ) C 3 (x K 2 )2 C q (x K 2 )3 2%x and x%3
P
R 4
Specifying that S0 0 x1 1 = S1 0 x1 1 and that S1 0 x2 1 = 0 gives
1CPCQ = 1 (4)
7
CpCq =0 (5)
4
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
Solution to Exercise 5.13 167
Specifying that S0 ' 0 x1 1 = S1 ' 0 x1 1 gives
PC3 Q = p (6)
Specifying that S0 '' 0 x1 1 = S1 '' 0 x1 1 gives
3
6Q= (7)
2
Using Eqs. (4), (5), (6) and (7) to solve for P, Q, p and q gives
1
4
abc := Q = , P =
4
K1
4
1
,p= ,q=
2
K9
4
5 (8)
1
p := (9)
2
K9
q := (10)
4
K1
P := (11)
4
1
Q := (12)
4
Plotting the spline and data points to make sure that everything is OK.
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
168 Solution to Exercise 5.13
2.0
1.8
1.6
1.4
1.2
y 1.0
0.8
0.6
0.4
0.2
0.0
1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0
x
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
Solution to Exercise 8.6 169
We are required to outline the numerical method to solve
d2 x
2
+ ω 2 x = 0. (8.40)
dt
Let
x1 = x
dx
x2 =
dt
Thus Eq. 8.40 can be written as
dx2
= −ω 2 x1 .
dt
So the original problem can be written in terms of the new variables x1 and x2 as a
system of ODE i.e.
d x1 0 1 x1
= (8.41)
dt x2 −ω 2 0 x2
Applying Eulers method to Eq. (8.41) gives
n+1 n ! n
1 x1 x1 0 1 x1
− =
h x2 x2 −ω 2 0 x2
n+1 n n
x1 1 0 x1 0 h x1
= +
x2 0 1 x2 −ω 2 h 0 x2
n+1 n
x1 1 0 0 h x1
= +
x2 0 1 −ω 2 h 0 x2
n+1 n
x1 1 h x1
= (8.42)
x2 −ω 2 h 1 x2
Equation (8.42) involves matrix-vector multiplication. It is given that at time t = 0,
0
x1 p
=
x2 0
So, using Eulers method, the estimated solution at at time t = nh can be written as
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
170 Solution to Exercise 8.6
n+1 n
x1 1 h p
=
x2 −ω 2 h 1 0
If it is desired to solve Eq. (8.40) using Trapezoidal method, then
n+1 n ! n n+1
1 x1 x1 1 0 1 x1 1 0 1 x1
− = +
h x2 x2 2 −ω 2 0 x2 2 −ω 2 0 x2
n+1 n
1 −h/2 x1 1 h/2 x1
2 = 2
ω h/2 1 x2 −ω h/2 1 x2
n+1 n
x1 1 1 h/2 1 h/2 x1
=
x2 1 + ω 2 h2 /4 −ω 2 h/2 1 −ω 2 h/2 1 x2
n+1 n
x1 1 1 − ω 2 h2 /2 h x1
=
x2 1 + ω 2 h2 /4 −ω 2 h 1 − ω 2 h2 /2 x2
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
lOMoARcPSD|51243766
Bibliography
[1] R. Burden and J. Faires. Numerical Analysis. Thomson, 1997.
[2] S. Chapra and R. Canale. Numerical Methods for Engineers. McGraw-Hill, 2002.
171
Downloaded by Emrys Emrys (wayofa2693@kvegg.com)
You might also like
- Numerical Analysis - I. Jacques and C. JuddNo ratings yetNumerical Analysis - I. Jacques and C. Judd110 pages
- An Intuitive Guide To Numerical Methods Heinold PDFNo ratings yetAn Intuitive Guide To Numerical Methods Heinold PDF121 pages
- Numerical Methods For Graduate School: JP Bersamina October 11,2018No ratings yetNumerical Methods For Graduate School: JP Bersamina October 11,201867 pages
- Numerical Analysis - I. Jacques and C. Judd100% (2)Numerical Analysis - I. Jacques and C. Judd110 pages
- An Intuitive Guide To Numerical Methods Heinold100% (1)An Intuitive Guide To Numerical Methods Heinold100 pages
- PIB - Numerical Analysis - Chua, Moore (2016) 58pgNo ratings yetPIB - Numerical Analysis - Chua, Moore (2016) 58pg58 pages
- Numerical Methods: Radostin Simitev Simon CandelaresiNo ratings yetNumerical Methods: Radostin Simitev Simon Candelaresi127 pages
- LEYKEKHMAN 2019 Numerical Analysis Lecture NotesNo ratings yetLEYKEKHMAN 2019 Numerical Analysis Lecture Notes87 pages
- Lecture Notes For Math-CSE 451: Introduction To Numerical Computation100% (2)Lecture Notes For Math-CSE 451: Introduction To Numerical Computation102 pages
- Introduction To Scientific Computing: Using MatlabNo ratings yetIntroduction To Scientific Computing: Using Matlab8 pages
- 11MMET - 25 - Semester 1 - Exam 1 - Final - SolnsNo ratings yet11MMET - 25 - Semester 1 - Exam 1 - Final - Solns8 pages
- Diagnostic Test 1 - Numerical Reasoning and MathematicsNo ratings yetDiagnostic Test 1 - Numerical Reasoning and Mathematics9 pages
- MUSI20149 Music and The Brain Interview Podcast - Audio Transcript-1No ratings yetMUSI20149 Music and The Brain Interview Podcast - Audio Transcript-18 pages
- Answers: Working With Unfamiliar Problems: Part 1No ratings yetAnswers: Working With Unfamiliar Problems: Part 1129 pages
- MAST10006 Calculus 2 Lecture Slides 2023S1No ratings yetMAST10006 Calculus 2 Lecture Slides 2023S1384 pages
- MAST10006 Calculus 2 Exercise Sheets 2023s1No ratings yetMAST10006 Calculus 2 Exercise Sheets 2023s134 pages
- Finite Difference Buckling Analysis of Non Uniform ColumnsNo ratings yetFinite Difference Buckling Analysis of Non Uniform Columns8 pages
- Chapter 9 Uniform Convergence and Integration100% (1)Chapter 9 Uniform Convergence and Integration27 pages
- Optimality Conditions: Unconstrained Optimization: 1.1 Differentiable ProblemsNo ratings yetOptimality Conditions: Unconstrained Optimization: 1.1 Differentiable Problems10 pages
- Vector Through An Angle of 60 Degrees in The Counterclockwise DirectionNo ratings yetVector Through An Angle of 60 Degrees in The Counterclockwise Direction11 pages
- Math10 - PPT - Week 5 - Day 3 - Intro and Graphs of Polynomial Functions100% (4)Math10 - PPT - Week 5 - Day 3 - Intro and Graphs of Polynomial Functions51 pages
- Maths Class Xii Chapter 12 Linear Programming Practice Paper 13No ratings yetMaths Class Xii Chapter 12 Linear Programming Practice Paper 135 pages
- Problems With Solutions in The Analysis of AlgorithmsNo ratings yetProblems With Solutions in The Analysis of Algorithms271 pages
- Solution of Algebraic & Transcendental Equation: Graeffe's Root Squaring MethodNo ratings yetSolution of Algebraic & Transcendental Equation: Graeffe's Root Squaring Method3 pages
- Greedy Algorithms: Comp 122, Spring 2004No ratings yetGreedy Algorithms: Comp 122, Spring 200429 pages
- (Ebook) Calculus: Graphical, Numerical, Algebraic by Franklin D. Demana Daniel Kennedy Bert K. Waits Ross L. Finney ISBN 9780133178579, 0133178579No ratings yet(Ebook) Calculus: Graphical, Numerical, Algebraic by Franklin D. Demana Daniel Kennedy Bert K. Waits Ross L. Finney ISBN 9780133178579, 013317857958 pages
- 1 CIE IGCSE Additional Mathematics Topical Past Paper FunctionsNo ratings yet1 CIE IGCSE Additional Mathematics Topical Past Paper Functions21 pages
- Maths Methods 2023 Unit 1 Trial Exam 2 SolutionsNo ratings yetMaths Methods 2023 Unit 1 Trial Exam 2 Solutions12 pages
