Newton Forward and Backward

BSc Computing Mathematics
BSc Computing Maths with Business Studies
Year 2
MATHEMATICAL TECHNIQUES:
NUMERICAL METHODS
CONTENTS
TOPIC
Interpolation
Difference Tables
Newton-Gregory Forward Interpolation Formula
Newton-Gregory Backward Interpolation Formula
Central Differences
Numerical Differentiation
Numerical Solution of Differential Equations
Euler's Method
Improved Euler Method (IEM)
Runge-Kutta Method
Page
4
6
8
13
16
21
26
26
33
39
NUMERICAL ANALYSIS
When handling problems using mathematical techniques it is usually necessary to establish a
model, and to write down equations expressing the constraints and physical laws that apply.
These equations must now be solved and a choice presents itself. One way is to proceed using
conventional methods of mathematics, obtaining a solution in the form of a formula, or set of
formulae. Another method is to express the equations in such a way that they may be solved
computationally, ie by using methods of numerical analysis. This will lead directly to
quantitative results, however if enough such results are obtained then qualitative results may
emerge. Using these methods, large and complex physical systems may be modelled, and
otherwise intractable problems can be solved with the help of modern computer systems.
Interpolation
We sometimes know the value of a function f(x) at a set of points, without knowing the analytic
expression for f(x) that would allow us to calculate its value at some arbitrary point. Often the
points at which we know the value of f(x) are equally spaced, but not always. The task is to
estimate f(x) for a given value of x by, in some sense, drawing a smooth curve through (and
perhaps beyond) the known values. If the desired value is in between the smallest and largest
x-values already know, then the process is called interpolation. If x is outside that range, it is
called extrapolation, which is considerably more hazardous. As a simple illustration, let us
consider linear interpolation.
Before the advent of computers, if it was required, for example, to find the square root of a
number x, a table of such numbers was consulted. If the number did not appear in the table,
then the two numbers above and below x were used, and interpolation provided the solution.
For example, suppose we wanted the square root of 2.155 and had a table such as shown
below:
n
2.140
2.150
2.160
2.170
n
1.46287
1.46629
1.46969
1.47309
The difference between the square roots of 2.15 and 2.16 is 0.0034. Since 2.155 is half-way
between 2.15 and 2.16 we could assume that its square root was half-way between the square
roots of 2.15 and 2.16. This gives 2.155 = 1.46799.
Had we wanted 2.153 we would have added 0.3 times 0.0034 to 2.15. This is the process of
LINEAR INTERPOLATION, in which we assume a linear variation between the two known
values to predict intermediate values.
Suppose we have a set of m numbers x1 , x2 , x3 ,... , xm and that corresponding to each xi there is
a yi , where y and x are related by some function f: y = f ( x ) . If we are given a value of x (not
equal to one of the xi ) we want to find a value of y obtained by linear interpolation. Naturally
we assume some functional relationship does indeed exist between x and y, and it is our
responsibility to ensure that this is a valid assumption. Suppose that the required value of the
variable x falls between the kth and k+1th known values of x [we assume that the x's are
ordered in an ascending sequence so that x1 < x2 < x3 < ... < xm ]
x k x xk + 1
BSc 2 Computing Mathematics / Computing Maths with Business Studies - Mathematical Techniques
and we want to calculate the corresponding value of y. To do this, we draw a line segment
joining the points ( x k , y k ) and ( x k +1 , y k +1 )
k+1
y
yk
xk
x k+1
We then find the point on that line with an x-coordinate of x .

The equation of the line through ( x k , y k ) and ( x k +1 , y k +1 ) is found by equating the slope of AB
to that of BC:
i. e.
y y k y k +1 yk
=
x xk
x k +1 x k
y = yk +
( x x k )( y k+1 yk )
xk +1 xk
Therefore, to find the required value of y we insert x for x:

y = yk +
( x x k )( y k+1 yk )
xk +1 xk
(*)
Consider a simple example:

Suppose we have tabulated the values of y = x 2 for all integer x. Find 6. 52 .
x k = 6, yk = 36, x k +1 = 7 , y k +1 = 49
Using (*)
y = 36 +
and x = 6. 5
(6.5 6)(49 36) = 36 + 05. 13 = 42.5

76
The correct value of 6.5 is 42.25 the error ET = -0.25.

[Error = true value - calculated value].
A General Approach to Interpolation

We could investigate methods of interpolation using higher order polynomials, but although
quadratic interpolation, for example, is likely to be more accurate than linear interpolation, it
is by no means certain to provide sufficient accuracy all the time. We need to develop a
general method which will also enable extra accuracy to be attained without having to resort to
a new set of calculations. In order to do this we need to introduce the idea of a difference
table.
Difference Table
As before, suppose that we have a set of values of an unknown function y1 , y2 , y3 ,..., yn
corresponding to a set of values of the (known) independent variable x1 , x2 , x3 ,..., xn where
x1 < x2 < x3 < ... < xn .
A reliable formula for interpolation, to obtain the value of y corresponding to any x (where
x1 < x2 < x3 < ... < xn ) will employ all the given information - ie it will be based on the values
of the n+1 data entries. Instead of using these directly it is more convenient to use
relationships dependent on these values. The relationship
yr +1 yr
is called the first difference of yr , denoted by yr :
yr = yr +1 yr
Similarly, the relationship yr +1 yr is the second difference:
2 yr = yr +1 yr
For example, 2 y0 = y1 y0 = ( y2 y1 ) ( y1 y0 ) = y2 2 y1 + y0
Similarly
2 y1 = y2 y1 = ( y3 y2 ) ( y2 y1) = y3 2 y2 + y1
In general terms:
2 yr = yr +2 2 yr +1 + yr
The third differences are
3 yr = yr +3 3 yr +2 + 3 yr +1 yr
r = 0, 1, 2 ,..., n 2
r = 0,1, 2,..., n 3
Clearly, successive differences can be calculated and put into a table such as that shown
below:
x
x0
y0
x1
y1
x2
x3
y2
y3
x4
y4
x5
y5
y
y0
y1
y2
y3
y4
2 y
2 y0
2 y1
2 y2
2 y3
3 y
3 y0
3 y1
3 y2
4 y
4 y0
4 y1
5 y
5 y 0
The polynomial interpolation formula, dependent on the n+1 entries, can be expressed in terms
of these differences.
x
-2
2 y
3 y
1
1
-1
2
3
14
23
34
0
2
0
2
0
2
0
2
11
In this example, all the differences after the second are zero. This is no coincidence - the nth
differences of an n-degree polynomial are always constant, higher differences being zero.
Errors in a Difference Table

