Numer 2017
Numer 2017
Numer 2017
UECM1693
MATHEMATICS FOR PHYSICS II
1 Preliminaries 3
1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.2 Error Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
2 Numerical Differentiation 5
3 Numerical Integration 7
3.1 The Trapezoidal Rule . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
3.2 Simpsons Rule . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
4 Roots Of Equations 10
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
4.2 Bracketing Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
4.2.1 The Bisection Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
4.2.2 The Method of False-Position . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
4.3 Open Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
4.3.1 Fixed-Point Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
4.3.2 Newton-Raphson Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
4.3.3 Finite Difference Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
1
5.3.4 Shifted Inverse Power Method with Scaling . . . . . . . . . . . . . . . . . . . . . 30
6 Optimization 32
6.1 Direct Search Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
6.1.1 Golden Section Search (Maximum) . . . . . . . . . . . . . . . . . . . . . . . . . 32
6.1.2 Gradient Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
2
Chapter 1
Preliminaries
1.1 Introduction
Numerical methods are methods for solving problems on a computer or a pocket calculator. Such
methods are needed for many real life problems that do not have analytic solutions or in other cases,
the analytic solutions may be practically useless.
3
Definition 1.2.2. (Errors In Numerical Methods)
There are two major sources of errors in numerical computations.
(a) Round-off error occurs when a computer or calculator is used to perform real-number calcula-
tions.
Remark. The error arises because for machine computation, the number must be represented
by a number with a finite number of digits. These errors become important as the number of
computations gets very large. To understand the nature of round-off errors, it is necessary to
learn the ways numbers are stored and arithmetic operations are performed in a computer. The
effect of the round-off errors can be illustrated by the following example.
Example 1. (The effect of round-off error)
x2 19
Let f (x) = .
x 13
0.1112 0.1111
(i) If a 4-digit calculator is used to find f (0.3334), we will obtain 1.
0.3334 0.3333
1 1
(ii) On the other hand, by using the formula f (x) = x + , we will obtain 0.3334 + 0.6667.
3 3
(b) Truncation errors are those that result from using an approximation in place of an exact
mathematical procedure.
Example 2. Recall that for all x,
x2 x3
ex = 1 + x + + + .
2! 3!
In particular, if we let x = 1, then
1 1
e=1+1+ + + .
2! 3!
1 1 1 1 1
+ ++
If we use 1 + 1 + to approximate e, then the truncation error is + + .
2! 3! 10! 11! 12!
Remark. (Truncation Errors)
Many numerical schemes are derived from the Taylor series
h2 00
y(x + h) = y(x) + hy 0 (x) + y (x) + .
2!
If the truncated Taylor series used to approximate y(x) is the order n Taylor polynomial
h2 00 hn
pn (x) = y(x) + hy 0 (x) +
y (x) + + y (n) (x),
2! n!
then the approximation is called an nth order method since it is accurate to the terms of order hn .
The neglected remainder term
hn+1 (n+1) hn+2 (n+2)
y (x) + y (x) +
(n + 1)! (n + 2)!
is called the (local) truncation error (TE).
We say the TE per step is of order hn+1 i.e TE= O(hn+1 ).
4
Chapter 2
Numerical Differentiation
We replace the derivatives of a function f by their corresponding difference quotients based on the
Taylor series:
1 1
(i) f (x + h) = f (x) + hf 0 (x) + h2 f 00 (x) + h3 f 000 (x) +
2! 3!
1 1
(ii) f (x h) = f (x) hf 0 (x) + h2 f 00 (x) h3 f 000 (x) +
2! 3!
f (x + h) f (x) 1
(iii) (i) f 0 (x) = hf 00 (x)
h 2!
f (x + h) f (x)
i.e. f 0 (x) = + O(h)
h
(the forward difference formula for f 0 )
f (x) f (x h)
(iv) (ii) f 0 (x) = + O(h)
h
(the backward difference formula for f 0 )
f (x + h) f (x h)
(v) (i) - (ii) f 0 (x) = + O(h2 )
2h
(the central difference formula for f 0 )
f (x + h) 2f (x) + f (x h)
(vi) (i) + (ii) f 00 (x) = + O(h2 )
h2
(the central difference formula for f 00 )
5
Example 3. Given the following table of data:
x 1.00 1.01 1.02 1.03
f (x) 5 6.01 7.04 8.09
(a) Use forward and backward difference approximations of O(h) and a central difference approxima-
tion of O(h2 ) to estimate f 0 (1.02) using a step size h = 0.01.
(b) Calculate the percentage errors for the approximations in part (a) if the actual value f 0 (1.02) =
104.
Answer .
6
Chapter 3
Numerical Integration
Rb
The ideal way to evaluate a definite integral a f (x) dx is, of course, to find a formula F (x) for an
anti-derivative of f. But some anti-derivatives are difficultor impossible to find. For example, there
2
are no elementary formulas for the anti-derivatives of sinx x , 1 + x4 and ex . When we cannot evaluate
a definite integral with an anti-derivative, we turn to numerical methods such as the Trapezoidal Rule
and Simpsons Rule. The problem of numerical integration is the numerical evaluation of integrals
Z b
I= f (x) dx
a
where a and b are constants and f is a function given analytically or empirically by a table of values.
Geometrically , if f (x) 0 x [a, b], then I is equal to the area under the curve of f between a and
b.
7
Theorem 3.1.1 (Error Bound for Trapezoidal Rule). If f 00 is continuous and |f 00 (x)| M for
3
all x [a, b], then ET (ba)
12n2
M.
Z 2
Example 4. Estimate x2 dx by Trapezoidal Rule with n = 4.
1
Z 2
1
Example 5. How many subdivisions should be used in the Trapezoidal Rule to approximate dx
1 x
with an error less than 104 ?
Answer . Since |f 00 (x)| = | x23 | 2 for x [1, 2], we have from Theorem 3.1.1, ET 6n1 2 .
We choose n so that 6n1 2 < 104 :
q
1 4 104 4
6n2
< 10 n > 6 n > 106 40.82. So we can use any n 41. In particular, n = 41 will
2
y0 = Ah2 Bh + C (3.1)
y1 = C (3.2)
2
y2 = Ah + Bh + C (3.3)
Solving these simultaneous equations, we obtain C = y1 , 2Ah2 = y0 +y2 2y1 , and Ap = h3 (y0 +4y1 +y2 ).
An approximation
h for the area under the curve y = f (x) from x = a to x = b is i
h
S = 3 (f (x0 ) + 4f (x1 ) + f (x2 )) + (f (x2 ) + 4f (x3 ) + f (x4 )) + + (f (xn2 + 4f (xn1 ) + f (xn )) .
Z b
hh i
f (x) dx f (x0 ) + 4f (x1 ) + 2f (x2 ) + 4f (x3 ) + + 2f (xn2 ) + 4f (xn1 ) + f (xn ) , where
a 3
ba
h= n .
Theorem 3.2.1 (Error Estimate for Simpsons Rule). If f (4) is continuous and |f (4) (x)| M
5
for all x [a, b], then the Simpsons rule error |ES | (ba)
180n4
M.
8
Z 1
Example 6. Approximate 4x3 dx by Simpsons Rule with n = 4.
0
Z 1
2
Example 7. Use Simpsons Rule with n = 10 to approximate the integral ex dx.
0
Estimate the error involved in this approximation.
2
x
Answer h . (i) f (x) = e , h = 0.1, xi = 0.1i. i
0.1
S = 3 f (x0 ) + 4f (x1 ) + 2f (x2 ) + + 2f (x8 ) + 4f (x9 ) + f (x10 )
h i
= 0.1
3
f (0) + 4f (0.1) + 2f (0.2) + + 2f (0.8) + 4f (0.9) + f (1)
h i
= 0.1
3
e 0
+ 4e 0.01
+ 2e 0.04
+ 4e 0.09
+ 2e0.16
+ 4e0.25
+ 2e0.36
+ 4e 0.49
+ 2e0.64
+ 4e0.81
+ e1
1.462681
(ii) 0 x 1
2
0 f (4) (x) = (12 + 48x2 + 16x4 )ex max f (4) (x) = f (4) (1) = 76e
0x1
5
76e(1 0)
ES 0.000115.
180(10)4
Z 4
170
Exercise 1. Estimate dx by using
2 1 + x2
(a) Trapezoidal Rule
(b) Simpsons Rule
with n = 6 subintervals. Calculate the absolute error in each case.
Answer . (a) 413, 0.6043 (b) 400, 13.6043
Answer . 20.7.
9
Chapter 4
Roots Of Equations
4.1 Introduction
Definition 4.1.1. Any number r for which f (r) = 0 is called a solution or a root of that equation
or a zero of f.
b
Example 8. The root of the linear equation ax + b = 0 is x = , a 6= 0.
a
Example 9. The roots of the quadratic equation ax2 + bx + c = 0 are given by
b b2 4ac
x= .
2a
Example 10. Find the roots of the equation x3 4x = 0.
Example 11. Show that the equation x = cos x has at least one solution in the interval (0, /2).
10
Reading Assignment 4.2.1. Show that the equation x3 4x + 1 = 0 has at least one solution in
the interval (1, 2).
We must have a stopping criterion s to terminate the computation. In practice, we require an error
estimate that is not contingent on foreknowledge of the root . In particular, an approximate percentage
a can be calculated as
present approximation previous approximation xnew xold
|a | = 100% = r new r 100%
present approximation xr
11
Example 12. Use bisection to find the root of f (x) = x10 1. Employ initial guesses of xl = 0 and
xu = 1.3 and iterate until the estimated percentage error a falls below a stopping criterion of s = 8%.
Answer .
(c) f (0.65)f (0.975) > 0 the root is in (0.975, 1.3). So set xl = 0.975, xu = 1.3. Hence
0.975 + 1.3
x3r = = 1.1375.
2
The estimated error is 1.1375 0.975
|a | = 100% = 14.3%.
1.1375
(d) f (0.975)f (1.1375) < 0 the root is in (0.975, 1.1375). So set xl = 0.975, xu = 1.1375. Hence
0.975 + 1.1375
x4r = = 1.05625.
2
The estimated error is 1.05625 1.1375
|a | = 100% = 7.7%.
1.05625
After 4 iterations, the estimated percentage error is reduced to less than 8%.
12
4.2.2 The Method of False-Position
In this method, we use the graphical insight that
if f (c) is closer to zero than f (d), then c is closer to the root than d (which is not true in general).
13
Example 13. Use the Method of False-Position to find the zero of f (x) = x ex . Use initial guesses
of 0 and 1.
Answer .
Reading Assignment 4.2.2. Use the Method of False-Position to find the zero of f (x) = x ex .
Use initial guesses of 0 and 1. Iterate until two successive approximations differ by less than 0.01.
Answer . Let ea = |xnew
r xold
r |.
14
4.3 Open Methods
In contrast to the bracketing methods, the open methods are based on formulas that require a single
starting value or two starting values that do not necessarily bracket the root. Hence, they sometimes
diverge from the true root. However, when they converge, they tend to converge much faster than the
bracketing methods.
Examples of open methods:
(a) Fixed-Point Method
(b) Newton-Raphson Method
(c) Secant Method
xn+1 = g(xn ), n = 0, 1, 2, . . .
Remark. This method is also called the successive substitution method , or one-point iteration.
Example 14. The function f (x) = x2 3x + ex 2 is known to have two roots, one negative and one
positive. Find the smaller root by using the fixed-point method.
Answer . The smaller root is the negative root.
f (1) = 2 + e1 > 0 and f (0) = 1 < 0 f (1)f (0) < 0 the negative root lies in (1, 0).
x2 + ex 2
f (x) = 0 can be written as x = = g(x), so the algorithm is
3
x2n + exn 2
xn+1 = g(xn ) = , n = 0, 1, 2, . . .
3
If we choose the initial guess x0 = 0.5, we obtain
x2 + ex0 2 (0.5)2 + e0.5 2
x1 = 0 = 0.3811564468 and so on. The results are summarized in
3 3
the following table :
k xk
0 0.5(initial guess)
1 0.381156446
2 0.390549582
3 0.390262048
4 0.390272019
5 0.390271674
6 0.390271686
7 0.390271686
Hence the root is approximately 0.390271686.
15
Note 1:
(a) There are many ways to change the equation f (x) = 0 to the form x = g(x) and the speed of
convergence of the correspondingiterative sequences {xn } may differ accordingly. For instance,
if we use the arrangement x = 3x ex + 2 = g(x) in the above example, the sequence might
not converge at all.
(b) It is simple to implement but in this case slow to converge.
(c) Even in the case where convergence is possible, divergence can occur if the initial guess is not
sufficiently close to the root.
16
4.3.2 Newton-Raphson Method
This is a method used to approximate a root of an equation f (x) = 0 assuming that f has a continuous
derivatives f 0 . It consists of the following steps:
Step 1. Guess a first approximation x0 to the root.
(A graph may be helpful.)
Step 2. Use the first approximation to get the second, the second to get the third, and so on, using
the formula
f (xn )
xn+1 = xn 0 , n = 0, 1, 2, 3, . . .
f (xn )
where xn is the nth approximation.
Stop when |xn+1 xn | < s , the pre-specified stopping criterion s .
xn+1 xn
Remark. You may also stop when 100% < s %
xn+1
Note 2: The underlying idea is that we approximate the graph of f by suitable tangents.
If you are writing a program for this method, dont forget also to include an upper limit to the number
of iterations in the procedure.
Example 16. Use Newton-Raphson Method to approximate the root of f (x) = x ex = 0 that lies
between 0 and 2. Continue the iterations until two successive approximations differ less than 108 .
(xn exn ) xn + 1
Answer . The iteration is given by xn+1 = xn x
=
1+e n exn + 1
Starting with x0 = 1, we obtain
x0 + 1 1+1
x 1 = x0 = 0.537882842
e +1 e+1
x1 + 1
x 2 = x1 0.566986991
e +1
x3 0.567143286
x4 0.56714329
Since |x4 x3 | = 0.000000004 = 0.4 108 < 108 , we stop the process and take x4 0.56714329 as
the required root.
17
Reading Assignment 4.3.2. The equation 2x3 + x2 x + 1 = 0 has only one real root. Use the
Newton-Raphson Method to find this root. Continue the iterations until two successive approximations
differ less than 0.0001.
Answer . Let f (x) = 2x3 + x2 x + 1.
(a) First, use the Intermediate Value Theorem to local the root.
Since f (2)f (1) = (9)(1) = 9 < 0, the root is in (2, 1).
(b) Using the iterative formula
4
Exercise 3. Use Newton-Raphson method to approximate 22 correct to six decimal places.
Answer . 2.165737
Advantages
(a) It needs only 1 initial guess.
(b) It converges very rapidly .
Disadvantages
(a) It may not converge if the initial guess is not sufficiently close to the true root.
(b) The calculations of f 0 (xn ) may be very complicated.
(c) Difficulties may arise if |f 0 (xn )| is very small near a solution of f (x) = 0.
18
4.3.3 Finite Difference Method
Consider a second order boundary value problem (BVP)
Suppose a = x0 < x1 < < xn1 < xn = b with xi xi1 = h for all i = 1, 2, . . . , n. Let
yi = y(xi ), Pi = P (xi ), Qi = Q(xi ) and fi = f (xi ). Then by replacing y 0 and y 00 with their central
difference approximations in the BVP, we get
yi+1 2yi + yi1 yi+1 yi1
2
+ Pi + Qi yi = fi , i = 1, 2, . . . , n 1.
h 2h
or, after simplifying
h h
1 + Pi yi+1 + (h2 Qi 2)yi + (1 Pi )yi1 = h2 fi .
2 2
The last equation, known as a finite difference equation, is an approximation to the differential
equation. It enables us to approximate the solution at x1 , . . . , xn1 ..
y 00 4y = 0, y(0) = 0, y(1) = 5.
That is,
y2 2.25y1 + y0 = 0
y3 2.25y2 + y1 = 0
y4 2.25y3 + y2 = 0
y2 2.25y1 = 0
y3 2.25y2 + y1 = 0
2.25y3 + y2 = 5
Notes : We can improve the accuracy by using smaller h. But for that we have to pay a price, i.e.
we have to solve a larger system of equations.
19
Chapter 5
That is, A is strictly diagonally dominant if the absolute value of each diagonal entry is greater than
the sum of the absolute values of the remaining entries in the same row.
20
2 7 4
Example 18. Show that A = 8 1 6 is not strictly diagonally dominant.
3 5 9
Answer . Since in the first row, |a11 | = 2 > 6 |a12 | + |a13 | = 7 + 4 = 11 and in the second row,
|a22 | = 1 6> |a21 | + |a23 | = 8 + 6 = 14, matrix A is not strictly diagonally dominant.
8 1 6
However, if we interchange the first and second rows, the resulting matrix B = 2 7 4 is strictly
3 5 9
diagonally dominant since
|a11 | = 8 > |a12 | + |a13 | = 1 + 6 = 7,
|a22 | = 7 > |a21 | + |a23 | = 2 + 4 = 6,
|a33 | = 9 > |a31 | + |a32 | = 3 + 5 = 8.
Theorem 5.1.1 ( Convergence of the Iterative Methods). If the square matrix A is strictly
diagonally dominant, then the Gauss-Seidel and Jacobi approximations to the solution of the linear
system Ax = b both converge to the exact solution for all choices of the initial approximation.
x1 10x2 + x3 = 13
20x1 + x2 x3 = 17
x1 + x2 + 10x3 = 18
Iterate until |xi [p+1] xi [p] | < 0.0002, for all i. Prepare all the computations in 5 decimal places.
Answer .
(i) To ensure the convergence of this method, we rearrange the equations to obtain a strictly diago-
nally dominant system :
20x1 + x2 x3 = 17
x1 10x2 + x3 = 13
x1 + x2 + 10x3 = 18
21
(ii) Make each diagonal element the subject of each equation:
1
x1 = (17 x2 + x3 )
20
1
x2 = (13 + x1 + x3 ) ()
10
1
x3 = (18 + x1 x2 )
10
(iii) Start with a reasonable initial approximation to the solution:
[0] [0] [0]
e.g. x1 = 0, x2 = 0, x3 = 0.
(iv) Substitute this initial approximation into the RHS of (*), and calculate the new approximation
1 [p] [p]
x1 [p+1] = (17 x2 + x3 )
20
1 [p] [p]
x2 [p+1] = (13 + x1 + x3 )
10
1 [p] [p]
x3 [p+1] = (18 + x1 x2 )
10
[p]
where xi is the pth iteration of the approximation to xi .
22
Example 20. Use Gauss-Seidel Method to solve
x1 10x2 + x3 = 13
20x1 + x2 x3 = 17
x1 + x2 + 10x3 = 18
Iterate until |xi [p+1] xi [p] | < 0.0002, for all i. Prepare all the computations in 5 decimal places.
Answer .
(iv) Substitute this initial approximation into the RHS of (*), and calculate the new approximation
1 [p] [p]
x1 [p+1] = (17 x2 + x3 )
20
1 [p+1] [p]
x2 [p+1] = (13 + x1 + x3 )
10
1 [p+1] [p+1]
x3 [p+1] = (18 + x1 x2 )
10
[p]
where xi is the pth iteration of the approximation to xi .
That is, the new approximation is
1 [0] [0]
x1 [1] = (17 x2 + x3 ) = 0.850
20
1 [1] [0]
x2 [1] = (13 + x1 + x3 ) = 1.215
10
1 [1] [1]
x3 [1] = (18 + x1 x2 ) = 2.0065
10
23
(v) To improve the approximation, we can repeat the substitution process. The next approximation
is
1
x1 [2] = (17 (1.215) + 2.0065) = 1.01108
20
1
x2 [2] = (13 + 1.0111 + 2.0065) = 0.99824
10
1
x3 [2] = (18 + 1.0111 (0.99824)) = 2.00093
10
Note 3:
(a) The Gauss-Seidel technique is not appropriate for use on vector computer, as the set of equations
must be solved in series.
(b) The Gauss-Seidel technique requires less storage than the Jacobi technique and leads to a conver-
gent solution almost twice as fast.
4 1 1
Example 21. Let A = . Verify that v1 = is an eigenvector of A and find the corresponding
3 2 1
eigenvalue.
4 1 1 1
Answer . Av1 = = = 5v1 is an eigenvector of A corresponding to the eigenvalue
3 2 1 1
= 5.
Note 4:
24
(b) To find the eigenvectors corresponding to , we solve the linear system
(A I)x = 0
for nonzero x.
4 1
Exercise 4. Find the eigenvalues and eigenvectors of A = .
3 2
1 1
Answer . 1 = 1, 2 = 5; v1 = , v2 =
3 1
Definition 5.2.2. A square matrix A is called diagonalizable if there is an invertible matrix P such
that P 1 AP is a diagonal matrix; the matrix P is said to diagonalize A.
2 0 0 1
Exercise 5. Let A = and P = . Verify that P diagonalizes A.
3 1 1 1
1 1 2 0 0 1 1 0
Answer . P 1 AP = = , a diagonal matrix, so P diagonalizes A.
1 0 3 1 1 1 0 2
Theorem 5.2.1.
Let A be an n n matrix.
(a) A is diagonalizable if and only if it has n linearly independent eigenvectors.
(b) These n linearly independent eigenvectors form a basis of Rn .
25
5.3 Approximation of Eigenvalues
Definition 5.3.1. An eigenvalue of an n n matrix A is called the dominant eigenvalue of A if its
absolute value is larger than the absolute values of the remaining n 1 eigenvalues. An eigenvector
corresponding to the dominant eigenvalue is called a dominant eigenvector of A.
Example 22.
1 = 5, 2 = 4, 3 = 2, 4 = 3
1 = 2, 2 = 5, 3 = 5
26
4 2
Example 23. Let A = . Use the power method to approximate the dominant eigenvalue and
3 1
1
a corresponding eigenvector with x0 = . Prepare 6 iterations.
0
Answer .
Remark.
is less than the pre-specified error criterion . Unfortunately, the actual value is usually unknown.
So instead, we will stop the computation at the ith step if the estimated relative error
(i) (i 1)
< .
(i)
The value obtained by multiplying estimated relative error by 100% is called the estimated
percentage error.
27
5.3.2 Power Method with Scaling
Theorem 5.3.2. The power method often generates a sequence of vectors {Ak x0 } that have incon-
veniently large entries.
This can be avoided by scaling the iterative vector at each step. That is, we
x1
.. 1
multiply Ax0 = . by and label the resulting matrix as x1 . We repeat the
max{|x1 |, |x2 |, . . . , |xn |}
xn
process with the vector Ax1 to obtain the scaled-down vector x2 , and so on.
Example 24. Repeat the iterations of Example 23 using the power method with scaling.
Answer .
4 1 1 4 1
(i) y1 = Ax0 = , we define x1 = y1 = = .
3 4 4 3 0.75
4 2 1 2.5 1 1 2.5 1
(ii) y2 = Ax1 = = , we define x2 = y2 = = .
3 1 0.75 2.25 2.5 2.5 2.25 0.9
Continue in this manner, we will obtain
1 1 1 1
x3 , x4 , x5 , x6 .
0.9545 0.9782 0.9894 0.9947
28
5.3.3 Inverse Power Method with Scaling
Theorem 5.3.3. If in addition, A is an invertible matrix, then
1 > 2 3 n > 0.
1
Hence will be the dominant eigenvalue of A1 . We could apply the power method on A1 to find
n
1
and hence n . This is called the inverse power method.
n
The algorithm for the Inverse Power Method with Scaling as follows:
1. Compute B = A1
2. Compute yk = Bxk1
1
3. Set xk = yk
max{|yk |}
Bxk xk
4. Dominant eigenvalue of matrix A1 = =
xk xk
1
5. Smallest eigenvalue of A =
8 8
Example 25. Let 1 and 2 be the eigenvalues of A = such that 1 > 2 > 0.
1 2
(a) Use an appropriate power method with scalingto approximate an eigenvector corresponding to
1
2 . Start with the initial approximation x0 = . Round off all computations to four decimal
1
places, and stop after two iterations.
(b) Use the result of part (a) to approximate 2 .
(c) Find the estimated percentage error in the approximation of 2 .
Answer .
(a) We use the inverse power method with scaling to find 2 , the smallest eigenvalue of A.
1 0.0833 0.3333
B=A =
0.0417 0.3333
0.4166 1.0000
Iteration 1 : y1 = Bx0 = , x1 = y1 /0.4166 =
0.2916 0.7000
0.3166 1.0000
Iteration 2 : y2 = Bx1 = , x2 = y2 /0.3166 = is an approximation to the required
0.1916 0.6052
eigenvector.
x2 x2
(b) An approximation to 2 is the second approximation 2 (2) = = 3.5782
Bx2 x2
x1 x1
(c) The first approximation of 2 is 2 (1) = = 3.3058
Bx1 x1
2 (2) 2 (1)
The percentage error is 100% = 7.6128%
2 (2)
29
5.3.4 Shifted Inverse Power Method with Scaling
Theorem 5.3.4. Let A be a diagonalizable n n matrix with eigenvalues 1 , . . . , n and assume that
This method can be used to find any eigenvalue and eigenvector of A. Let a be any number and let k
be the eigenvalue of A closest to a. The inverse power iteration with A aI will converge to |k a|1
and a multiple of vk .
The algorithm for the Shifted Inverse Power Method with Scaling as follows:
1. Compute C = (A aI)1
2. Compute yk = Cxk1
1
3. Set xk = yk
max{|yk |}
Cxk xk
4. Dominant eigenvalue of matrix (A aI)1 = =
xk xk
1
5. Eigenvalue closest to a = + a
7 2
Example 26. Apply the shifted inverse power method with scaling (2 steps) to A = to find
3 1
1
the eigenvalue nearest to a = 6. Start with the initial guess x0 = .
0
1 2
Answer . A aI = A 6I = .
3 7
7 2
C = (A aI)1 = = .
3 1
Iteration 1:
7 1 1
y1 = Cx0 = , x1 = y 1 = .
3 7 0.4286
Iteration 2:
6.1428 1 1
y2 = Cx1 = , x2 = y2 = is an approximation to an eigenvector correspond-
2.5714 6.1428 0.4186
ing to the required eigenvalue,say, .
Let be the dominant eigenvalue of C.
Cx2 x2 7.2433 6.1628
Then = = 6.1634 since Cx2 = .
x2 x2 1.1752 2.5814
1 1
Hence, = a + = 6 + = 6.1622.
6.1634
30
Remark.
(a) The Power Method converges rather slowly. Shifting can improve the rate of convergence.
(b) If we have some knowledge of what the eigenvalues of A are, then this method can be used to find
any eigenvalue and eigenvector of A.
(c) We can estimate the eigenvalues of A by using the Gerschgorins Theorem.
That is, the eigenvalues of A lie in the union of the n discs with radius ri centered at aii .
Furthermore, if a union of k of these n discs form a connected region that is disjoint from all the
remaining n k discs, then there are precisely k eigenvalues of A in this region.
D1 : Center 7, radius 8
31
Chapter 6
Optimization
Definition 6.1.1. A function f (x) is unimodal on [a, b] if it has exactly one maximum (or minimum)
on [a, b].
32
Example 28. Find the maximum value of f (x) = 3 + 6x 4x2 on the interval [0, 1] using the golden
section search method with the final interval of uncertainty having a length less than 0.25.
Answer . Solving
rk (b a) < 0.25
for the number k of iterations that must be performed, we obtain
k > 2.88.
Iteration 1 : xL = 0, xR = 1
51
x1 = xR (xR xL ) = 0.381966,
2
51
x2 = xL + (xR xL ) = 0.618034
2
f (x1 ) = 4.708204 < f (x2 ) = 5.180340 take xL = x1 = 0.381966 with the same xR .
Iteration 2 : xL = 0.381966, xR = 1
51
x1 = xR (xR xL ) = 0.618034,
2
51
x2 = xL + (xR xL ) = 0.763932
2
f (x1 ) = 5.180340 < f (x2 ) = 5.249224 take xL = x1 = 0.618034 with the same xR i.e., xR = 1 .
Iteration 3 : xL = 0.618034, xR = 1
51
x1 = xR (xR xL ) = 0.763932,
2
51
x2 = xL + (xR xL ) = 0.854102
2
f (x1 ) = 5.249224 > f (x2 ) = 5.206651 take xL remains the same with the same xR = x2 = 0.854102
i.e., xL = 0.618034, xR = 0.854102
Thus, I3 = [0.618034, 0.854102] and L3 = 0.854102 0.618034 = 0.236068 < 0.25 as wanted.
So the required maximum point lies within I3 = [0.618034, 0.854102]
33
6.1.2 Gradient Method
Review. Recall that the gradient vector of a function f at a point x,
f (x),
points in the direction of maximum increase (or the direction of steepest ascent). And,
f (x),
Theorem 6.1.1 (Method of Steepest Ascent/ Gradient Method). An algorithm for finding
the nearest local maximum of a twice continuous differentiable function f (xk ), which presupposes that
the gradient of the function can be computed. The method of steepest ascent, also called the gradient
method, starts at a point x0 and, as many times as needed, moves from xk to xk+1 by maximizing
along the line extending from xk in the direction of f (xk ), the local uphill gradient. That is, we
determine the value of t and the corresponding point
xk+1 = xk + tf (xk )
Remark : This method has the severe drawback of requiring a great many iterations (hence a slow
convergence) for functions which have long, narrow valley structures.
Example 29. Use the method of steepest ascent to determine a maximum of f (x) = x2 y 2 starting
from the point x0 = (1, 1).
34
Reading Assignment 6.1.2. Starting at the point (0, 0), use one iteration of the steepest-descent
algorithm to approximate the minimum value of the function
Answer .
fx = 4(x + y) + 2(x y) + 3 = 6x + 2y + 3
fy = 4(x + y) 2(x y) + 2 = 2x + 6y + 2
f (x, y) = (fx , fy ) = (6x + 2y + 3, 2x + 6y + 2)
Start with the initial approximation x0 = (0, 0),
we obtain the next approximation
x1 = x0 + tf (x0 ) = (0, 0) + t(3, 2) = (3t, 2t)
since f (x0 ) = f (0, 0) = (3, 2).
Find t that minimizes g(t) = f (x1 ) = f (3t, 2t) = 2(5t)2 + (t)2 + 3(3t) + 2(2t) = 51t2 + 13t
Solve g 0 (t) = 102t + 13 = 0
t = 13/102 = 0.127
Hence x1 = (0.382, 0.255) and f (x1 ) = f (0.382, 0.255) = 0.828 minimum value of f
35
Chapter 7
yn+1 = yn + hf (xn , yn )
xn+1 = xn + h = x0 + (n + 1)h, n = 0, 1, 2, . . . .
(a) As the truncation is performed after the first term of the Taylor series of y(x), Eulers method is
a first order method with TE = O(h2 ). This error is a local error as it applies at each and every
step as the solution develops.
(b) The global error , which applies to the final solution is O(h), since the number of operations
1
would be proportional to .
h
36
Example 30. Consider the IVP
y 0 = x + y, y(0) = 1.
(a) Use Eulers method to obtain a five-decimal approximation to y(0.5) using the step size h = 0.1.
(b) (i) Estimate the truncation error in your approximation to y(0.1) using the next two terms in
the corresponding Taylor series.
(ii) The exact value for y(0.1) is 1.11034 (to 5 decimal places) . Calculate the error between
the actual value y(0.1) and your approximation y1 . How does this error compare with the
truncation error you obtained in (i)?
(c) The exact value for y(0.5) is 1.79744 (to 5 decimal places) . Calculate the absolute error between
the actual value y(0.5) and your approximation y5 .
37
Example 31. Given that the exact solution of the initial-value problem
y 0 = x + y, y(0) = 1
is y(x) = 2ex x 1. The following table shows
(a) the approximate values obtained by using Eulers method with step size h = 0.1
(b) the approximate values with h = 0.05, and the actual values of the solution at the points
x = 0.1, 0.2, 0.3, 0.4, 0.5 :
38
Example 32. Use Improved Eulers method with step size h = 0.1 to approximate the solution of the
IVP
y 0 = x + y, y(0) = 1
on the interval [0, 1]. The exact solution is y(x) = 2ex x 1. Make a table showing the approximate
values , the actual values together with the absolute errors.
39
7.1.3 Taylor Series Method of Order p
h2 00 h3 (3) hp
yn+1 = yn + hyn0 + yn + yn + + yn(p)
2! 3! p!
Example 33.
40
Example 34. Consider y 0 = x + y, y(0) = 1. Use Runge-Kutta method of order 4 to obtain an
approximation to y(0.2) using the step size h = 0.1.
41
7.1.5 Multi-step Methods
The methods of Euler, Heun, and RK are called one-step methods because the approximation for
the mesh point xn+1 involves information from only the previous mesh point, xn . That is, only the
initial point (x0 , y0 ) is used to compute (x1 , y1 ) and in general, yn is needed to compute yn+1 .
Methods which use more than one previous mesh points to find the next approximation are called
multi-step methods.
Example 35. Use the Adams-Bashforth/Adams-Moulton Method with h = 0.2 to obtain an approx-
imation to y(0.8) for the solution of
y 0 = x + y, y(0) = 1.
Answer . With h = 0.2, y(0.8) y4 . We use the RK4 method as the starter method, with x0 =
0, y0 = 1 to obtain
(a) y1 = 1.242800000
(b) y2 = 1.583635920
(c) y3 = 2.044212913
(d) fn = f (xn , yn ) = xn + yn
f0 = x0 + y0 = 0 + 1 = 1
f1 = x1 + y1 = 0.2 + 1.242800000 = 1.442800000
f2 = 0.4 + 1.583635920 = 1.983635920
f3 = 0.6 + 2.044212913 = 2.644212913
0.2 0.2
(e) y4 = y3 + (55f3 59f2 + 37f1 9f0 ) = 2.044212913 + (55(2.644212913) 59(1.983635920) +
24 24
37(1.442800000) 9(1)) = 2.650719504
(f) f4 = f (x4 , y4 ) = x4 + y4 = 0.8 + 2.650719504 = 3.450719504
0.2
(g) y4 = y3 + (9f4 + 19f3 5f2 + f1 ) = 2.651055757 y(0.8)
24
42
7.1.7 Summary:Orders Of Errors For Different Methods
Method order Local TE Global TE
Euler 1 O(h2 ) O(h)
Heun 2 O(h3 ) O(h2 )
Second-Order Taylor Series 2 O(h3 ) O(h2 )
Third-Order Taylor Series 3 O(h4 ) O(h3 )
Fourth-Order Runge-Kutta 4 O(h5 ) O(h4 )
Adams-Bashforth/Adams-Moulton 4 O(h5 ) O(h4 )
43
7.2 Higher-Order Initial Value Problems
The methods use to solve first order IVPs can be applied to higher order IVPs because a higher IVP
can be replaced by a system of first-order IVPs. For example, the second-order initial value problem
can be decomposed to a system of two first-order initial value problems by using the substitution
y0 = u :
y0 = u , y(x0 ) = y0
u = f (x, y, u) , u(x0 ) = y00
0
Each equation can be solved by the numerical techniques presented earlier. For example, the Euler
method for this system would be
yn+1 = yn + hun
un+1 = un + hf (xn , yn , un )
Remark. The Euler method for a general system of two first order differential equations,
u0 = f (x, y, u) , u(x0 ) = u0
y 0 = g(x, y, u) , y(x0 ) = y0
is given as follows
un+1 = un + hf (xn , yn , un )
yn+1 = yn + hg(xn , yn , un )
Reading Assignment 7.2.1. Use the Euler method to approximate y(0.2) and y 0 (0.2) in two steps,
where y(x) is the solution of the IVP
y 00 + xy 0 + y = 0, y(0) = 1, y 0 (0) = 2.
y0 = u
0
u = xu y
yn+1 = yn + 0.1un
un+1 = un + 0.1(xn un yn )
44
7.2.1 The Linear Shooting Method
Consider the linear second order BVP
with boundary conditions y(a) = , y(b) = . The shooting method replaces the BVP by two IVPs
in the following way:
(a) Let y1 be the solution of the IVP
Replace this second-order IVP by a system of two first-order IVPs. Then solve each of the two
first-order IVPs using, say, the fourth-order Runge-Kutta method to obtain y1 .
(b) Let y2 be the solution of the IVP
45
Example 36. Use the shooting method (together with the fourth-order Runge-Kutta method with
1
h = ) to solve the BVP
3
y 00 = 4(y x), 0 x 1, y(0) = 0, y(1) = 2.
Answer .
y0 =
00 0 u , y(0) = 0 [1]
(a) y = 4(y x), y(0) = 0, y (0) = 0 0
u = 4(y x) , u(0) = 0 [2]
(i) Use RK4 to solve [1] and [2] : we obtain
y1 ( 13 ) = 0.02469136, y1 ( 23 ) = 0.21439821 and y1 (1) = 0.80973178
y 0 = u , y(0) = 0 [3]
00 0
(b) y = 4y, y(0) = 0, y (0) = 1
u0 = 4y , u(0) = 1 [4]
(i) Use RK4 to solve [3] and [4] : we obtain
y2 ( 31 ) = 0.35802469, y2 ( 23 ) = 0.88106488 and y2 (1) = 1.80973178
2 y1 (1) 2 + 0.80973178
(c) y(x) = y1 (x) + y2 (x) = y1 (x) + y2 (x)
y2 (1) 1.80973178
y( 31 ) = 0.53116634, y( 23 ) = 1.15351499, y(1) = 2
Note 9:
An advantage of the shooting method is that the existing programs for initial value problems
may be used.
However, the shooting method is sensitive to round-off errors and it becomes rather difficult to
use when there are more than two boundary conditions. For these reasons, we may want to use
alternative methods. This brings us to the next topic.
46
7.2.2 Finite Difference Method
Consider a second order BVP
Suppose a = x0 < x1 < < xn1 < xn = b with xi xi1 = h for all i = 1, 2, . . . , n. Let
yi = y(xi ), Pi = P (xi ), Qi = Q(xi ), and fi = f (xi ). Then by replacing y 0 and y 00 with their central
difference approximations in the BVP, we get
yi+1 2yi + yi1 yi+1 yi1
2
+ Pi + Qi yi = fi , i = 1, 2, . . . , n 1.
h 2h
or, after simplifying
h h
1 + Pi yi+1 + (h2 Qi 2)yi + (1 Pi )yi1 = h2 fi .
2 2
The last equation, known as a finite difference equation, is an approximation to the DE. It enables
us to approximate the solution at x1 , . . . , xn1 ..
Notes : We can improve the accuracy by using smaller h. But for that we have to pay a price, i.e.
we have to solve a larger system of equations.
47
Chapter 8
48
8.2 Numerical Approximation To Derivatives :1-Variable Func-
tions
We replace derivatives by their corresponding difference quotients based on the Taylor series :
1 1
(i) u(x + h) = u(x) + hu0 (x) + h2 u00 (x) + h3 u000 (x) +
2! 3!
1 1
(ii) u(x h) = u(x) hu0 (x) + h2 u00 (x) h3 u000 (x) +
2! 3!
u(x + h) u(x) 1
(iii) (i) u0 (x) = hu00 (x)
h 2!
u(x + h) u(x)
i.e. u0 (x) = + O(h) (the forward difference formula for u0 )
h
u(x) u(x h)
(iv) (ii) u0 (x) = + O(h) (the backward difference formula for u0 )
h
u(x + h) u(x h)
(v) (i) - (ii) u0 (x) = + O(h2 ) (the central difference formula for u0 )
2h
u(x + h) 2u(x) + u(x h)
(vi) (i) + (ii) u00 (x) = + O(h2 )
h2 00
(the central difference formula for u )
Then
u ui+1,j ui,j
= + 0(h) the FD for ux at P
x P h
u ui,j+1 ui,j
= + 0(k)the FD for uy at P
y P k
u ui+1,j ui1,j
= + 0(h2 ) the CD for ux at P
x P 2h
2 u ui+1,j 2ui,j + ui1,j
2
= 2
+ 0(h2 )the FD for uxx at P
x P h
49
8.4 Methods for Parabolic Equations
k
where r = .
h2
Note 10: We use the above formula to estimate the values of u at time level j + 1 using the values at
level j. This is known as the FTCS explicit finite difference method.
u 2u
= , 0 x 1, t0
t x2
subject to the initial condition
u(x, 0) = sin x, 0x1
and boundary conditions
u(0, t) = u(1, t) = 0, t 0.
Using an x step of h = 0.1 and a t step of k = 0.0005, use the FTCS scheme to estimate u(0.1, 0.001).
Answer .
k 0.0005
r= 2
= = 0.05
h (0.1)2
(a) u(0, t) = u(1, t) = 0 u0,j = u10,j = 0
(b) u(x, 0) = sin x ui,0 = sin 0.1i
(c) ui,j+1 = (1 2r)ui,j + r(ui+1,j + ui1,j )
ui,j+1 = 0.9ui,j + 0.05(ui+1,j + ui1,j )
50
(d) j = 0 : ui,1 = 0.9ui,0 + 0.05(ui+1,0 + ui1,0 )
(i) i = 1 : u1,1 = 0.9u1,0 + 0.05(u2,0 + u0,0 ) = 0.3075
since ui,0 = sin 0.1i u0,0 = 0, u1,0 = sin 0.1 = 0.3090, u2,0 = sin 0.2 = 0.5878
(ii) i = 2 : u2,1 = 0.9u2,0 + 0.05(u3,0 + u1,0 ) = 0.5849
since u3,0 = sin 0.3 = 0.8090
(e) j = 1 : ui,2 = 0.9ui,1 + 0.05(ui+1,1 + ui1,1 )
(i) i = 1 : u1,2 = 0.9u1,1 + 0.05(u2,1 + u0,1 ) = 0.3060 u(0.1, 0.001).
Note 11: The FTCS scheme is stable for a given space step h provided the time step k is restricted
by the condition
k 1
0<r= 2 .
h 2
The restriction recommends that we do not move too fast in the t-direction.
Note 12: Suppose there is an initial condition : u(x, 0) = f (x). Then the scheme starts with (i, j) =
(0, 0), the left-hand corner of the 2D grid. Along the horizontal line j = 0, where t = 0, we have
Note 13: Suppose we include the boundary conditions : u(0, t) = u(a, t) = 0. Then we have
k
where r =
h2
Note 14: For each time level j, we obtain an (M 1) (M 1) tridiagonal system which can be
solved using iterative methods.
Note 15: The CNI scheme has no stability restriction. It is more accurate than the explicit scheme.
51
Example 42. Consider the one-dimensional heat equation
u 2u
= , 0 x 1, t0
t x2
subject to the initial condition
u(x, 0) = sin x, 0x1
and boundary conditions
u(0, t) = u(1, t) = 0, t 0.
Using an x step of h = 0.2 and a t step of k = 0.001, use the Crank-Nicolson method to estimate
u(0.4, 0.001).
Answer .
Since the initial temperature distribution is symmetric with respect to x = 0.5, we only need to
consider the grid points over 0 x 0.5.
52
Example 43.
The exact solution of the 1-D heat equation.
ut = uxx
u(0, t) = u(1, t) = 0 , t>0
u(x, 0) = f (x) = sin x , 0 < x < 1
2
is u(x, t) = e t sin(x).
(a) By using the FTCS explicit scheme with h = 0.2 and k = 0.008 , we obtain the following table
(b) By using the FTCS and the CNI scheme with h = 0.2 and k = 0.04 , we obtain the following table
53
8.5 A Numerical Method for Elliptic Equations
Example 44 (Difference Equation for the Laplace Equation). Using C.D for both uxx and uyy ,
the finite difference approximation for the 2-D Laplace equation uxx + uyy = 0 is
Example 45. The four sides of a square plate of side 12 cm made of uniform material are kept at
constant temperature such that
Using a mesh size of 4 cm, calculate the temperature u(x, y) at the internal mesh points A(4, 4), B(8, 4),
C(8, 8) and D(4, 8) by using the Gauss-Seidel method if u(x, y) satisfies the Laplace equation uxx +
uyy = 0, start with the initial guess uA = uB = 80, uC = uD = 50. Continue iterating until all
[p+1] [p] [p]
|ui ui | < 103 where ui are the pth iterates for ui , i = A, B, C, D.
Answer .
54
8.6 CTCS scheme for Hyperbolic Equations
Example 46 (CTCS scheme for Hyperbolic Equations). Consider the BVP involving the 1-D
wave equation
2u 2u
=
t2 x2
u(0, t) = u(1, t) = 0, t>0
u(x, 0) = f (x), 0x1
ut (x, 0) = g(x), 0x1
Use a 2-D grid with x = h and t = k. By replacing both uxx and utt with the central difference
quotients, we obtain a CTCS scheme
ui,j+1 2ui,j + ui,j1 ui+1,j 2ui,j + ui1,j
=
k2 h2
i.e. ui,j+1 = ui+1,j + 2(1 )ui,j + ui1,j ui,j1
where = (k/h)2 .
Note 16:
55
Example 47. If the string (with fixed ends at x = 0 and x = 1) governed by the hyperbolic partial
differential equation
utt = uxx , 0 x 1, t 0
with boundary conditions
u(0, t) = u(1, t) = 0
starts from its equilibrium position with initial velocity
g(x) = sin x
Answer .
56
Bibliography
[1] Koay Hang Leen, Lecture Note for Numerical Methods and Statistics, 2007.
57