CHE 306 NOTES 1 COURSE OUTLINES A.
B.
ITERATIVE METHODS – Solution Solution of Equation of f(x) = 0
1.
Bisection Method
2.
Fixed-Point Iterative (FPI) Method
3.
Jacobi Iteration
4.
Gauss-Seidel Iteration
INTERPOLATION INTERPOLATION AND POLYNOMIAL APPROXIMATION APPROXIMATION
1.
Linear Interpolation
2.
Quadratic Interpolation
3.
Lagrange interpolation
4.
Divided Difference Interpolation – Interpolation – Divided Divided differencesand polynomials
5.
Equispaced Interpolations - Difference operators and difference Tables - Forward, backward and central differences
C.
MUMERICAL INTEGRATION AND DIFFERENTIATION
1.
Numerical Differentiation
2.
Difference Notation and Operators
3.
Numerical Integration - Trapezoidal Rule - Simpson‟s Rule - Mid-Point Rule - Romberg Integration
D.
NUMERICAL SOLUTION OF INITIAL VALUE PROBLEMS
1.
Euler Method
2.
Runge-Kutta Methods
3.
Predictor-Corrector Methods.
INTRODUCTION
In the process of solving problems in Science, Engineering, Economics, etc., a physical situation is first converted into a mathematical model. This is often called call ed formulation of the problem. This mathematical model often gives rise to mathematical problems which are too difficult to solve in a neat closed form e.g. 1
(i.)
Integration: Find
e
x2
dx
0
(ii.) (iii.) (iv.)
cos x x Nonlinear Equation: Solve cos Linear Algebra: Find the eigenvalues of a large matrix. Differential equations: Solve a system of nonlinear differential equations.
When such problem arises, numerical analysis is then used for developing techniques to find a solution or approximate solution of the mathematical equations describing the model. A numerical method (or a combination of numerical methods) which can be to solve a problem is often called an algorithm. An algorithm is algorithm is a complete and unambiguous set of procedures leading to the solution of a mathematical problem. The results obtained for the solution of a problem will be affected by various source of error. Numerical analysts must consider how much accuracy accurac y is required, estimate the magnitude of round-off and discretization errors, determine an appropriate step-size or the number of iterations required, provide for checks on the accuracy and make allowance for corrective action in cases of non-convergence. The efficiency of any numerical method (or algorithm) must also be considered. An algorithm would be of no practical use if it required the largest computer error built to obtain a useful answer. The final phase in solving a problem is programming. Programming is the transformation of the algorithm into a set of unambiguous step-by-step instructions for the computer. In this segment of the course, we will look at the design (formulation) and analysis of various numerical methods and assess them in terms of accuracy, efficiency and computer effort. This will involve some mathematical analysis and a nd some practical work using MATLAB.
Mathematical Mathematical Modelling: Example 1.1 0
Cannon
Target d Figure 1.1: Aiming a Cannon
Where = angle of elevation. v0= muzzle speed. From kinematic and projectile motion, the distance travelled by the canon ball is obtained as follows. x
V
0
cos t
(1.1)
y V0 sin t 1 / 2 gt 2
(1.2)
. V0 cos t x x y V0 sin t 1 / 2 gt 2
(1.3)
x 0 x g y
(1.4)
.
..
..
..
Analysis: When y 0,
V0 sin t 1 / 2 gt 2 0 2
t (V0 sin t 1 / 2 gt )
t
2
0,V0 si sin t 1 / 2 gt 0
V0 sin t 1 / 2 gt 2 0
t
2V 0 sin g
Distance travelled to hit the target, x v0 cos t
(1.5)
v0 cos t .
2v0 sin
g
2
x
2v0 sin cos
g
(1.6) 2
The distance travelled by the canon ball is
2v0 sin cos
, where g= acceleration due to
g
gravity. In order to find the correct elevation to hit the target requires satisfying the equation d
2v0 2 sin cos g
So, we require to solve f ( ) f ( )
(1.7)
0 . Where
2v0 2 sin cos
g
d
(1.8)
Note the following observations: (i)
(ii)
The modelling process gives an idealisation: Some features have been ignored e.g. air resistance, length of the muzzle. They may be significant. The nonlinear equation may not have a solution. The maximum range is
v0
2
g
when (iii)
/
4 .So
if d
v0 2 g
, the target is out of range.
The nonlinear equation may have many solutions. If so is
is a solution, then
* 2k for any integer k. These are trivial rotated solutions. If a
solution / 2 then so is / 2 * . The equation can be rearranged and solved easily as
(iv)
*
2
2v0 sin cos
g
d
v0 2 (2sin cos ) g sin(2 )
1 2
(1.9)
v0 2 (sin 2 ) g
d
(1.10)
dg v0 2
dg 2 v 0
sin 1
(1.11)
Normally, in such problems, a closed form solution is not possible and so an approximate solution is sought from a numerical method.
Example 1.2
A sphere of radius r and density
is
floating in water of unit density. It is required to
determine the depth h of the sphere that is below the waterline. Solution Mass of the sphere = Volume x Density 4
3
3
r
(1.12)
r
h
Displaced water mass: Applying Archimede‟s Principle, a body partially submerged in a fluid is buoyed up by a force equal to the weight of the displaced fluid. Mass of displaced water 1/ 3 h2 (3r h) 1 equating the two weight masses, we have 4 3
r3 3
1 3
h 2 (3r h) 1
2
4r h (3r h) 4
3rh
2
r3
h3 r 3
2
h h 4 3 r r
3
Define x = r 4 3 x 2 x3 x
3
3x
2
4 0
(1.13)
Equation (1.13) is a cubic polynomial equation with three zeros. If 0.6 , to what depth does the sphere sink as a fraction of its radius? This example can be visualized as finding the values of x for which the graph f ( x) = 0 i.e. touches the x-axis.
y 3
2 1
x -1
1
2
3
-1 -2
SOLUTION OF EQUATIONS OF ONE VARIABLE We have the problem of finding values of x that satisfy f ( x)
0 for a given function f . A
solution is called a zero of f or a root of f ( x) 0 . The need for numerical methods arises since equation of the form f ( x) 0 rarely provide closed-form expressions for the roots. A well-known equation having a closed form expression for its roots is the quadratic equation. ax
2
bx c
0
(1.14)
The roots of (2) are defined explicitly by x
b
b 2 4 ac 2a
,a 0
(1.15)
And simply require the substitution of values for a, b and c into (1.15).
LOCATING A ROOT
To develop a numerical method for finding the rots of f ( x) 0 , it is useful to first determine an interval a x b that contains at least one root. So, we need to ask. What properties must a function „ f of x‟ satisfy on the interval a x b to guarantee at least one root in this interval? Example 1.3:
Consider the polynomial curve, x 3 3x 2 2 x (shown below). The curve cuts the x-axis at x 0, x 1 and
solutions.
x
2,
indicating that the cubic equation
x
3
2
3x 2 x
0 has three
y 1.5
1.0 0.5
x -1
1
2
3
-0.5 -1.0
A root of f ( x ) 0 will be denoted by p or x*. So f ( x*) has value 0 or f ( p ) . An important application of iteration is in the problem of finding zero(s) of a function of one variable x, the problem is written as: Find all x: f ( x ) 0 . A root of equation f ( x ) 0 is a value of x which when substituted into the L.H.S of the equation gives zero. The value is then a zero of the function f . We consider several numerical methods for this problem.
1.1 BISECTION METHOD Is a numerical method for solving f ( x) 0 . The underlying mathematics for this methodis
contained Intermediate Value Theorem (IVT). If the function f is continuous on [ a,b] and k is any number lying between f ( a ) and f (b) then there is a point C somewhere in (a,b) such that f (c ) k . y y = f ( x)
f (a) k f (b)
a
c
Figure 2: Initial value Theorem
b
x
For the equation f ( x)
0, we use k
0. then
the IVT tells us that if f is continuous on [ a,b]
and f ( a ) and f (b) have different signs, then there is a solution of f ( x )
0 between a and b.
however, there might be more than one. y
y
y = f ( x)
solution
a
b
solution
y =f ( x)
x a
b solution
solution solution
solution
Figure 3: Solutions using IVT methods So, f (a ), f (b ) opposite signs, then at least one solution exists. If f ( a ) and f (b) have the same sign, then the theorem does not apply, so we cannot say whether is a solution or not. y
y
a
b
x a
b
x
Figure 4: IVT with signs changes PROCEDURES 1. Starting with an interval [a,b] (obtained by lucky or informed guess) with
f (a ) f b 0 then there is at least one root a or b, in the interval [a,b]. Note: i. That if f (a )
0, then a is a root or if f (b)
0, then b is a root.
ii. An interval [a,b] with f (a ), f (b) 0 is called a Nontrivial Bracket for a root of f ( x ) 0 .
2. Suppose [a,b] is a nontrivial bracket for a root of f ( x) 0, then the value at the midpoint c 1 / 2( a b) gives 3 possibilities. (i)
f (c) 0, then c is a root and the problem is solved (unlikely in general).
(ii)
f (c ) 0 and f (c). f (b) 0, then [c, b] is a nontrivial bracket for f ( x) = 0.
(iii)
f (c ) 0 and f (c). f (b) 0, then [a, c] is a nontrivial bracket for f ( x) = 0.
With case (ii) or (iii), the resulting nontrivial bracket is half the size of the original.
x
3. This process is repeated until the interval containing the root is small enough or any convergence is decided upon (used). Example 1: (Use your calculator here) Find a root of x 2 0.5 0 , using the Bisection method starting with the interval [0, 1]. Use 4 steps (i.e. calculate 4 subintervals).
n 0 1 2 3 4
a 0
c 0.5
f (a)
b 1
f (c)
f (b)
Some Notes with using University Calculator 1.
REPLAY FUNCTION: is used to safe retyping
2.
ANS: holds the last number calculated.
3.
STO: stores a number in a register A – F, X, Y, M. e.g. ANS stores the last number calculated in A.
4.
RCL: recalls a number from a register. e.g.
5.
Changing mode from rad to degree etc. in the display.
n 0 1 2 3 4
a 0 0.5 0.5 0.625 0.6875
c 0.5 0.75 0.625 0.6875 0.71875
RCL
Mode
STO
A the number in register A. recalls then
Mode
f (a) -
b 1 1 0.75 0.75 0.75
So
small D
f (c) + -
The root lies in [0.6875, 0.75] a4 ) 1/
note
1
Approximate root c4= 0.71875
The error in the approximation is 1 / 2(b4
A
2(0.75 0.6875)
0.03125
p c4 0.03125
The root cannot be further from C 4 than 1/ 2(b4 a4 ) but it may be closer.
f (b) + + + + +
y y = f ( x)
p x c0
a0
b0
a1
c1
a2
c2 a3 c3
b1
b2 b3
Figure 5: Bisection Method through IVT solution The root is at x = p. Thus the bisection method generates subintervals [an, bn], n = 1, 2, 3, …all containing the root p with [a0, b0] = [a, b]. The corresponding midpoints c0, c1,c2, c3, …are a sequence of approximations for p. BISECTION METHOD ALGORITHM
INPUT:
Nontrivial Bracket - [a, b], Tolerance – ε.
OUTPUT:
Approximate solution, c.
STEP 1:
WHILE
a b
do
STEPS 2-4.
STEP 2:
SET
STEP 3:
IF f (c) = 0,
c
1 2
(a b) .
THEN OUTPUT (c), STOP
STEP 4:
IF f (a) x f (c) < 0
THEN SET
b = c
SET
a = c
ELSE
STEP 5:
1
( a b) , 2 OUTPUT (c), STOP
SET
c
BISECTION METHOD ALGORITHM (MATLAB PROGRAM) function c = bisect(a, b, epsilon) %BISECT
- Bisection method for f(x) = 0.
%
BISECT (A, B, EPSILON) is an approximate solution on the interval [a, b].
%
The final subinterval has length < = EPSILON
NOTES on Bisection MATLAB Program 1.
LINE 1: begins with keyword function followed by output argument c and the “=” symbol. This is followed by the input arguments within round brackets.
2.
LINE 2: called the H1 (help 1) line. It should be a comment line of special form: beginning with %, followed without space by the function name in capitals, followed by one or more spaces and then a brief description. The description should begin with a capital, end with a period “.” and omit the words „the‟ and „a‟.
3.
All comments lines from H1 up to the first non-comment line (usually a blank line for readability) are displayed when “ helpfunction_Name” is typed in the command window after the “ >> “ symbol. These lines should describe the function and its arguments (capitalized by convention).
Note that the function name must be the same as the name of the m-file in which it is stored. So here “ bisect” and “ bisect.m” for the m-file.
SAMPLE OF MATLAB CODES
function c = bisect(a,b,epsilon) %BISECT
Bisection method for f(x) = 0.
%
BISECT(A,B,EPSILON) is an approximate solution on the
%
interval [A,B]. The final subinterval has length <= EPSILON. disp('
a
b
c
fa
fc
fb');
fa = f(a); fb = f(b); % check initial bracket if fa*fb>= 0 error('fa*fb is non-negative! The bracket is out of range!!!') end % bisection loop
SAMPLE OF A FUNCTION m-FILE while abs(b-a) > epsilon c = (a + b)*0.5; fc = f(c); disp([a,b,c,fa,fc,fb]); if fc == 0 return end if fa*fc < 0
% root is to left of mid-point
b = c; fb = fc; else
% root is to right of mid-point a = c; fa = fc;
end end
function y = f(x) %F(X)
Function of X.
%
Y = F(X) computes given function of X.
%
F(X) = x^2 - 0.5. y = x^2 - 0.5;
INTERVAL LENGTH
If the original bracket has length
L0 i.e L0
(calculation of k brackets after the original) is
b
Lk
a
, then the bracket length after k steps L0
2
k
. This formula allows us to calculate
the number of steps required to make the bracket length Lk E . Thus, L0
Lk E implies
E
2 k
2 k
L0 E
L0 E
ln k
Example 2: If
L0
1
ln(2 )
and E 10 6 , calculate the number of iterations required to make
Lk E . This requires
1 6 10
ln k
ln(2)
ln 10 6 ln(2)
19.93
Thus k 20 (k is an integer) iterations are required. Note: useful numbers: 210 ~ 103 , 2 20 ~ 106 .
Note: the „error‟ in the mid-point as an approximation to the root introduces an extra ½, Ck P 1 / 2 Lk .
ERROR ANALYSIS
The Bisection algorithm was stopped when
b a E .
This gave us an error estimate for the
approximate root. If the approximate root is ≈1, then an error estimate of 10 6 , say, is reasonable. If however, the approximate root is ≈ 107 , then an error estimate of bad.
10
6
is very
ABSOLUTE ERROR
If y approximates x, then
y x
is the absolute error.
RELATIVE ERROR y x
If y approximates x, then
is the relative error. This definition requires x 0 . Often
x
relative error is multiplied by 100 to give percentage error. Class Work
If y approximates x, then compute the (i) Absolute Error, (ii) Relative Error and (iii) Percentage Error for the following: 1.
x 0.3000 x101, y 0.3100 x101 .
2.
x 0.3000 x10 3 , y 0.3100 x10 3
3.
x 0.3000x104 , y 0.3100x104
Solutions
1.
x
1
0.3000 x10 , y
Absolute Error
1
0.3100 x10
0.1 .
Relative error 0.3333x101 2.
x
3
0.3000 x10 , y
0.3100 x10
Absolute Error 0.1x10 4 Relative error 0.3333x101
3
3.
x
4
0.3000x10 , y
4
0.3100x10
Absolute Error 0.1x103 . Relative error
0.3333x10
1
Observation: Widely varying absolute errors for the same relative error. Generally, relative error is more meaningful while absolute error can be misleading. STOPPING CRITERIA
In the Bisection algorithm, the computation is „stopped‟ when
b a E .
Then the current
midpoint c is an approximation to the root p with Absolute Error
c P
1 / 2 a b
i.e. half the current interval length. Thus this stopping criterion gives an approximation with absolute error 1/ 2 E . ADVANTAGES AND DISADVANTAGES OF THE BISECTION METHOD Disadvantages
1. Slow convergence, a large number of iterations may be required to make the absolute error, c P , small. 2. Good intermediate approximations may be discarded when the root is close to the midpoint. Advantages
1. Always converge to a solution. This makes it an attractive technique for finding a „starting value‟ (initial guess) for faster methods. 2. Does not require derivative of a given function. 1.2 FIXED-POINT ITERATIVE METHODS 1.
First-Order Fixed-Point Iteration.
In this section, we explore a class of numerical methods for solving the equation f ( x) = 0, which given a first term x0, generates a sequence { xn} such that: (a)
xn+1 is obtained from xn (first-order recurrence relation) and
(b)
xn as
n
, where α is a root of the equation f ( x) = 0.
The requirements are to (i) generate a sequence based on the equation f ( x) = 0 and (ii) to ensure that xn
( xn approaches α) as
n
(n tends to ∞) i.e. convergence is reached.
Notation of a Fixed-Point Iteration
It is first necessary to develop the notation of a fixed-point. Example 1.2: Let a function, f , be defined by f ( x)
2x
2
x . The roots of f ( x) = 0 are x = 0
and x = 0.5. Consider the sequence generated by the recurrence relation: x
n xn
0.1,
x
1
n
2 x (1 x ),
0 0.100
n
n
n
1 0.180
0,1, 2,...
2 0.295
(1.2.1) 3 0.416
4 0.486
5 0.5000
6 0.5000
From the Table, the sequence approaches the root x = 0.5 where “it gets stuck”. On dropping subscripts from equation 1.2.1, the equation x 2 x (1 x ) is obtained.
y
y
y= x
y= 2 x(1-x)
0.25 0.5 -0.1
x
0.75 1.0
x
-0.2
a:Roots
b: Fixed-Point
In the figure (b), it is seen that the graph of y = x and y 2 x(1 x) intersect at x = 0.5. Thus when RHS of equation 1.2.1 is evaluated at x = 0.5, the value of xn+1 = 0.5 is obtained. All subsequent terms in the sequence,
x
n
2
,x
n
3
, ... etc. will equal to fixed value 0.5.
Definition
The value
x
is
a fixed point of the function g if α =
The above example suggest that if a sequence generated by a recurrence relation of the form xn+1 = g ( xn) has a limiting value, α, then α is a fixed point of the function g. Formulation of Fixed-Point recurrence Relation
We are set here to find a value of x : f ( x ) 0 is satisfied. To utilize recurrence relation, the equation must be written in the form x
g ( x)
.
This can be accomplished by adding x to each side of f ( x) = 0. i.e. Given
f ( x ) 0
(1.2.2)
Adding x to both sides of equation (1.2.2), we have x f ( x ) x
(1.2.3)
Rearranging (1.2.3), we have x x f ( x )
(1.2.4)
x = g ( x)
(1.2.5)
There are therefore infinity ways of arranging (1.2.2) in the form of (1.2.5). Example 1.2.2: For the equation f ( x) 0 , where f ( x) x3 4 x 2 10
The following are the three possible arrangements in the form of x = g ( x): (i)
x
(ii)
x
(iii)
x
3
10 4 x 2
1 2
10
10 4 x
x
3
g1 ( x)
g 2 ( x)
g 2 ( x)
The systematic iterative approach required to determine a zero of f follows the following steps: 1.
Write f ( x) = 0 in the form of x = g ( x).
2.
Define the first-order recurrence relation, xn+1 = g ( xn)
3.
Given x0, use the recurrence relation to generate the sequence { xn}.
Example 1.2.3. Find the roots of x 3 4 x 2 10 0 .
Applying steps 1-3 above, we have Step 1: f ( x) x 3 4 x 2 10 Step 2: defining the 1 st-orde recurrence relation, we have xn 1
xn 1
xn 1
3
1 2
10 4 xn 2 10 xn
10 4 xn
3
g1 ( xn )
(1.2.6)
g 2 ( xn )
(1.2.7)
g 3 ( xn )
(1.2.8)
Step 3: Starting with x0 = 1.5, we obtain the following sequences: n 0 2 3 4 5 6 7 8 9 10 11 12 13 14
g 1( xn) 1 1.8171 -1.4748 1.0914 1.7364 -1.2725 1.5215
g 2( xn) 1.2870 1.4025 1.3455 1.3752 1.3601 1.3678 1.3639
g 3( xn) 1.3484 1.3674 1.3656 1.3653 1.3652 1.3652 d e gr
t e n e .
v
e . n .
C
o y
gr v o
n d a
C er o l N A
1.3652
The sequences show a range of behaviour. The first sequence does not approach a limiting value. The second sequence attains a limiting value after 14 iterations. The third sequence attains a limiting value after just 5 iterations. Then the question to ask is: What property does g ( x) require near a fixed- point α for α to be the limiting value of a sequence generated by xn+1 = g ( xn)? The answer leads us to another topic. Fixed-Point Convergence
Convergence of the sequence { xn} generated by xn+1 = g ( xn) occurs if g ( x ) 1
Smaller g ( x ) implies faster convergence. Classification of Fixed-Points Class Attracting
Condition if g ( ) 1
Bahaviour Convergent
Indifferent
if g ( )
1
Uncertain
Repelling
if g ( )
1
Divergent
So if the slope of g is positive, then monotonic convergence or divergence res ults. If the slope of g is negative, then oscillatory convergence or divergence results.
Initial Guess
Fixed-Points methods require a good initial value x0. A simple solution is to find an interval [a, b] in which the original function f satisfies the sign property and then use the mid-point 1 x0 ( a b ) as the initial value. In other words, bisection method can be used in evaluating 2 an approximate starting value. Example 1.2.4:
Find a root of f ( x) x2 5 x 4 0 . The possible arrangements are: x 2 4
(i)
x g1 ( x)
(ii)
x g 2 ( x ) 5 x 4
(iii)
x
g3 ( x )
5
x
2
4x 4
Obtaining convergence condition, we have 2 x
(a)
g1 ( x )
(b)
g 2 ( x)
(c)
g3 ( x )
5
For (a): Solving
5 2 5 x 4 2x 4 g1 ( x) 1 ,
2 x
1
5
we have
4 x 2
1
25 x 2
25 4
0
5 5 ( x )( x ) 0 2 2
5
2
x
5 2
(a guess of x0
1.5 will be appropriate)
So the derivative is less than 1 in magnitude for all x such that
5
2
x
5 2
. Moreover, if
g ( x) is near zero in the entire region, the iteration converges quickly; if the derivative is
near 1 in magnitude, the iteration converges quite slowly.
For (b): Solving
g 2 ( x) 1 ,
we have
25
4(5 x 4)
1
25 4(5 x 4)
20 x 9
0.45
x
(a guess of For (c): solving
x
0
1.0 will
g3 ( x) 1 ,
2 x 4
(2 x 4)2 2 x 4
i.e.
2
(a guess of x0
1.3
5 2
we have
1
1
1 or
x 3
be appropriate)
2x 4
or x
1
3 2
5
x
2.0 will be appropriate)
2
ACCELERATED CONVERGENCE – AITKEN’S METHOD
From the Fixed-Point processes seen so far it is clear that the choice of g is instrumental in defining the convergence properties of the sequence { xn} generated by
xn 1
g ( xn ) .
We now
investigate how to accelerate the convergence of { xn}. The first approach uses the sequence { xn}itself and simply modifies the term xn to improve convergence. The Aitken‟s acceleration scheme is given as
* n 1
x
xn 1
( xn ) 2
2
xn 1
(1.2.9)
Where x n
2
xn 1
x
n 1
x
n
xn 1 2 xn
xn 1
Equation (1.2.9) is sometimes referred to as Aitken‟s Delta Squared method. Example 1.2.5: In Example 1.2.3, fixed-point iteration was used to find the real root of the
equation x3 4 x 2 10 0 , specifically, for the sequence
1 x
n 1
2
10
x
n
3
, the process took
14 iterations to converge to 4 d.p. i.e. x14 1.3652 . The Table below shows the associated Aitken‟ sequence { xn*} obtained by applying Equation (1.2.9). Convergence to 4 d.p. is now achieved after just 6 iterations. n 0 1 2 3 4 5 6
xn* 1.3619 1.3643 1.3650 1.3652 1.3652
xn 1.5 1.2870 1.4025 1.3455 1.3752 1.3601 1.3678
For n = 1, * * x x11 x2 x2 n 1
( x2 x1 )
*
* x 2
1.4025
2
x
2
2 x1 x0
(1.4025 1.2870)2 1.4025 2(1.2870) 1.50
1.3619
For n =2, x 1.3455 * 3
(1.3455 1.4025)2 1.3455 2(1.4025) 1.2870
1.3643
For n =3, x 1.3752 * 4
For n =4,
(1.3752 1.3455)2 1.3752 2(1.3455) 1.4025
1.3650
(1.3601 1.3752)
2
x 1.3601 * 5
1.3652
1.3601 2(1.3752) 1.3455
For n =5, x 1.3678 * 6
1.4
(1.3678 1.3601)2
1.3652
1.3678 2(1.3601) 1.3752
ACCELERATED CONVERGENCE – The ‘g ’- Factor
A second approach to accelerating convergence is based upon the choice of the function g . Earlier several ideas were established regarding the convergence of a sequence { xn} generated by
xn 1
g ( xn ) ,
if g ( )
to a fixed- point α of g , that
(a)
xn
(b)
Convergence is accelerated for smaller values of g ( ) .
1
and
Thus this suggests the following questions: (i)
Can g always be chosen such that
(ii)
Can the value of
g ( ) be
g ( )
< 1?
optimised (minimised)?
To answer the first question, let us add the quantity
x to
each side of the equation x
g (x) ,
where 1 , to obtain x x g ( x) x
Rearranging, we have x(1 )
x
g ( x ) x
g ( x) x (1 )
G ( x)
Equation (1.30) has the same root(s) as x
Example 1.2.6: The quadratic equation x
form x g ( x) , i.e. g (1)
g (2)
1 2 3
x
42 3
( x2 2) / 3 , then
1
2
(1.30)
2
g ( x)
.
3x 2
0 has 2 real roots, x =1 and x = 2. In the
And x = 1 and x = 2 are fixed-points of the expression
x
2
2
3
.
If, for example, 7 x is added to both sides and then arranged to the form equation (1.2.6), we obtain x x
2
2
3
7 x
1 7
x
2
21x 2
24
G x
G(1) = 1 and G(2) = 2. Thus, x = 1 and x = 2 are fixed point of G, and hence the roots of x = G(2). To find a fixed point α of G, the relation xn + 1 = G( xn) could be used to generate a sequence { xn} with limiting value α. To be efficient G ( ) should be small, ideally close to
zero. G x
g x x 1
G x
G
g x 1
,
1
(1.31)
Equals zero if g 0, i.e.
g
Typically, if a root is sought in the vicinity of x0, then λ is taken as
g
x . The better the 0
initial guess, the better the value of λ. Example1.2.7: consider again the equation x 3 4 x 2 10 applied to all three iterative processes of example 1.6.
i.e. xn
1
gi xn , I = 1, 2, 3 where
0.
Then acceleration scheme is
1 g i 15 .
Compare the convergence of the sequences in the table below (with 11, 5 and 3 iterations respectively) with the original table under example 1.6 (divergence, 4 and 6 iterations).
λ
4.000
0.6556
n
0.1226
xn
0
1.5000
1.5000
1.5000
1
1.4000
1.3713
1.3653
2
1.3785
1.3657
1.3652
3
1.3705
1.3653
1.3652
4
1.3674
1.3652
5
1.3661
1.3652
6
1.3656
7
1.3654
8
1.3653
9
1.3653
10
1.3652
11
1.3652
1.2.5 Optimized Accelerated Convergence – Newton Newton Raphson Iteration.
Let the sequence { xn} generated by the modified relation xn+1 = G( x) x) converge to a limiting value equal to the fixed point α of the equation x = x = g ( x). x). It is reasonable to assume that each term x term xn+1 will be an improved estimate to α when compared with the estimate xn.
To further accelerate the sequence, we might update the value of λ at each step using λ = fixed value λ = g ( xn ) in place of the fixed
xn should be closer to α than x0 for a g ( x0 ) ( x
convergent sequence). The effect of this modification is now investigated and it results in the new iterative scheme – scheme – Newton-Raphson Newton-Raphson iteration.
A recurrence relation needs to be developed:
xn
1
G xn
g xn xn 1
g xn g xn xn 1 g xn
(1.28)
The function g function g is is related to function f function f appearing appearing in the original equation f equation f ( x) x) = 0. So by simple manipulation, adding x adding x to to each side yields.
x f x x x x f x g x Hence, g(x) = x – x – f(x) f(x)
g1 x
f i x
1
Then equation (1.28) becomes
xn
xn
f xn
1
1
xn
1
f xn x n i
f f i xn xn
1
1
xn
f
i
f
i
xn
n
n
x x n
n
x
f xn
f xn i
n
n
f
i
f
f i xn
xn
x x x
xn
f xn f i xn
(1.29)
Equation (1.2.9) is called Newton-Raphson iteration scheme. Example 1.2.8: For the equation f equation f ( x) x) = 0 with
f x x 3 4x 2 10 ; f x 3 x 2 8 x .
The recurrence relation (1.29) becomes:
xn
1
xn
f xn f
i
x
xn
n
xn 3 4 xn 2 10 3 xn
2 xn 3 4 xn 2 10 3 xn 2 8 xn
Stating with x with x0 = 1.5 n = 0; x0 = 1.500
x0
1
x1
2 x0 3 4 x12 3 x0
2
8 x1
10
13733 .
2
4 xn
n =1; xi = 1.3733
x1
x2
1
2 x13 4 x12 10 3 x12 8 x1
13653 .
n = 2; x2 = 1.3653
x2
1
x3
2 x2
3
3 x2
4 x2 2
2
10
8 x2
13652 .
n = 3; x3 = 1.3652
x
31
2 x3
x4
3
4 x3
3 x3
2
2
10
8 x3
13652 .
Convergence of Newton-Raphson Method.
The general relation is xn 1
g xn
So it could be deduced from equation (1.2.9) that g x
x
f x
f i x
Recall that the derivative of g is the key descriptor of convergence points. Here g x i
1
f i x f i x
f
1
i
f x f ii x
x
f i x f i x 2
f i x
2
f x f ii x
f
i
2
x
f x f ii x
f
i
2
x
At a root x root x = α, the value f(α) is zero by definition. Hence, for the Newton -Raphson scheme, g x 0 - „optimal‟ convergence. So the condition for convergence:
g x
1 holds.
f x f x
f x
2
< 1
2.0
INTERPOLATION AND POLYNOMIAL APPROXIMATION
Engineers and scientists commonly assume that relationships between variables in a physical problem can be approximately reproduced from the data given by the problem. The ultimate goal might be to determine the values at intermediate points to approximate the integral or derivative of the underlying function, or to simply give a smooth or continuous representation of the variables in the problem.
Interpolation refers to determining a function that exactly represents a collection of data. The most elementary type of interpolation consists of fitting polynomials, so they are a natural choice for approximating derivatives and integrals. We will see here how polynomials are easily constructed, the result of which implies that there are polynomials that are arbitrarily close to any continuous function.
2.1
Langrange Polynomials
Here, we find approximating polynomials that can be determined simply by specifying certain points in the plane through which they must pass. Determining a polynomial of degree that passes through the distinct points ( x0, y0) and ( xi, yi) in the same as approximating a function f for which f ( x0) = y0 and f ( x1) = y1 by means of a first – degree polynomial interpolating, or agreeing with, the values of f at the given points. We first define the functions L0 x
x
x0
x1 x1
,
L1 x
x
x1
x0 x0
and then define
P x L0 x f x0 L1 x f x1 since
L0 x0 1,
L0 x1 0,
L1 x0 0, and L1 x1
1,
we have
P x0 1. f x0 0. f x1 f x0 y0 and P x1
0.
f x0 1. f x1 f x1 y1
So, P is the unique linear function passing through ( x0, y0) and ( x1, y1),
y y = f ( x ) y 1 = f ( x 1)
y 0 = f ( x 0)
y = P( x )
x
x 1
x 0
To generalize the concept of linear interpolation to higher degree polynomials, consider the construction of a polynomial of degree at most n that passes through the n + 1 points.
x , f x , x , f x ,....., x , f x . See the figure below. 0
0
1
n
1
n
y f
P
x 0
x n - 1
x 2
x 1
x n
x
In this case, we need to construct, for each k = 0, 1, … , n, a function Ln, k ( x) with the point that Ln, k ( xi) =0 when i k and Ln, k ( xk ) = 1. To satisfy Ln, k ( xi) =0 for each i = k requires that the numerator of L n, k (x) contains the term.
x x x x ...... x x x x ...... x x 0
k 1
1
k 1
n
To satisfy Ln, k (xk ) =1, the dominator of L n, k (x) must be equal to this term evaluated at x = xk . Thus Ln ,
k
( x) =
x - x0 ... x x 1 x x 1 ... x x x x0 ... x x 1 x x 1 ... x x k
k
k
k
k
k
k
n
k
A sketch of the graph of a typical L n, k is shown below
n
Ln, k ( x )
x k -1
x 1
x 0
x k
x k +1
x
x n -1 x n
The interpolating polynomial is easily described now that the form of Ln, known. This polynomial is called the nth Lagrange interpolating polynomial. Pn x
f x0 Ln , o x ...... f xn Ln ,n x
k
( x) is
n
f xn Ln k x ,
k 0
Where Ln, k x
x x x x ..... x x x x x ...... x
k
0
k
0
1
k
1
xk
1
xk
x x
1
k
xk
...... x x x ...... x x
1
k 1
n
k
n
For each k = 0, 1, ………., n. If x0, x1,..., xn are (n + 1) distinct numbers and f is a function whose values are given at these numbers, then P n( x) is the unique polynomial of degree at most n that agrees with f(x) at x0, x1, …, xn. The notation for describing Lagrange interpolating polynomial P n( x) is rather complicated. To reduce this somewhat, we will write Ln, k ( x) simply as Lk ( x). Example 2.1: Using the numbers, or node, x0 = 2, x1 = 2.5, and x 2 = 4 to find the second interpolating polynomial for f ( x) = 1/ x.
Solution: this requires that we first determine the coefficient polynomials L0, L1 and L2. L0 x
x 2.5 x 4 2 2.5 2 4
L1 x
x 2 x 4 2.5 2 2.5 4
x 2 x 2.5 4 2 4 2.5
x 6.5 x 10
4 x 24 x 32
3
And L2 x
x 4.5 x 5
Then f(x) = 1/x So, f x0 f 2 1 2 0.5
3
f x1 f x2
f 2.5
f 4
1 2.5 0. 4 and
14
0.25
We have 2
f xk Lk x
P x x
k
0
f x0 L0 x f x1 L1 x f x2 L2 x
05 . x 6.5 x 10 0.4
4 x 24 x 32 3
0.25
x 4.5 x 5 3
0.05 x 0.425 x 115 . An approximation to f(3) = 1/3 is
f 3 P 3
0325 .
Example 2.2: The table below lists values of a function at variable points. compare the approximation to f(1.5) obtained by Lagrange polynomials of define 2 and 3.
x
f ( x)
1.0
0.7651977
1.3
0.6200860
1.6
0.4554022
1.9
0.2818186
2.2
0.1103623
Since 1.5 is between 1.3 and 1.6, then the most approximate linear polynomial uses x0 = 1.3 and x1 = 1.6. The value of interpolating polynomial at 1.5 is
15 P . 1
15 15 . 1.6 . 1.3 0.6200860 0.4554022 1.3 1.6 1.6 1.3
=0.5102968 Two x
0
polynomials 1.3,
x
1
of 1.6,
degree x
2
2
could also 1.9 , which gives
be
used,
one
by
letting
P 2 1.5
1.5 1.61.5 1.9 0.6200860 1.5 1.31.5 1.9 0.4554022 1.3 1.6 1.3 1.9 1.6 1.3 1.6 1.9 1.5 1.31.5 1.6 0.2818186 1.9 1.3 1.9 1.6
= 0.5112857 In the third – degree core these are two reasonable choices for the polynomial. One uses: (i)
x0 = 1.3, P 3 15 .
x1 = 1.6,
15 . 1.3
x2 = 1.9,
1.6 15 .
1.9 15 .
2.2
1.6 13 .
1. 913 .
2.2
x3 = 2.2, which gives
0.6200860
15 . 13 . 15 . 19 . 15 . 2.2 0.4554022 16 . 13 . 16 . 19 . 16 . 2.2
15 . 13 . 15 . 16 . 15 . 2.2 0.2818186 19 . 13 . 19 . 16 . 19 . 2.2
15 . 13 . 15 . 16 . 15 . 19 . 0.403623 2.2 13 . 2.2 16 . 2.2 19 .
= 0.5118302 Although
2.2
P 3 15 . is
the most accurate approximate .
Langrange Polynomial Error Formula.
It has the same form as the error formula for the Taylor Polynomials. f x Pn x
f
n 1
x
n 1 !
x
x0 x x1 ... x xn
For some x between x0, x1, …, xn and x. A practical difficulty with Langrage interpolating is that since the error term in difficult to apply, the degree of polynomial needed for the desired accuracy is generally not known until computations are determined. The usual practical is to compute the results given from various polynomials until appropriate agreement is obtained, as was done in the previous example. Recursively generated Langrage Polynomials. Let f be defined at x 0, x1, …xk and x j, xi, be two numbers in this set. If
P x
x
x j P0,1,..., j
1 , j 1 ,...k
x
x
i
x xj
xi P0 ,1 ,...,i
1 ,i 1 ,...k
x
Then P ( x) is the k th Lagrange polynomial that interpolates f at the k th points x0, x1, …, xk . To test the correctness of the recursive formula, Let Q
P 0,1,...,i
Q P 0,1,..., j
P x
x
1,i 1,...k
1, j 1,...k
since Q( x) and Q( x) are polynomials of degree at most k – 1
x j Q x
and
x
i
x xj
xi Q x
Must be of degree at most k if 0 ≤ r ≤ k with r ≠ j, then Q( xr ) = Q( xr ) = f ( xr ), so P xr
x
r
x j Q xr
x
r
xi Q xr
x
i
xi
xj
xj
f
x r
f xr
This result implies that the approximations from the interpolating polynomials can be generated recursively in the manner shown in the table below: x0
P 0 = Q0, 0
x1
P 1 = Q1, 0
P 0,1 = Q1,1
x2
P 2 = Q2, 0
P 1, 2 = Q2, 1 P 0,1, 2 = Q2, 2
x3
P 3 = Q3, 0
P 2, 3 = Q3, 1 P 1,2, 3 = Q3, 2
P 0,1, 2,3 = Q3, 3
x4
P 4 = Q4, 0
P 3, 4 = Q4, 1 P 2, 3,4 = Q4, 2
P 0,1, 2,3,4 = Q4, 3
P 0,1, 2,3,4 = Q4, 4
This procedure is referred to as Neville‟s method.
Example 2.3: In Example 2.2, values of various interpolating polynomials at x = 1.5 were obtaining using the data shown in the first two columns of the table below. Suppose that we want to use Neville‟s method to calculate the approximate to f (1.5). If x0 10 . , x1 13 . , x2 16 . , x 3 19 . and x 4 2.2, then f (1.0) =
Q0,0, f (1.3) = Q1,0, f (1.6) = Q2,0, f( 1.9) = Q3,0, and f (2.2) = Q4,0, so these are the 5 polynomials of degree zero (constant) that approximate f(1.5). Calculating:
i
xi
f i = Qi, 0
0
1.0
0.7651977
1
1.3
0.6200860
Q1,1
2
1.6
0.4554022
Q2,1
Q2,2
3
1.9
0.2818186
Q3,1
Q3,2
Q3,3
4
2.2
0.1103623
Q4,1
Q4,2
Q4,3
Q4,4
Calculating the approximation Q1,1(1.5) gives Q1,1 15 .
15 . 10 . Q1,0 15 . 13 . Q0 ,0 13 . 10 .
05 . 0.6200860
0.20.7651977
0.3
0.523349
Similarly, Q2,1 15 .
1.5 1.3 0.4 554022 1.5 1.6 0.6200860 16 . 13 .
= 0.5102968 Q3,1 15 . 05132634 . . and Q4,1 15
05104270 .
The higher degree approximates are generated in a similar manner and are shown in the table below: 1.0
0.7651977
1.3
0.6200860
0.5233449
1.6
0.4554022
0.5102968
0.5124715
1.9
0.2818186
0.5132634
0.5112857
0.5118127
2.2
0.1103623
0.5104270
0.5137361
0.5118302
0.5118200
The best linear approximate is expected to be Q2,1 since 1.5 is between x1 = 1.3 and x2 =1.6
3.0
THE FINITE DIFFERENCE CALCULUS
In conventional calculus the operation of differentiation of a function is a well – defined formal procedure which the operations highly dependent on the form of the function involved. Thus we need a technique for differentiating functions by employing only arithmetic operations. The finite difference calculus satisfies this need. 3.1
Divided differences
Iterated interpolating was used in the previous section to generate successful higher degree polynomial approximations at a specific point. Divided – difference methods introduced here are used to successfully generate the polynomials themselves. We first need to introduce the
notation. The zeroth divided difference of the function f w.r.t. xi f xi is simply the value of f at xi
f xi f xi The remaining divided differences are defined inductively. The first divided difference of f w. r. t. x1 and xi+1 is denoted by f xi , xi f xi , xi 3.2
1
f xi
xi
1
1
1
and is defined as
f xi xi
Forward and Backward Differences.
Consider a function f ( x) which is analytic (can be expanded in a Taylor series) in the neighbourhood of a point x. We find f ( x + h) by expanding f ( x) in a Taylor series about x: f x h f x hf x
h2 2!
f x
h3 3!
f x ...
(3.1)
Solving equation (3.1) for f ( x) yields f x
f x h f x h
h
f x
2
h2 6
f x ....
(3.2)
Using the Error notation, f x
f x h f x h
h
(3.3)
In words, equation (3.3) states that we have found an expression for the 1 st derivative of f w. r. t. x which is accurate to within an error of order of h. we shall employ the subscript notation. f x h
f x f j
f j i
(3.4) (3.5)
Using the notation, (3.3) becomes f x
f j 1 f j h
h
(3.6)
We define the first forward difference of f at j as f j
f j 1 f j
(3.7)
The expression for f ( x) may now be written as f x The term
f j / h is
f j
h
(3.8)
h
called a first forward difference approximate of error order h to f ( x) .
Gradually, the expression f j 1 f j
h approximates the slope of the function f at the pt x by
the slope of the straight line passing through f ( x +h) and f ( x).
We now use the Taylor series expansion of f ( x) about x to determine f ( x – h): f x h f x hf
x
h2 2!
f
x
h3 3!
f x ...
(3.9)
Solving for f ( x) , we have f x
Or
f x f x h h
f x
h 2
f x
f x f x h h
h2 6
f x ...
h
(3.10)
(3.11)
Using the subscript notation, f x
f j
f j 1
h
h
(3.12)
The first backward difference of f at j is defined as f j
fj
f j 1
(3.13)
So that the expression (3.12) for f ( x) may be written as f x
f j
h
h
(3.14)
The term f j h is called a first background difference approximation of error order h to
f ( x) . The geometric interpretation of the approximate is that of the slope of the straight line
connectivity f ( x) and f ( x – h). 3.3
Approximations to Higher Order Derivatives
Considering Taylor series expansion (3.1) for f ( x + h), f x h f x hf x
h2 2!
f x
h3 6
f x ...
(3.15)
Performing a similar expansion about x, f ( x + 2h) is found as 2
f x 2h f x 2hf x 2h f
x
4 h3 3
f
x ...
(3.16)
Multiplying (3.15) by 2, and subtracting (3.15) from (3.16), we have f x
f x 2h 2 f x h f x h2
hf
x ...
(3.17)
Or, employing the subscript notation, f x
f j 2 2 f j 1 f j h2
h
(3.18)
We have now found an expression for the second derivative of f w.r.t. x which is accurate to within an error order of h. The 2nd forward difference of f at j is defined as
f j
f j 2 2 f j 1 f j
(3.19)
We may rewrite (3.18) for f ( x ) as f x
2
f j
h2
h
(3.20)
By using the backward expansion (3.9) to obtain f ( x - h) and a similar expansion about x to obtain f ( x – 2h), we can find backward difference expression for f ( x) which is accurate to h :
f x
f j
2 f j 1
f j 2
h2
h
(3.21)
The 2nd backward difference of f at j as defined as
2
f j 2
fj
2
f j 1 f j 2
Equation (3.21) may then be written as
(3.22)
note very well f x
2
f j
h2
you may be ask to obtain the expression h
(3.23)
We may now define the procedure for finding any higher forward and backward differences of order, say n. Any forward or backward difference may be obtained starting from the 1 st forward and backward difference (3.1) and (3.13) by using the following recurrence formulas.
n
f j
n
f j
n 1
f j
(3.24)
f j
(3.25)
n 1
e.g. the second backward difference of f at j may be found as
2
f j
f j
f
fj
j
f j 1 f j 1 f j 2
f j 1
fj
f
j
2 f j1
f j 1
f j 2
Forward and backward difference of expression for derivatives of any order are given by
note very well
d n f dx n
and
xj
n
f j
hn
h
(3.26)
find expression for delta, d n f dx n
xj
n
f j
hn
h
(3.27)
Forward and backward difference expressions of
h
are tabulated below for derivatives of
th
up to 4 order. It may be a convenient memory aid to note that the coefficients of the forward difference experience expressions for the nth derivative starting form j and proceeding forward are given by the coefficient of (- 1) n(a - b)n in order, while those for the backward difference expressions starting from j and proceeding backward are given by the coefficients of (a - b)n in order.
hf ( x)
f j
f j +1
-1
1
f j + 2
f j + 3
f j +4
h
h f ( x)
1
-2
1
h3 f ( x)
-1
3
-3
1
h4 f iv ( x)
1
-4
6
-4
2
(a) Forward difference representations
1
f j - 4
f j - 3
f j - 2
f j - 1
f j
-1
1
1
-2
1
-1
3
-3
1
-4
6
-4
1
hf ( x)
h2 f ( x) h3 f ( x) h4 f iv ( x)
1
h
(b) Backward difference representations 3.4
Higher Order Forward and Backward Difference Expressions
The difference expressions for derivatives which we have thus far obtained are of h . More accurate expressions may be found by simply taking more for example, the series in equation (3.1) for f ( x + h).
f x h
h2 h3 f x hf x f x f x ... 2! 3!
(3.28)
As before, solving for f ( x) yields f x
f x h f x h
h 2
f x
h2 6
f x ...
(3.29)
From (3.17) we have a forward difference expression for f ( x ) complete which its error term. Substituting this expression into (3.29), we obtain f x
f x h f x h
h f x 2 h 2 f x h f x
2
h2
h2 f x ... hf x ...... 6 (3.30)
Collecting terms, f x
f
x 2h 4 f x h 3 f x 2h
h2 3
f x ...
(3.31)
Or in subscript notation, f x
f j 2 4 f j 1 3 f j
2h
h
2
(3.32)
We have thus found a forward difference representation for the first derivative which is accurate to h
2
.
Note that the expression is exact for a parabola since the error term
involves only 3 rd and higher derivatives. A similar backward difference expression of h 2 could be obtained by using the backward Taylor series expansion of f ( x - h) and replacing
f ( x ) by the backward difference expression of
backward difference expressions of
for
h
2
h
form equation (3.21). Forward and
derivatives of up to 4 th order are tabulated
below. f j
f j + 1
f j + 2
2hf ( x j )
-3
4
-1
h 2 f ( x j )
2
-5
4
-1
-5
18
-24
14
-3
3
-14
26
-24
11
-2
f j - 2
f j - 1
f j
1
-4
3
-1
4
-5
2
3
-14
24
-18
5
11
-24
26
-14
3
3
2h f ( x j ) 4
h f iv ( x j )
f j + 3
f j + 4
f j + 5
h
2
h
2
(a) Forward difference representations
f j - 5
f j - 4
f j - 3
2hf ( x j ) h 2 f ( x j ) 2h3 f ( x j ) h 4 f iv ( x j )
-2
(b)Backward difference representations
Higher order forward and backward difference representations, although rarely used in practice, can be obtained by replacing successively more terms in the Taylor series expansions by difference representations of h .
3.5
Central Differences
Consider again the analytic function f ( x); the forward and backward Taylor series expansions about x are respectively. f x h f x hf x
f x h f x hf
x
h2 2!
h2 2!
f x
f
x
h3 3!
h3 3!
f x ...
(3.33)
f x ...
(3.34)
Equation (3.34) – equation (3.33), yields h3
f x h f x h 2hf x
3!
f x ...
(3.35)
Or solving for f x , f
x
f x h f x h
2h
h2 6
f x ...
(3.36)
Or f x
f x h f x h
h ...
2h
2
(3.37)
Employing subscript notation, f x
f j 1 f j 1 2h
h ... 2
(3.38)
Thus, difference representation is called a central difference representation and is accurate to
h 2 . Note that the expression is exact for polynomials of degree 2 (parabolas) and lower. 2 An expression of h for f x is readily obtainable from (3.33) and (3.34) by adding
these equation and solving for f x to yield. f x
f j 1 2 f j h
f j 1
2
h 2
The central difference expressions of
(3.39)
for
h
2
derivatives up to 4 th order are tabulated
below. A convenient memory aid for these central difference expressions of h2 in terms of ordinary forward and backward differences is given by d n f dx
n
d n f dx
n
n
f j n 2 n f j n 2 2h
n
f j n 1
2
2
2h
n
n
f j n 1
h , n
2
2
even
h , n 2
odd
(3.40)
(3.41)
f j – 2
f j – 1
f j
f j + 1
-1
0
1
2hf ( x j )
f j + 2
h
1
-2
1
-1
2
0
-2
1
1
-4
6
-4
1
2
h f ( x j )
2hf ( x j ) 4
h f iv ( x j )
Representation of
3.6
h
2
2
Difference and Polynomials
Difference expressions for derivatives and polynomials have some distinct relationships which can be very useful. Thus, if we consider a polynomial of order n, the nth difference representation taken anywhere along this polynomial will be constant and exactly equal to the nth derivative regardless of the mesh spacing h (since all of the error terms will be zero). This knowledge may be used to get some idea of how well a given polynomial will fit data obtained at a series of equally-spaced points on the independent variable. For example, if the 3rd differences taken at various values of the independent variable are approximately equal and the 4th differences are close to zero, then a cubic polynomial should fit the date relatively well. Example 3.1: find the fifth backward difference representation which is of h .
Solution From the recurrence scheme for difference, the fifth backward difference can be expressed as
5
f j
4 f j 1
f j
5 f j 1 10 f j 2 10 f j 3 5 f j 4 f j 5
f j
4
f j
6 f j 2
4 f j 3
f j 4
f
f 1
4
f j 2
6 f j 3
4
f j 4
f j 5
and d 5 f dx 5
5 f j h5
h
Example 3.2: Given the function tabulated at the points j, j +1, and j + 2 shown in the figure
below, find a three points difference representation for f j .
f ( x )
h
2h
j
j +1
x
j +2
Solution: We pass the parabola f x Ax 2 Bx C through the points x j = 0,
x j + 1 = h, x j + 2 =3h and solve for f (0) : f j f j
1
f j
2
(a)
C
Ah 2 Bh C
9 Ah
2
3
b c
Bh C
f x 2 Ax B f 0 B
Now f 0 B , so solving for B yields 2
f j 1
Ah
f j
9 Ah
2
2
Bh f j
(d )
3Bh
f j
(e)
Eqn (e) – 9 x Eqn (d ) yields
f j
2
9 f j 1 6Bh 8 f j
B
f 0
9 f j 1
fj
2
8 f j
6h
8 f j 9 f j 1
f j 2
6h
Example 3.3: 2 Find a central difference representation of h for d 5 f / dx 5 . From example 3.1,
5
f j
fj
5
fj
1
10 f j
2
10 f j
3
5 fj
4
f j
5
5
f j can be found in a similar manner as
5
f j
f j
5
f j
5
4
f j
4 f j4
6 f j 3
5 f j4
10 f j 3
4 f j 2
fj
10 f j 2
1
f
5 f j1
j 4
4
fj
3
6 f j 2
4
fj
1
f j
f j
Applying equation (3.41) directly d 5 f dx
5
5
f j
2
2h
5
f j
2
2
h 2
Or d 5 f dx5
{ f j 2 5 f j 1 10 f j 10 f j 1 5 f j 2 f j 3 [ f j 3 5 f j 2 10 f j 1 10 f j 5 f j 1 f j 2 ]} / 2h5 (h)2
f j 3
4 f j2
5 f j 1
5 f j 1
2h
4 f j2
f j 3
5
h
2
Example 3.4: Given the following equally spaced data:
x
0
1
2
3
4
f ( x)
30
33
28
12
-22
Find f (0) , f (2) , f (4) and f (0) using difference representations which are of
h
2
.
Solution: At x = 0, a forward difference must be used since no points are available in the backward direction. f 0
f 0
f 2 4 f 1 3 f 0 2 1
28 4 33 3 30
2
7
1
2
to 1
2
At x = 2, we have a choice of several representations. We arbitrary select a central difference representation of f 2
h
2
f 3 f 1 2
1
2
12 33 2
10.5
to
1
2
At x = 4, a backward difference representation must be employed: f 4
3 f 4 4 f 3 f 2 2 1
3 22 4 12 28 2
43
1
2
to 1
2
At x = 0, the forward difference representation of f 0
2 f
0 5 f 1 4 f 2 f 3 2
1
3 30 5 33 4 28 12 1
5 to 1
2
h
2
h
h
2
2
is given as
4.0
DIFFERENCE EQUATION
Definition
Suppose U 1, U 2, U 3,…, U n are the terms of a sequence. We now introduce a difference operator such that the result of operating with on U r , is defined by U r Ur 1 U r
(4.1)
Where r = 1, 2, 3, …, (n - 1), hence
(4.2)
U 2 U3 U 2
(4.3)
U1 U 2 U 1
And so on. There expression are usually called the first finite difference of U 1, U 2, U 3, . . ., U n. In the same way the second finite differences
2
U r Ur
U
r 2
Ur
2
Ur
2U r
1
U
U
1
r 1
r 1
Ur
Ur
2
U r are
defined by
(4.4)
U r
Consequently,
2
2
U1 U3 2U 2 U1 U2
(4.5)
U 4 2U 3 U 2
Higher order finite differences
(4.6)
3
4
n
U r , U r ,..., U r may
be defined in a similar way.
Any equation which expresses a relation between finite differences 2
m
U r , U r ,..., U r is
called a difference equation, the order of the equation being
equal to the higher order finite difference term contained in it. e.g. (is of order 1)
U r 5U r 3 2
Ur
3U r 2U r
3
Ur
r U r
r2
2r U r 0
(is of order 2) (is of order 3)
The general mth order linear difference equation therefore has the form.
a0 mU r
a1
m 1
Ur
a2
m 2
Ur
......am 1 U r amUr
f r
(4.7)
Where f (n) and a0, a1, …, am are given functions of r . by analogy which the terminology of o.d.e (4.7) is said constant coefficient type if a0, a1, ………am are numbers in dependent of r .
Difference equations may always be written in a form involving successive values of U r by substituting the definition of
U r ,
2
U r ,...
e.g. using equation (4.1) the equation, Ur
5U r 3
becomes
U r 1 U r
5U r
U r 1 6r
3
3
(4.8)
Also the equation, 2
U r 3 Ur 2Ur r
2
becomes
U
r2
2U r 1
Ur
3U
r 1
Ur
2U r
U r 2 5U r 1 6U r r 2
r
2
(4.9)
When written in this form, difference equations are often called recurrence relations. 4.2
Formation of Difference Equations The following examples illustrate the formation of difference equations from expressions defining U r as a function of r .
Example 4.1:
If U r A4 r where A is an arbitrary constant, then U r 1 A4r 1
A.4.4r
4.A 4r 4U r
Hence we see that the elimination of the arbitrary constant A leads to a first order difference equation.
In general, if the expression for
U r contains
m arbitrary constants, it must satisfy mth order
difference equations.
Example 4.2: if r
U r A B(3)
(4.10)
where A and B are arbitrary constants, then r 1
U r 1
A B (3)
U r 2
A B(3)
(4.11)
And r 2
(4.12)
Hence eliminating B from (4.11) and (4.12) using (4.10), we have
U r A 3U r 2 A r 3
(4.13)
U r A r 9U r 8 A 3
(4.14)
r 1 U r 1 A 3
And U r 2
r 2
A3
Finally, eliminating A from (4.13) and (4.14) we get
U r 2 4U r 1 3U r 0
(4.15)
Consequently, (4.10) which contain 2 arbitrary constants must satisfy a 2 nd order difference equation. It should be clear that the general solution of an mth order linear difference equation must contain m arbitrary constants.
4.3
Solution of Difference Equations
1.
Homogeneous Linear 1 order Equation
st
Ur
1
arU r 0
(4.16)
Where ar is a given function of r . By putting r = 1, 2, 3, …, in (4.16), we have U2
a1U
U3
a2U 2 a2a1U1
(4.17) (4.18)
U4
a3U 3
a3a2a1U1
(4.19)
And so on, giving finally Ur
ar 1ar 2... a3a2 a1U1
(4.20)
Hence, if U 1 has a known value, say c, the solution of (4.16) may be written as r 1
U r c ap
(4.21)
p 1
2.
st
Inhomogeneous Linear 1 order equation
The equation U r 1 arUr br
(4.22)
Where ar , br are given functions of r , we define the general solution of (4.16) as Ur
Ur
wr
(4.23)
Where V r is the several solution of the homogeneous equation
Vr 1 ar vr 0
(4.24)
And where wr is any particular solution of the inhomogeneous equation Wr 1 ar wr
br
(4.25)
Example 4:3: Solve the homogeneous equation r
U r 1 4 U r
Given that
U 1
0
(4.26)
2
We write U2
4U 1,
U3
4
U2
4
U4
43 U 3
43 42 4U 1
2
(4.27) 2
4 U 1
Consequently, U r 4r 1 4r 2 ...43 42 4 U 1
r 1
2 4
p
p 1
Example 4:4: To solve the inhomogeneous equation U r 1 2Ur r
(4.28)
The homogeneous equation is given by
U r 1 2V r 0
(4.29)
and the inhomogeneous equating is Wr 1 2Wr
r
(4.30)
The general solution, therefore, is Ur
Vr
W r
Equation (4.29) is similar to (4.16), so the solution is given by
V2
2V 1
V3
2V2
2 2V 1
V4
2V3
2 2 2 V 1
Vr
2V r
U r 1
2
V 1 2
r 1
V 1
r 1
where V 1 is arbitrary. Hence the solution of (4.28) now depends on our ability to find any one solution of (4.30). To do this we adopt a trial solution of the form Wr
r
(4.31)
where α and are to be determined such that (4.31) satisfies (4.30) we find
r 1 2
v r
which gives comparing coefficients 0,
(coefficient of r 0)
(coefficient of r )
Hence 1,
1,
(4.32)
Wr
1 r
Consequently the general solution is U r V1 2 r 1 1 r
If in addition we are given that
U1 C
U1
(4.33)
C then
from (4.33)
C V 1 20 1 1
V 1
2
V1
U r
C 2
2 c 2r 1 1 r
5.0
NUMERICAL INTEGRATION
Here we are concerned with obtaining estimates of the value of the definite integral b
I
f x dx
(5.1)
a
Where f is a function of one real variable x; a and b define the lower and upper limits of the domain of integration on the x – axis. A formal evaluation of equation (5.1) requires finding a primitive F which, when differentiated, gives the original function f , i.e. dF dx
f
In this case, I
b
dF
a
dx
dx F b f a
(5.2)
A major problem with integration is finding a primitive F , and for this reason, the numerical integration is required in order to obtain numerical estimates to the value of I . A useful first step is to develop a visual appreciation of integration. It is often stated that the integral of a function f between limits x = a and x = b is equal to the area under the graph of y = f ( x) on the interval a x b . The value of the integral (5.1) equals the area under the graph y = f ( x). y
y = f ( x )
a
b
x
b
a
f x dx
Numerical integration can be thought of as the development of numerical methods for estimating the area of the shaded region.
5.1
THE TRAPEZOIDAL RULE
We wish to evaluate the integral I
b
a
f x dx
For an integreable function f ( x) on the interval a n equal subintervals each of width x where x
x b.
We divide the interval
a x b into
b a n
Thus, the trapezoidal rule is presented as n 1 x 2 I f 0 f n 2 f j x 2 j 1
(5.3)
The rule is thus termed, a second order method of numerical integration.
For most reasonably well – behaved functions, it is possible to obtain a much improved integration technique by estimating the error term involved in (5.3) and this result in 2
n x x I f f fj 2 n f b f a 2 12 j 1
0
1
Equation (5.4) is called the trapezoidal rule with end correction.
Example 5.1: use trapezoidal rule to evaluate I
sin xdx
0
using 3 panels and then include the end correction. Solution: f ( x) = sin x n = 3;
x
b a n
0
3
3
5.4
i
xi
sin( xi)
0
0
0.00000
1
3
0.866025
2
2 3
0.866025
3
0.00000
The trapezoidal rule yields x
I
2
3
2
6
f
0
f 3 2 f 1 f 2
. . 0866025 0.00000 0.0000 2 0866025
1813799 . . 4 0866025
The end correction is
x 12
2
f b f a
f x cos x
3
2
12
f f 0
3
cos
12
3
2
12
cos 0
2
2 0182770 .
Adding the end correction to the previously obtained trapezoidal rule value yields I
0182771 199659 1813799 . . .
This value is indeed closer to the analytical solution ( I = 2) than the previous results.
5.2
SIMPSON’S RULE
Simpson‟s rule is a numerical integration technique which is based on the use of parabolic arcs to approximate f ( x) instead of the straight lines employed as the interpolating polynomials with the trapezoidal rule. The method is given as n1 n 2 x 6 I f0 f n 4 f j 2 f j x 3 j 1 j 2
(5.5)
Equation (5.5) is Simpson‟s one – third rule for the entire interval. It is a fourth order method. The Simpson‟s 3/8 th rule is obtained by substituting the factor
x
3
with
3 x 8
in (5.5) with
which we obtain
3 x
I
8
f 0 fn 4
n 1
f
n2
j
2
j 1
j 2
j odd
j even
f j
(5.6)
With the ends correction taken into considering equation 5.5 becomes
n n x 1 I 14 f f n f j 16 f j x f a f b 15 j j 2 2
1
0
2
5.7
1
j even
j odd
Example 5.2: Using the same integrand in example 5.1, evaluate I using Simpson‟s 1/3 rd and 3/8th Rule with 3 panels and compared the results.
Solution: I
0
f x dx
sin xdx 0
x
3
Applying the Simpson‟s 1/3 rd rule. x
I
3
3
3
9
f
0
f 3 4 f 1 2 f 2
0 0 4 0.866025 2 0.866025
6 0866025 .
3
1.732050
18145292 .
The 3/8th Rule: I
3 x 8
8
f
0
f3
4 f 1 2 f 2
6 0866025 .
2.041344
5.3
ROMBERG INTEGRATION
This powerful and efficient integration technique is based on the use of trapezoidal rule combines with Richardson extrapolation. In order to describe the algorithm in detail, we adopt a new notation. The trapezoidal rule estimates of the integral will be denoted as l x 2 Tl k f a f b f a j x 2 j
(5.8)
.
1
Where x b a 2
k 1
and
l
2
k 1
1.
The number of panels, n, involves in T l, k is
2k – 1,
i.e. n = 2k -1. Thus T 1,1
T 1,2
T 1,3
b a 2
f a
f b
b a
b a
b a
f a f b 2 f a
4
b a
f a f b 2 f a
8
2
3 b a b a 2 f a 2 f a 2 2 4
e.t.c.
NOTE THAT :
T 1,2
T 1,3
b a b a f a 2 2
b a b a 3 b a f a f a 4 4 4
T 1,1 2
T 1,1 2
e.t.c. The extrapolating is carried out according to Tl ,k
1
4
l 1
4 1
l 1
Tl 1,k 1 T l 1,k
For example, for l = 2, T2 ,1
T2 ,2
Now for l = 3,
1
3
1 3
4T
T 1,1
1,3
T 1,2
1,2
4T
(5.9)
T3,1
1
15
16T
2 ,2
T 2 ,1
These results can conveniently be arranged in tabular form. T 1,1 I n c r e a s in
T 1,2
T 2 ,1
T 1,3
T 2 ,2
T 3,1
T 1,4
T 2 ,3
T 3, 2
T 4 ,1
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
T l 1,1
T 1,l
T 2 ,l 1
T 3,l 2
.
.
.
T l 1,2
g ly a c c u r a t e T r a p e z o id a l r u le v a lu e
T l ,1
Example 5.3: Consider the integral
5 x 4 I 4 x 3 2 x 1 dx 0 8 8
The analytical solution of the integrand is I = 72. Romberg extrapolation should yield this exact answer in only a few extrapolations. f x
5 x 8
4
4 x 3 2 x 1
So, f (0) = 1 f (8) = 2560 - 2048 + 16 + 1= 529 and b – a = 8 – 0 = 8 Trapezoidal approximate with 1 and 2 panels are T 1,1
8 2
1 529 2120
T1,2
2120
2
8 2
f 4
1060 4 160 250 8 1
712
Extrapolating these 2 values yields, T 2 ,1
1 3
4 712 2120
242
2 3
The trapezoidal rule with four panels gives T1,3
712 2
8 4
f 2
f 6 356 2 17 41 240
Extrapolating T 1,2 and T 1,3 yields T 2 ,2
1 3
4 240 712 82
2 3
By extrapolating T 2,1 and T 2,2 according to (5.9) we obtain T 3,1
1
2 2 16 82 242 72 . (which is the exact answer) 15 3 3
The Romberg table obtained so far, is 2120 2 242
712 240
82
3
2
72
3
The best available trapezoidal rule value of 240 using four panels is still very far from correct and the greatly accelerated convergence along the diagonal should be apparent. In general of course, we would not know that the exact answer had been obtained, so another line of the table would have to be computed. After this computation, the table would be 2120 2
712
242
240
82
128
1 2
72
2 3 2 3
3
72
72
72
6.0
NUMERICAL SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONS
All problem involving ordinary differential equations fall into two categories: the initial value problem and boundary value problems. Initial value problems are those for which conditions are specified at only one value of the independent variable. These conditions are termed initial conditions. A typical IVP might be of the form. A
d2y
B
2
dt
dy
dt
Cy g t ,
y 0
y0 ,
dy dt
0 V 0
On the other hand, Boundary value problems are those for which conditions are specified at two values of the independent variable. a typical BVP might be of the form.
d2y dx
2
dy
D
dx
Ey h x ,
y 0
y0 ,
y L
Y L
The problem is a boundary value problem if any conditions are specified at two different values of the independent variable. Thus, d4y
A y
dx 4 y 0
f x dy
y0 ,
dx
0 W 0,
d2y dx 2
0 V0 , y L
y L
is a boundary value problem.
6.1
The General IVP
Any IVP can be represented as a set of one of more coupled first – order ordinary differential equations, each which an initial condition. For example, the simple harmonic oscillator described by 2
A
d y dt
2
y 0
dy dt
B
Cy
dt
g t
y0
0
dy
V 0
(6.1) (6.2) (6.3)
Can be related by making the substitution Z
dy dt
The differential equation (6.1) can now be written as
(6.4)
A
dz dt
Bz Cy g t
(6.5)
With some rearrangement, the problem separated by equations (6.1) – (6.3) can now be written as dy dt
dz
z
dt
B A
(6.6)
z
t C y g A A
(6.7)
With initial conditions (i.cs.): y 0 y0
(6.8)
z 0 V 0
(6.9)
Any nth order differential equations can similarly be reduced to a system of n first – order differential equations.
The general form of any IVP can thus be stated as dy1 dt
dy2 dt
f1 y1 , y2 ,..., yn , t
f 2 y1 , y2 ,..., yn , t
.
.
.
.
.
.
dyn
dt
(6.10)
f n y1 , y2 ,..., yn , t
Subject to the initial conditions:
y 0 y1 0 2
y10 y20
. .
yn 0 yn0
(6.11).
Since any IVP can be expressed as a set of first – order ordinary differential equations, our concern now will be to consider numerical methods for the solution of first – ordinary differential equations. We will be concerned with two classes of methods. The first of these consists of the formulas of the Runge – Kutta type. In these formulas, the desired solution y j 1
is obtained in terms of y j between t j
t j
,
f y j t j and f y t evaluated for various estimated values of y ,
,
and t j 1 . Thesemethods are self starting because solution is carried directly from
to t j 1 without requiring values of y or f y t for t ,
t j .
y
Solution Known to here
Δt
y j ●
Solution desired here
●
t
t j + 1
t
A second class of methods consists of formulas of the multistep types. These formulas, in general, requires information for t t j . (see Fig below)
y
Δt
y j - 3
●
Δt
●
Δt
Δt
y j - 2 ●
y j - 1 ●
y j ●
t j - 3
t j - 2
t j - 1
t j
t j + 1
y j + 1 t
The solution for y j +1 might require the value of y j and values of f ( y, t ) at each of the points t j, t j -1, t j – 2 and t j – 3. These multi step formulas are obviously not self – starting.
6.2
THE EULER METHOD
Consider again the 1 st – order IVP dt
dy
y 0
f y, t
(6.12)
y0
The numerical solution could be obtained by replacing
dy
by a simple forward difference dt
representation. Thus (6.12) can be approximate by y j 1 y j t
f y j , t j
(6.13)
Solving for y j 1 yields
y j 1 y j
t x
f y j ,t j
(6.14)
Equation (6.14) is called the Euler formula.
Example: 6.1: stating at t = 0 and Δt = 0.1, solve the ordinary differential equation dy dt
y 2
(6.15)
The problem can be solved numerically by applying Euler recurrence formula (6.14). y j 1 y j
t y j
2
Stating at t = 0 ( j = 0) as Δt = 0.1, we can find y at t = 0.1. y1
1 0.9 2
1 0.1
The exact solution at these points is yexact 0.1
1
1 0.1
0.9090909
At t = 0.2
y2
0.9 0.1
0.9
2
0.819
The exact solution is yexact 0.2
1 1 0.2
0833333 .
(6.16)
After 10 steps of Δt = 0.1, we find y10 0.4627810
yexact
1 11
0.5
Although this numerical method is very simple, it is obviously not extremely accurate. However, because it is so simple, it is convenient to use as an introduction to numerical technique for ordinary difference equations.
6.3
RUNGE – KUTTA TYPE FORMULAS
These formulas are among the most widely used formula for the numerical solution of ordinary differential equations.
Their advantages include:
they are easy to program the step size can be changed as desired without any complication. they are self – starting they have good stability characteristics
Their disadvantages are:
they requires significantly more computer time than other methods of comparable accuracy. local error estimates are somewhat difficult to obtain.
For a n – step Runge – Kutta process, we may write n ki f x Ci h, y h aij k j , i = 1, 2, …, n j 1
(6.17)
and n
y y h bi k i l 1
The coefficients of the common procedures are listed in the table below
(6.18)
c
a
a
c
a
a
1
11
21
2
12
22
.
.
.
a1n
.
.
.
a
2n
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
a
.
.
.
bn
c
a
a
b1
b2
n1
n
n2
nn
(a)
0
0
0
0
1
1
1
0.5
0.5
0
0
0.5
0.5
Modified Euler
0.1
method
improved Euler method
Euler method
(d) (c)
(b)
0
0
0
0
0.5
0.5
0.5
0.5
1
-1
0.5
0
0.5
1
0
0
2
2
1
2
1
6
3
6
1
1
1
1
1
6
3
3
6
Runge-Kutta Method
Runge-Kutta Method
(e)
(f)
Employing equations (6.17) and (6.18), we can formulate the formulas presented in table ( a) to ( f ), e.g. (a)
Euler Method k1
f x 0.h, y h 0.k1
f x y
,
y y h 1. k 1
yi
(b)
y hf x y ,
1
yi
hf xi , yi
Improved Euler method
k1
k2
f x, y f x h, y k1h
y y h0.5k1 0.5k 2
(c)
y
h 2
k
1
k 2
Modified Euler method
k1
f x, y
k2 f x
1 2
h, y
1 2
k 1h
y y hb1k1 b2 k 2
(d)
y hb2 , k 2
y hk 2
Standard Runge – Kutta method
k1
f x, y
k2 f x
1 2
h, y
1 2
hk 1