The first step in any interpolation procedure is the formation of a difference table from the
given data. Since this often involves a large number of additions and subtractions it is an
inevitable source of error. It is of interest to note how errors propagate in a difference table.
Suppose y is subject to an error :
x
x0
y0
x1
y1
x2
x3
x4
y2
y3 +
y4
x5
y5
x6
y6
2 y
y0
2 y0
y1
y2 +
y3
2 y1 +
2 y2 2
2 y3 +
y4
2 y4
y5
3 y
3 y0 +
3 y1 3
3 y2 + 3
3 y3
The error in the third difference column, for example, appears in four terms, with error
3
coefficients given by the binomial coefficients of (1 ) .
The Newton-Gregory Interpolation Formula

The unique polynomial satisfying any n+1 values is of degree n:
p( x ) = A0 + A1 x + A2 x 2 + ... + An x n
The n+1 constants Ai being found from the n+1 linear equations obtained by substituting in the
data.
The Forward Interpolation Formula
Suppose the (n+1) values are (xi,yi) i=0,1,2,...,n. We can write the polynomial as:
p( x ) = a 0 + ( x x0 ) a1 + ( x x0 )( x x1 ) a2 + ...
If the equal intervals between successive values is h, then xr xs = ( r s) h .

We then have the n+1 equations for the ai:
yr = a0 + rha1 + r (r 1)h 2 a2 + ... + r (r 1)(r 2 )...(r n + 1)h n an

n
= a0 + r (r 1)( r 2 ) ...( r t + 1)h t at

t =1
with r = 0,..., n
Forming successive differences gives:

n
y r = yr+1 yr = ha1 + t{r (r 1)( r 2 )...( r t + 2 )}ht at

t =2
yr = yr +1 yr = 2h a2 + t (t 1){r ( r 1)(r 2 )...(r t + 3)}ht at

2
t =3
yr = yr +1 yr = 3 2 h a3 + t (t 1)( t 2 ){r( r 1)(r 2)...(r t + 4 )}ht a t

3
t =4
n 1 yr = ( n 1)!h n1a n1 + {n( n 1)( n 2)... 2}rhn an
n yr = n !h n an
From the above set of equations we may show that:
a0 = y 0
a1 = y0 h
a2 = 2 y0 ( 2! h 2 )
a3 = 3 y0 ( 3! h3 )
an = n y0 ( n! h n )
and hence the polynomial is p( x ) = y0 +

if we write x x0 = kh
x x0
( x x0 )( x x1 ) 2 y
y0 +
+ ...
h
h2
2!
where 0 k n
p( x ) = y0 + k y0 + k ( k 1)
then
2 y0
n y0
+ ... + k ( k 1)( k 2) ...( k n + 1)
2!
n!
By looking at the difference table we can see that this formula uses the values along the
diagonal of the differences of y - it is a FORWARD DIFFERENCE formula. It is therefore
used for interpolation near the beginning of a table where k is small.
Since the coefficients are binomial the formula may be written
n
y y0 + k Cr r y 0
r =1
Example (i)
(i) In the following table, use the Newton-Gregory Forward Interpolation formula to find
(a) f(2.4)
(b) f(8.7).
2
9.68
x
f(x)
Solution
4
10.96
6
12.32
8
13.76
10
15.28
Form a difference table and note that all differences > 2 are zero.
x
y=f(x)
9.68
10.96
12.32
13.76
10
15.28
2 y
1.28
0.08
1.36
0.08
1.44
0.08
1.52
(a)
x = 2.4; x = 2; h = 2; k = 0.2
we get f ( 2. 4) 9. 68 +
so
(b)
2. 4 2
(2. 4 2)(2. 4 4 ) 0. 08
1. 28 +
2
4
2
f ( 2. 4) 9. 68 + 0. 2 1. 28 + 0.1 ( 1. 6) 0. 04 = 9. 9296
x = 8.7; x = 2; h = 2; k = 3.35
we get f (8. 7) 9. 68 + 3. 35 1. 25 + 3. 35 2. 35 0. 04 = 14. 2829
Example (ii)
In the following table of ex use the Newton-Gregory formula of forward interpolation to
calculate
(a) e 0.12 ,
(b) e 2.00 .
x
ex
Solution:
0.1
1.1052
0.6
1.8221
1.1
3.0042
Form a difference table.
1.6
4.9530
2.1
8.1662
Note that in this case there is no difference column that is constant. This is to be expected
since ex cannot be represented by a polynomial function of finite degree.
x
y=e x
0.1
1.1052
2 y
3 y
4 y
0.7169
0.6
1.8221
0.4652
1.1821
1.1
3.0042
1.6
4.9530
2.1
8.1662
0.3015
0.7667
1.9488
0.1962
0.4977
1.2644
3.2132
(a)
With x = 0.12, xo = 0.1, h = 0.5, k = 0.04

0. 4652
2
0. 3015
0.1962
+0. 04 ( 0. 96) ( 1. 96)
+ 0. 04 ( 0. 96) ( 1. 96) ( 2. 96)
6
24
e0.12=1.1269
(b)
x = 2, xo=0.1, h = 0.5, k = 3.8
e 0.12 1.1052 + 0. 04 0. 7169 + 0. 04 ( 0. 96)
(correct value to 5 d.p. is 1.12750)
e 2 1.1052 + 2. 72422 + 2. 47486 + 0. 96239 + 0.12525

e 2 7. 3919 (to 4 dp).
(correct value 7.3891 to 4 dp)
In Example (i) the interpolation formula is identical with f(x), which is a quadratic function,
and the results for f(2.4) and f(8.7) will therefore be correct to the number of decimal places
retained.
In example (ii) the function ex is replaced by a 4th degree polynomial which takes the value of
ex at the five given entries. Because the successive difference decrease, higher differences are
relatively small and the value of the estimate converges. From direct calculation it turns out
that the error in the estimate for e 0.12 is about 0.05 percent and for e 2.00 it is about 0.04 percent.
It is important to note that interpolation may not always yield a valid result. Consider the
following example:
10
Example (iii)
Use the Newton-Gregory forward interpolation formula to estimate f(2.5) from the following
data:
2
1
x
f(x)
Solution:
3
2
4
6
5
24
6
120
Form the forward difference table:

x
y=f(x)
24
120
2 y
3 y
4 y
1
3
4
11
14
18
53
64
78
96
With x = 2.5, xo= 2, h = 1 then k = 0.5 and we obtain

f(2.5) 1 + 0.5 - 0.375 + 0.6875 - 2.0703 = -0.2578.
The function f(x) in this example happens to be the GAMMA FUNCTION written
( x ) = ( x 1)!
providing x is a positive integer.
In fact the value of f(2.5) is known to be 1.3293, and the value estimated by interpolation is
obviously incorrect. The reason for this is that the successive differences are increasingly
rapidly and the series obtained from the formula is divergent. The values of f(x) are increasing
so rapidly that it is not possible for a fourth-degree polynomial to represent the function
accurately even though it is exact at the five given values. In other words, there is not enough
information given - it is necessary to fix the value of the function at more closely spaced
intervals of x since it is varying so rapidly.
11
The Newton-Gregory Backward Interpolation Formula

The formula derived in the previous section is appropriate when the required value (xo,yo)
lies near the beginning of the tabulated data. When (xo,yo) is near the end of the table, it is
necessary to find a polynomial satisfying the (n+1) values (x-i,y-i), i=0,1,2,...,n. This may be
achieved in a parallel manner to the previous section. It follows that the backward formula is
2 y0
3 y0
n y0
+ k ( k + 1)( k + 2)
+ ... + k ( k + 1) ...( k + n 1)
2!
3!
n!
which may alternatively be written as
y = f ( x ) y0 + k y0 + k ( k + 1)
y y0 + k + r 1Cr r y0
k + r 1
where
Cr =
r =1
k ( k + 1)( k + 2 )...( k r + 1)
r!
The following table shows the paths of the forward and backward interpolation formulae
through the difference table for n = 5.
x 5
y5
x 4
y4
x 3
y3
Backward
y0
5
x 2
y2
y0
4
y0
3
x 1
y1
x0
y0
x1
y1
x2
y0
2
y0
y0
2 y0
3 y0
y2
x3
y3
x4
y4
x5
y5
4 y0
5 y 0
Forward
12
Example
Apply the backward formula to find e 2.00 in Example (ii) of the previous section.
x
y=e 2.00
0.1
1.1052
0.6
1.8221
1.1
3.0042
1.6
4.9530
2 y
3 y
4 y
0.7169
0.4652
1.1821
0.3015
0.7667
1.9488
0.1962
0.4977
1.2644
3.2132
2.1
8.1662
We have: x = 2.0, xo= 2.10, h = 0.5 and k = -0.2

e 2.00 8.1662 0. 2 3. 2132 0. 2 0. 8
1. 2644
0. 4977
0. 2 0.8 1. 8
+...
2
6
= 7.3920.
It has been observed that the two Newton-Gregory formulae use the information in the
difference table obtained by proceeding through the table along different paths. Both formulae
will give accurate results when the interval h in the given data is not too large (compared with
the variability in the driving function) and providing always that higher differences decrease in
value and eventually become negligible.
13
Test Exercises
1.
Form a difference table from the following data and deduce that f(x) is a
quadratic
function.
Use the table to find f(0).
Use the Newton-Gregory forward interpolation formula based on x = 0 to find
f(x).
x
f(x)
2
3.00
3
4.28
4
5.88
5
7.8
6
10.04
2.
Use an appropriate interpolation formula to estimate f(16.4) and f(23.5) from
following data table:
x
f(x)
16
261.3
18
293.7
20
330.0
14
22
372.2
24
422.3
the
Central differences
We have already defined the forward difference operator by
yr = y ( xr + h) y ( xr ) = yr +1 yr
and the backwards difference operator by

yr = y ( xr ) y ( xr h ) = yr yr 1
We now define the central difference operator, , by
h
h
yr = y xr + y xr = yr + 1 yr 1
2
2
2
2
The forwards and backwards difference operators may be written in terms of the central
difference operator
h h
h h
yr = y ( xr + h) y ( xr ) = y xr + + y xr +
2 2
2 2
= y xr + = yr + 1
2
Similarly it may be shown that yr = y xr = yr 1

2
2
The difference table, expressed in terms of central differences is as follows x
x
x0
y0
x1
y1
2 y
y 12
y 12
y1
3 y3
2
2 y2
y2
y 52
x3
3 y5
4 y2
2 y3
y3
y 72
x4
4 y
y 23
x2
3 y
3 y7
4 y3
y4
We may develop an interpolation formula in terms of central differences. We first define the
shift operator E by Ey r = yr +1
Now,
so
yr = yr + 1 yr
Ey r = (1 + ) yr
yr +1 = yr + yr = (1 + ) yr
E 1+
15
1
Also
yr = yr + 12 = E 2 yr
Note that
E 2 yr = E 2 ( yr+ 12 yr ) = E 2 yr + 12 E 2 yr = yr +1 yr
so
E 2 yr = yr = E 2 y
E 2
and
E 1+ 1+ E 2
2 2 + E 2 3
and
3 2 E 2 ( 2 + E 2 3 ) E 2 3 + E 4 E 2 3 + (1 + E 2 ) 4
3 E 2 3 + 4 + E 2 5
the operators E and commute
2 ( E 2 )( E 2 ) E 2 (1 + E 2 ) 2
1
We may substitute in the Newton-Gregory forward formula

y ( xr + kh ) = (1 + ) k yr
2
3
= 1 + k + k ( k 1) + k ( k 1)( k 2 ) +... yr
2!
3!
1
1
1
k ( k 1) 2
k ( k 1)( k 2 ) 12 3
= 1 + kE 2 +
( + E 2 3 ) +
( E + 4 + E 2 5 )+... yr
2!
3!
1
k ( k 1) 2 k ( k 1)( k 2) 12 3
= 1 + kE 2 +
+
E +... yr
2!
3!
We now consider the relationship between backward differences and central differences.
yr +1 = yr +1 yr
yr = yr +1 yr +1 = (1 ) yr +1
By definition
Ey r = yr +1 and yr = E 1 yr+1
Also
E 1 (1 )
1
Also by definition yr +1 = yr +1 yr = yr + 12 = E 2 yr +1
or alternatively
12
12
2 = ( E 2 )( E 2 ) = E 1 2 = (1 ) 2 = (1 E 2 ) 2
1
2 = 2 E 2 3
3 E 2 ( 2 E 2 3 ) = E 2 3 E 1 4 = E 2 3 (1 E 2 ) 4
1
3 E 2 3 4 + E 2 5
We may substitute in the Newton-Gregory backward formula
16
y ( xr + kh ) = (1 ) k yr
2
3
= 1 + k + k ( k + 1)
+ k ( k + 1)( k + 2) +... y r
2!
3!
k ( k + 1) 2
k ( k + 1)( k + 2) 12 3
1
1
1
( E 2 3 ) +
( E 4 + E 2 5 )+... yr
= 1 + kE 2 +
!
!
2
3
k ( k + 1) 2 k ( k + 1)( k + 2 ) 12 3
1
+
E +... yr
= 1 + kE 2 +
2!
3!
We now take the average of the forward and backward formulae and obtain
1
1
1
1
E 2 + E 2
k 2 2 k ( k 1)( k + 1) E 2 + E 2 3 k 2 ( k 2 1) 4
+... yr
y( xr + kh ) = 1 + k
+
+
+
2
2!
3!
2
4!
(The 4 term may be derived in the same way as the other terms were derived).
E
We introduce a final operator =
12
+ E2
, called the averaging operator.
2
1
The function of this averaging operator can be seen from the following:
E 2 + E 2
y + yr+1 yr 1 yr yr +1 yr 1
y r 1 ) =
( yr + 1 yr 1 ) = r
=
2
2
2
2
2
2
yr = ( yr+ 1
2
We have finally derived Stirling's formula for interpolating in terms of central differences.
Stirling's formula
k2
k ( k 1)( k + 1) 3 k 2 ( k 2 1) 4
y ( xr + kh) = 1 + k + 2 +
+
+... yr
2!
3!
4!
17
We may see how the terms relate to the difference table
x0
y0
x1
y1
xr 2
yr 2
xr 1
yr 1
2 y
2 y1
yr 32
2 yr 1
3 yr 1
2
4 yr
yr
yr
yr+ 1
3 y r+ 1
xr +1
4 y
y 12
yr 12
xr
3 y
yr + 1
yr +1
yr+ 3
2
xr + 2
yr +2
E 2 + E 2
1
1
1
1
Note that yr =
y r = 2 ( E 2 yr + E 2 yr ) = 2 (y r 12 + yr+ 12 )
2
ie yr is the average of the two differences highlighted in column 3 of the table.

Similarly 3 yr = 12 ( 3 yr 12 + 3 yr + 21 ) and is the mean of the two differences in boxes in column
5 of the table.
Stirling's formula is useful for interpolating near the middle of a table, as the following
example illustrates.
18
Example
Use Stirling's formula, with the difference table for ex to obtain e1.3 .
x
y=e x
0.1
1.1052
2 y
3 y
4 y
0.7169
0.6
1.8221
0.4652
1.1821
1.1
0.3015
3.0042
0.7667
1.9488
1.6
0.1962
0.4977
4.9530
1.2644
3.2132
2.1
8.1662
Solution
x = xr + kh
x = 1. 3, xr = 1.1, h = 0. 5
k=
x xr
= 0. 4
h
.
+ 1.9488 0.42
11821
e1.3 3.0042 + 0.4
0.7667
+
2
2
0.4( 0.4 1)( 0.4 + 1) 0.3015 + 0.4977 0.42 (0.4 2 + 1)
01962
.
+...
+
+
3!
2
4!
e1.3 3. 0042 + 0. 62618 + 0. 061336 0. 0223776 0. 0010987 = 3. 6682

(Exact result e1.3 = 3.6693, to 4 dp)
19
NUMERICAL DIFFERENTIATION
Differentiation is a very important mathematical process and a great deal of effort has been
devoted to the development of analytic techniques of finding the derivatives of various
mathematical functions. It often occurs, however, that it is not possible to utilise these
traditional methods. This happens when:
(i) The function is too complex
(ii) The function is unknown (when data are collected from some experiment).
In this section numerical techniques are described which provide an estimate of the derivative
of a tabulated function.
Imagine that you have a procedure which computes a function f(x) and you want to find its
derivative f'(x). Easy? The definition of the derivative is the limit, as h 0 of
f ( x)
f ( x + h) f ( x )
h
practically suggests an algorithm: Pick a small value h, evaluate f ( x ) and f ( x + h) and apply
the above formula. The above procedure, however, contains two main sources of error:
truncation error and roundoff error. The roundoff error depends upon the machine used to
carry out the computation - it is necessary to ensure that the value used for h is exactly
represented by the machine. Truncation error comes from the higher order terms in the Taylor
series expansion:
f ( x + h ) = f ( x ) + hf ( x ) + 12 h 2 f ( x ) + 16 h 3 f ( x ) +...
which leads to
f (x + h) f ( x)
= f ( x ) + 12 hf ( x ) +...
h
We will investigate methods of using difference tables to evaluate the derivatives of a function
numerically. From page 6, the Newton-Gregory forward polynomial is:
f ( xk ) Pn ( x k ) = f 0 + k f 0 + k C2 2 f 0 + .... + k Cn n f 0
Differentiating the above equation, remembering that
f ( xk ) Pn( xk ) =
1
h
{f
where k = ( x x0 ) h
dP dP dk
=
dx dk dx
+ 12 ( 2k 1)2 f 0 + 12 (3k 2 6k + 2 )3 f 0 +...
We can simplify this considerably if we take k = 0, giving a derivative corresponding to x = x0
f ( x0 ) 1h f 0 12 2 f0 + 13 3 f 0 14 4 f 0 + ... ( 1)n 1n n f0
20
(1)
Example
Use the data in the table below to estimate y'(1.7).
Use h = 0.2 and find the result using 1, 2, 3 and 4 terms of the formula.
x
y=e x
1.3
3.669
1.5
4.482
1.7
5.474
1.9
6.686
2 y
3 y
4 y
0.813
0.179
0.992
0.041
0.220
1.212
0.048
0.268
1.480
2.1
0.012
0.060
8.166
0.328
1.808
2.3
0.007
0.012
0.072
9.974
0.400
2.208
2.5
With one term
:
With two terms
:
With three terms :
With four terms
:
12.182
y (1. 7 ) = 01.2 (1. 212 ) = 6. 060

y (1. 7) = 01.2 (1. 212 12 0. 268) = 5. 390
y (1. 7 ) = 01.2 (1. 212 12 0. 268 + 13 0. 060) = 5. 490
y (1. 7 ) = 01.2 (1. 212 12 0. 268 + 13 0. 060 41 0. 012 ) = 5. 475
The data in the table are for y = ex, rounded to 3 dp.

The error in the derivative is rather small for 4 terms, which can be anticipated since the fourth
differences are approximately constant. This function is well-represented by a 4th degree
polynomial.
The above formula (1) for derivatives is a forward-difference approximation, for which the fit
provided by the interpolating polynomial is not symmetrical about x0 - interpretation is more
accurate near the centre of the range of fit.
For most purposes it is better to use central differences for the purpose of evaluating
derivatives. We use Stirling's formula as the basis of derivatives in an analogous manner to the
Newton-Gregory formula:
k2
k ( k 1)( k + 1) 3
k 2 ( k 2 1) 4
y ( xr ) = 1 + kyr + 2 yr +
yr +
yr +....
2!
3!
4!
k = ( x xr ) / h
and using the fact that
dy dy dk
=
dx dk dx
21
Differentiating gives y ( xr ) =
Using
1
3k 2 1 3
4k 3 2 k
2
+
+
+
y
4 yr +...
r
r
h
3!
4!
d 2 y d dy d dy dk d dy 1
=
=
=
dx 2 dx dx dk dx dx dk dx h
y ( xr ) =
we obtain
1 2
6k
12 k 2 2 4
3
+
+
y
y r +...
r
r
2
h
3!
4!
1 2
12k 2 2 4
3
+
+
y
k
y
yr
r
r
h2
4!
Summary of results
At a tabular point k = 0 and so
y ( xr ) =
1
1 3
1
5
yr yr + yr +...
h
6
30
y ( xr ) =
1 2
1
yr 4 y r +...
2
h
12
y ( xr )
1
3 yr
3
h
y iv ( xr )
1 4
yr
h4
Higher derivatives may be obtained in exactly the same way.

Writing the expressions out in full:
y ( x r )
1
1 yr 12 + yr + 12 1
y yr1
yr =
( y r yr1 + yr +1 yr ) = r +1
=
h
h
2
2h
2 h
y ( xr )
1 2
1
1
yr = 2 ( yr ) = 2 ( yr + 1 yr 1 )
2
2
2
h
h
h
=
1
y 2 yr + yr 1
yr +1 yr yr + yr 1 = r+ 1
2
h
h2
22
Example 1
Use the tabulated values of y = ex to approximate dy / dx and d 2 y / dx 2 at x = 1.1.
x
y=e x
0.1
1.1052 ( y0 )
2 y
3 y
4 y
0.7169
0.6
1.8221 ( y1 )
0.4652
1.1821
1.1
3.0042 ( y2 )
1.6
4.9530 ( y3)
2.1
8.1662 ( y4 )
0.3015
0.7667
1.9488
0.1962
0.4977
1.2644
3.2132
dy2 y3 y1 4. 9530 1. 8221
=
= 3.1309
dx
2h
1
(exact value 3.0042, to 4 d.p.)

d 2 y2
1
1
2 ( y3 2 y2 + y1 ) = 2 ( 49530
.
2 30042
.
+ 18221
.
) = 30668
.
2
dx
h
0.5
(exact value 3.0042, to 4 d.p.)
23
Test Exercises
1.
The following table is for y = 1 + log x. Determine estimates of y' at x = 0.15, 0.19 and
0.23 using
(a) one term, (b) two terms, (c) three terms of equation (1) at the beginning of this
section. Determine the errors by comparison with the analytical values.
x
y=(1+logx)
0.15
0.1761
0.17
0.2304
0.19
0.2788
0.21
0.3222
0.23
0.3617
0.25
0.3979
0.27
0.4314
0.29
0.4624
0.31
0.4914
2 y
3 y
2. A central-difference approximation for derivatives is more accurate than using a forwarddifference approximation. Repeat Exercise 1(b) but use appropriate central difference
formulae and compare the errors.
24
NUMERICAL SOLUTION OF DIFFERENTIAL EQUATIONS

Problems involving ordinary differential equations (ODEs) can always be reduced to the study
of sets of first order differential equations. For example, the second order equation:
d2 y
dy
+ q( x)
= r( x)
2
dx
dx
can be re-written as two first-order equations:
dy
= z( x)
dx
dz
= r ( x ) q ( x ) z( x )
dx
where z is a new variable. We shall therefore just consider first-order ODEs. Also, we shall
consider initial value problems in which all the y values are known at some starting x value.
The underlying idea of any method for solving the initial value problem is to rewrite the dy's
and dx's in the formula above as finite steps y and x and multiply the equations by x. This
gives algebraic formulas for the change in the functions when the independent variable x is
"stepped" by one "stepsize" x. In the limit of small stepsize, a good approximation to the
underlying differential equation is achieved. Literal implementation of this procedure results
in Euler's method, below, which is NOT recommended for practical use because it is
computationally inefficient. Euler's method is conceptually important, however, in that one
way or another all practical methods come down to the same idea: add small icrements to your
functions corresponding to derivatives (right-hand sides of the equations), multiplied by
stepsizes.
EULER's METHOD
Suppose we must solve the equation
dy
= f ( x , y ) , given that y = y0 when x = x0 .
dx
Divide the range of values of x into n equal steps x0 , x1 , x2 ,... , xn .

x x0
The step width is then given by: h = n
n
y
y1
x1
x0
y1 = y0 + h
dy
dx
= y0 + hf ( x0 , y0 )
x0
25
y2 = y1 + hf ( x1 , y1 )
yn = yn 1 + hf ( xn 1 , yn 1 )
dy
= y2 x
dx
Example 1
Take h=0.1 :
Find y(2.3) given y(2)=3.
y1 = y0 + 0.1( y0 2 x0 )
yi +1 = yi + 01
. y xi
2
i
x0 = 2 , y0 = 3
i
0
1
2
3
xi
2.0
2.1
2.2
2.3
yi
3.000
yi2
yi2-xi
0.1(yi2-xi)
yi+1
i
0
1
2
3
xi
2.0
2.1
2.2
2.3
yi
3.000
3.700
4.859
7.000
yi2
9.000
13.690
23.610
yi2-xi
7.000
11.590
21.410
0.1(yi2-xi)
0.700
1.159
2.141
yi+1
3.700
4.859
7.000
26
dy
= x+ y
dx
Example 2
Find y(0.05) correct to 5 d.p. given that y(0)=2
Take h=0.01:
i
0
1
2
3
4
5
xi
0.00
0.01
0.02
0.03
0.04
0.05
yi
2
2.02
2.0403
2.060903
2.08181203
2.1030301503
dy
= 2(1 + x ) y
dx
Example 3
y i+1 = yi + 0.01( x 0 + y0 ) , y0 = 2, x0 = 0
yi +1 = y i + 0. 01 yi '
yi'
2
2.03
2.0603
2.090903
2.12181203
yi+1
2.02
2.0403
2.060903
2.08181203
2.1030301503
Find y for x=2.0(0.2)3.0, given y(2)=5

yi +1 = yi + 0.2 y i '
Take h=0.2
i
0
1
2
3
4
5
0.01yi'
0.02
0.0203
0.020603
0.02090903
0.0212181203
xi
2
2.2
2.4
2.6
2.8
3.0
yi
5
5.2
5.44
5.712
6.0096
6.32768
27
yi'
1
1.2
1.36
1.488
1.5904
dy
= 1+ x y
dx
Example 4
Find y for x=1.0(0.2)2.0, given y(1)=2

yi +1 = yi + 0.2 y i '
Take h=0.2
i
0
1
2
3
4
5
xi
1.0
1.2
1.4
1.6
1.8
2.0
yi
2
2
2.04
2.112
2.2096
2.32768
yi'
0
0.2
0.36
0.488
0.5904
0.67232
EXERCISES
Solve the following ODEs to obtain the required values of y, under the given initial conditions:
1. If
dy
= x + y find y for x=0,(0.1),0.5, given that y(0)=1.
dx
(Solution: 1.0, 1.1, 1.22, 1.362, 1.5282, 1.72102)
2. If
dy
= 1 + xy find y for x=0,(0.1),0.5, given that y(0)=1.
dx
(Solution: 1.0, 1.1, 1.211, 1.3352, 1.4753, 1.6343)
3. Solve the above two equations by analytic methods and check your numerical results from
questions 1 and 2.
28
SECOND-ORDER ODEs
To solve first order ODEs, one boundary/initial condition was required; for second order ODEs
it becomes necessary to have two. We approximate the second derivative by the first term of the
central difference formula presented earlier in this booklet. We can see from the following
diagram how this can be derived geometrically:
C
B
A
y
-1
-1
Approximate slope at B may be given by:

(1)
BC:
dy y1 y0
dx
h
(Forward difference)
(2)
AB:
dy y0 y 1
dx
h
(Backward difference)
(3)
AC:
dy y1 y 1
dx
2h
(Central difference)
y1 y0 y0 y 1
d y
y 2 y0 + y 1
h
h
= 1
2
dx
h
h2
2
For the second derivative:
29
d2y
=y
dx 2
Example 1
Find y(0.6), given y(0)=1, y(0.1)=1.1

y i+1 2 yi + yi 1
= yi yi +1 = 2. 01 y i y i 1
h2
h=0.1
i
0
1
2
3
4
5
6
xi
0
0.1
0.2
0.3
0.4
0.5
0.6
d2y
= 4y
dx 2
Example 2
2.01yi-yi-1
1.211
1.33411
1.4705611
1.6217178
1.7890917
Find y(0.6), given y(0)=0.5, y(0.1)=6.0

y i+1 2 yi + yi 1
= yi yi +1 = 2. 01 y i y i 1
h2
h=0.1
i
0
1
2
3
4
5
6
yi
1
1.1
1.211
1.33411
1.4705611
1.6217178
1.7890917
xi
0
0.1
0.2
0.3
0.4
0.5
0.6
yi
0.5
6.0
11.56
17.2356
23.083556
29.162348
35.532762
30
2.01yi-yi-1
11.56
17.2356
23.083556
29.162348
35.532762
d 2 y dy
2y = 0
dx 2 dx
Example 3
h=0.01
i
0
1
2
3
4
5
Find y(0.05),given y(0)=3, y(0.01)=3.1
y i +1 2 yi + yi 1 y i+1 y i1
2 yi = 0
h2
2h
( yi+1 2 yi + yi1 ) ( yi+1 yi1 ) 2h 2 y i = 0
h
h
yi +1 (1 ) = 2 yi yi 1 + y i1 2 h2 y i
2
2
h
h
y i+1 = 2 2h 2 yi + 1 yi 1 1
y i+1 = 10050251
.
.
yi 0.995 yi 1 )
(19998
h
2
xi
0.00
0.01
0.01
0.03
0.04
0.05
yi
3
3.1
3.2305326
3.3928834
3.5886515
3.8197649
31
1.9998yi-0.995yi-1
3.21438
3.3759191
3.5707083
3.8006662
IMPROVED EULER METHOD (IEM)

The ordinary Euler method may be improved by estimating the slope over the interval between
points by the average of the slopes at the beginning and the end. We know the slope at the start
point, so we can use this to estimate the y-value of the end point. From this we can estimate the
slope at the end point. Taking the mean of the two slopes, we can then predict the true y-value
with greater accuracy.
y1
y0
x1
Estimate of y1 based on slope at x 0 :
y1 = y0 + hf ( x0 , y 0 )
Estimate of y1 based on slope at x1:
y1 = y0 + hf ( x1 , y1 )
Take y1 = 12 ( y1 + y1 ) = y0 + 12 h f ( x0 , y0 ) + f ( x1 , y1 )
This will set up an iterative scheme as follows:
y1(1) = y0 + hf ( x 0 , y0 )
{
}
h{ f ( x , y ) + f ( x , y )}
y1( 2 ) = y0 + 12 h f ( x0 , y0 ) + f ( x1 , y1(1) )
y1(3 ) = y 0 + 12
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
32
1( 2 )
Summary of the method

The steps in the method are as follows:
1.
We start with the given equation y = f ( x , y ), and with the initial condition that y = y 0 at
x = x 0 . We have to determine values of y for x = x 0 to xn step h .
2.
From the equation and the initial values we can determine y0 = f ( x 0 , y0 ) . This is the
ordinary Euler method.
3.
We now know x 0 , y0 , y 0 and h , and can therefore calculate

(a)
x1 = x0 + h
(b)
y1 = y0 + hy 0 (again, this is the same step as the ordinary Euler method)
(c)
We now calculate f ( x1 , y1 ) , and then we can estimate
y1 = y0 + 12 h y0 + f ( x1 , y1 )
(d)
4.
This is the first estimate of y1. We can, if we wish, repeat steps (c) and (d) in 3 above,
replacing y1 with the calculated value y1 and thereby finding an improved estimate of y1.
This is known as the recursive method. We can continue to employ this recursion until we
obtain enough accuracy (i.e. until the n th decimal place remains unchanged).
5.
Finally, we can calculate y1 = f ( x1 , y1 ) and repeat steps 3 and 4 to find x 2 and y 2 . This is
then repeated for all the values of y that we are required to find.
Example 1
Apply the Improved Euler Method, without recursion, to solve:

dy
= x+y
dx
for x=0.0 to 0.5 step 0.1, given that y(0)=1.
We proceed in a step-by-step routine: we have the equation y ' = x + y .

Initial conditions x0 = 0, y0 = 1
y0 = x0 + y0 = 1.
Also
h = 0.1
x 1 = x 0 + h = 0 .1
(a)
First calculate the value y1 = y0 + hy0 = 1 + 0.1 1 = 1.1
(b)
y1 = y0 + 12 h{ y0 + f ( x1, y1 )} = 1 + 0.05 {1 + ( 01
. + 11
. )} = 111
.
(c)
So we have x1 = 0.1, y1 = 1.1
(d)
We now take these last results in place of x0 , y0 and y0 in the above

calculations, and find the corresponding values for the next interval :
x2 , y2 and y2 . These values are shown in the following table:
y1 = 0.1 + 1.11 = 1. 21
33
xn
0.0
0.1
0.2
0.3
0.4
0.5
n
0
1
2
3
4
5
yn
1.0
1.21
1.44205
1.698465
1.981804
2.294893
yn
1.0
1.11
1.24205
1.398465
1.581804
1.794893
We now use recursion to obtain improved estimates of y.

As before, we have x0 = 0, x1 = 0.1, y0 = 1 and y0 = 1.
We next calculate, as before, y1 = y0 + hy0 = 1.1
and estimate y1 = y0 + 12 h{ y0 + f ( x1, y1 )} = 1 + 0.05 (1 + 12
. ) = 111
.
So far, this is all identical to the previous calculation, however, instead of stopping now and
using this value as the final estimate for y1 we substitute this value back into the last formula
above in place of y1 . For the value found so far we use the notation y1(1) .
We therefore proceed to find: y1( 2 ) = y0 + 12 h y0 + f ( x1 , y1(1) ) = 1 + 0.05 (1 + 01

. + y1(1) )
We now have a recursion relation which can be simplified by multiplying out the constants to
give:
y1( n) = 1. 055 + 0. 05 y1( n1)
and we can then tabulate the values of y1( n) for increasing values of n:
y1( n)
1.1
1.11
1.1105
1.110525
1.1105263
1.1105263
n
1
2
3
4
5
6
Note that after 5 iterations we have achieved 7 decimal place accuracy in the solution and that
the estimate without recursion is in error by 0.0005263, or 0.0474144%.
34
dy
= x2 y
dx
Example 1
x 0 = 1, x1 = 11
.,
y0 = 3,
f ( x, y ) = x2 y
y1(1) = y0 + hf ( x 0 , y0 )
1
y1( i +1) = y0 + h f ( x0 , y 0 ) + f x1 , y1( i )
2
)]
y1(1) = 3 + 0.1(3)
here
y1(i +1) = 3 + 0.05(3 + 11

. 2 y1( i ) ) = 315
. + 0. 0605 y1( i )
i
0
1
2
3
4
5
y1(i)
3
3.3
3.34965
3.35265
3.35283
3.35285
35
y1(i+1)=3.15+0.0605y1(i)
3.34965
3.35265
3.35283
3.35285
dy 1
=
dx y 2
Example 2
x0 = 1, x1 = 1.1, y0 = 1,
f ( x, y) =
Euler's Method:
1
. = 11
.
y1 = y0 + hf ( x0 , y0 ) = 1 + 01
1
Improved Euler Method:
y1(1) = y0 + hf ( x 0 , y0 )
1
y1( i +1) = y0 + h f ( x0 , y 0 ) + f x1 , y1( i )
2
1
y2
)]
y1(1) = 1.1
here
1
1
0. 05
y1( i +1) = 1 + 0.1(1 +
) = 1. 05 +
2
2
y1(i )
y1(i ) 2
i
0
1
2
3
4
Exact Solution:
dy 1
=
dx y 2
y1(i)
1
1.1
1.09132
1.09198
1.09193
y1(i+1)=1.05+0.05/y1(i)2
y 2 dy = dx
y3
= x+c
3
1.09132
1.09198
1.09193
y(1.1)1.091392
36
y = 3 3x 2
dy
= x+ y
x 0 = 0, x1 = 0.1, y0 = 1,
dx
y1(1) = y0 + hf ( x 0 , y0 )
1
y1( i +1) = y0 + h f ( x0 , y 0 ) + f x1 , y1( i )
2
Example 3
f (x, y) = x + y
)]
for y1
y1(1) = 1 + 0.1(1) = 11
.
y1(i +1) = 1 + 0.05(1 + 0.1 + y1( i ) ) = 1. 055 + 0. 05 y1( i )
here
i
0
1
2
3
4
5
y1(i)
1.1
1.11
1.1105
1.110525
1.1105263
1.1105263
y1(i+1)=1.055+0.05y1(i)
1.11
1.1105
1.110525
1.1105263
1.1105263
for y2:
y1 1.1105263
x1 = 0.1, x2 = 0. 2, y1 = 11105623
.
, f ( x, y ) = x + y
y1(1) = y0 + hf ( x 0 , y0 )
1
y1( i +1) = y0 + h f ( x0 , y 0 ) + f x1 , y1( i )
2
)]
f ( x1 , y1 ) = x1 + y1 = 1. 2105263
y2(1) = 11105263
.
+ 0.1(1. 2105263) = 1. 2315789
y2(i +1) = 11105263
.
+ 0.05(1.2105263 + 0. 2 + y2(i ) ) = 11810526
.
+ 0. 05 y 2( i )
here
i
0
1
2
3
4
5
y2(i)
1.2315789
1.2426315
1.2431842
1.2432118
1.2432132
1.2432133
y2(i+1)=1.1810526+0.05y2(i)
1.2426315
1.2431842
1.2432118
1.2432132
1.2432133
37
y2 1. 2432133
THE RUNGE-KUTTA METHOD

The Euler method, which has the formula
yn +1 = yn + hf ( xn , yn )
advances a solution from xn to xn +1 xn + h . The formula is non-symmetric - it advances a

solution through an interval h using derivative information only at the beginning of the interval,
as illustrated in the figure below:
3
2 l
1
x
x
There are several reasons why Euler's method is not recommended for practical use, among
them (i) it is not very accurate when compared with other methods run at the equivalent
stepsize and (ii) it is not very stable. Consider, however, the use of a step like in the above
illustration to take a "trial" step to the midpoint of the interval. Then use the value of both x
and y at that midpoint to compute the "real" step across the whole interval. The following
diagram illustrates the method:
4
1l
x
x
Mathematically, this can be written as
k1 = hf ( xn , yn )
k2 = hf ( xn + 12 h , yn + 21 k1 )
yn +1 = yn + k 2 + O(h 3 )
This procedure is called the second-order Runge-Kutta method (a method is conventionally

called n th order if its error term is O(h n +1 ) ).
We need not stop there. There are a number of ways of evaluating the right-hand-side f(x,y)
that all agree to first order, but that have different coefficients of higher-order error terms.
Adding up the right combination of these, we can eliminate the error terms order by order.
This is the basic idea of the Runge-Kutta method. By far the most commonly used of these
methods is the classical fourth-order Runge-Kutta formula:
38
Fourth-Order Runge-Kutta Formula

k1 = hf ( xn , yn )
k2 = hf ( xn + 12 h , yn + 21 k1 )
k3 = hf ( xn + 12 h , yn + 12 k2 )
k4 = hf ( xn + h , yn + k 3 )
yn +1 = yn + 16 k1 + 13 k 2 + 13 k 2 + 16 k 3 + O ( h5 )
The fourth-order Runge-Kutta method required four evaluations of the right hand side per step
h . This will be superior to the second-order method if at least twice as large a step is
possible for the same accuracy. Is this so? The answer is: often, perhaps even usually, but not
always. This emphasises the fact that high order does not always mean high accuracy. A
good system should exert some adaptive control over its own progress, making frequent
changes to its stepsize to maintain the error within predetermined limits. This, however, is
beyond the scope of this unit.
Example 1
Solve
dy
= x + y,
dx
given that y( 0) = 1
x0 = 0, y0 = 1, y0 = 1, h = 0.1
x1 = x0 + h = 0.1
k1 = hy0 = 0.1
k2 = hf ( 0 + 0. 05, 1 + 0. 05) = hf (0. 05, 1. 05)
k3 = hf ( 0 + 0. 05, 1 + 0. 055) = hf (0. 05, 1. 055)
k4 = 0.12105
y1 = 1.110342
39
k 2 = 0.11
k 3 = 0.1105
Exercises
1.
Using the 4th order Runge Kutta method,

(a)
compute an estimate of x(0.75) for the initial value problem:

dx
= x + t + xt ,
x ( 0) = 1
dt
using a step size of h = 0.15 .
(b)
compute an estimate of x(2) for the initial value problem

dx
1
,
=
x(1) = 2
dt x + t
using a step size of h = 0.1 .
(Solution: 3.2345)
(Solution: 2.2771)
2.
Solve y = x + y with y = 1 at x = 0.
Determine the function values for y for x=0(0.1)0.5.
(Solution: 1.0,1.110342,1.242806,1.399718,1.583649,1.797442)
3.
Solve y = y 2 xy with y = 0. 4 at x = 0.
(Solution: 0.4, 0.4259, 0.4374, 0.4319, 0.4085, 0.3689)
4.
yx
with y = 1 at x = 0.
y+x
Solve y =
(Solution: 1.0, 1.1679, 1.2902, 1.3817, 1.4497, 1.4983)
40

Newton Forward and Backward

Uploaded by

Copyright:

Available Formats

Newton Forward and Backward

Uploaded by

Document Information

Original Description:

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Newton Forward and Backward

Uploaded by

Copyright:

Available Formats

BSc Computing Mathematics

BSc Computing Maths with Business Studies

We then find the point on that line with an x-coordinate of x .

Therefore, to find the required value of y we insert x for x:

Consider a simple example:

(6.5 6)(49 36) = 36 + 05. 13 = 42.5

The correct value of 6.5 is 42.25 the error ET = -0.25.

A General Approach to Interpolation

is called the first difference of yr , denoted by yr :

Similarly, the relationship yr +1 yr is the second difference:

The third differences are

Errors in a Difference Table

The Newton-Gregory Interpolation Formula

If the equal intervals between successive values is h, then xr xs = ( r s) h .

yr = a0 + rha1 + r (r 1)h 2 a2 + ... + r (r 1)(r 2 )...(r n + 1)h n an

= a0 + r (r 1)( r 2 ) ...( r t + 1)h t at

Forming successive differences gives:

y r = yr+1 yr = ha1 + t{r (r 1)( r 2 )...( r t + 2 )}ht at

yr = yr +1 yr = 2h a2 + t (t 1){r ( r 1)(r 2 )...(r t + 3)}ht at

yr = yr +1 yr = 3 2 h a3 + t (t 1)( t 2 ){r( r 1)(r 2)...(r t + 4 )}ht a t

n 1 yr = ( n 1)!h n1a n1 + {n( n 1)( n 2)... 2}rhn an

and hence the polynomial is p( x ) = y0 +

Form a difference table.

With x = 0.12, xo = 0.1, h = 0.5, k = 0.04

x = 2, xo=0.1, h = 0.5, k = 3.8

e 0.12 1.1052 + 0. 04 0. 7169 + 0. 04 ( 0. 96)

(correct value to 5 d.p. is 1.12750)

e 2 1.1052 + 2. 72422 + 2. 47486 + 0. 96239 + 0.12525

(correct value 7.3891 to 4 dp)

Form the forward difference table:

With x = 2.5, xo= 2, h = 1 then k = 0.5 and we obtain

providing x is a positive integer.

The Newton-Gregory Backward Interpolation Formula

We have: x = 2.0, xo= 2.10, h = 0.5 and k = -0.2

and the backwards difference operator by

We now define the central difference operator, , by

Similarly it may be shown that yr = y xr = yr 1

the operators E and commute

We may substitute in the Newton-Gregory forward formula

We may substitute in the Newton-Gregory backward formula

We introduce a final operator =

We may see how the terms relate to the difference table

ie yr is the average of the two differences highlighted in column 3 of the table.

e1.3 3. 0042 + 0. 62618 + 0. 061336 0. 0223776 0. 0010987 = 3. 6682

Differentiating the above equation, remembering that

+ 12 ( 2k 1)2 f 0 + 12 (3k 2 6k + 2 )3 f 0 +...

We can simplify this considerably if we take k = 0, giving a derivative corresponding to x = x0

y (1. 7 ) = 01.2 (1. 212 ) = 6. 060

The data in the table are for y = ex, rounded to 3 dp.

and using the fact that

At a tabular point k = 0 and so

Higher derivatives may be obtained in exactly the same way.

dy2 y3 y1 4. 9530 1. 8221

(exact value 3.0042, to 4 d.p.)

NUMERICAL SOLUTION OF DIFFERENTIAL EQUATIONS

Divide the range of values of x into n equal steps x0 , x1 , x2 ,... , xn .

Find y(2.3) given y(2)=3.

Find y(0.05) correct to 5 d.p. given that y(0)=2

Find y for x=2.0(0.2)3.0, given y(2)=5