International Journal of Analysis and Applications
ISSN 2291-8639
Volume 3, Number 2 (2013), 104-111
http://www.etamaths.com

NUMERICAL DIFFERENTIATION AND INTEGRATION

THROUGH AITKEN-NEVILLE SCHEMES

RAMESH KUMAR MUTHUMALAI

Abstract. Some new formulas are given to approximate higher order deriva-
tives and integrals through Aitken-Neville iterative schemes for arbitrary s-

paced grids. An algorithm is given in MATLAB for numerical differentiation.
Also, numerical examples are provided to study error analysis of new formulas
for numerical differentiation and integration.

1. Introduction.

The problem of numerical differentiation is a long-standing issue. There are
plenty of published works devoted to estimate derivatives of a function numerically
for arbitrary spaced grids. Some of them are polynomial interpolation type [5, 9],
finite difference formulas [2, 4, 6] and Richardson extrapolation [2]. Most of them do
not aim to iterate derivatives upto kth order per addition of a new grid. The finite
difference formulas for the calculation of any order derivative in a one dimensional
grid with arbitrary spacing are discussed in Ref[4] (algorithm given in Fortran)
with the cost of O(kn2) operations. It also generates sequence of approximate
derivatives upto order k per addition of a new grid. It should be noted that there
exists a great deal of formulas and techniques of numerical integration. But, they
have been of considerable complexity and often been limited to lower order formulas
on equidistantly spaced grids.

Aitken-Neville schemes [5, 7, 8] are popular interpolation methods to iterate
interpolation when a new grid is added. An obvious advantage of these schemes
is that it gives a good idea of the accuracy of the result at any stage [8]. In the
present study we describe new formulas to approximate numerical derivatives and
integrals through Aitken-Neville iterative schemes for arbitrary spaced grids. Also,
we provide an efficient algorithm in MATLAB to iterate numerical differentiation
(upto kth order) per addition of a new grid with cost of O(k2n) operations.

2. Numerical differentiation

2.1. Formulae for numerical differentiation.

2010 Mathematics Subject Classification. 65D25; 65D30.
Key words and phrases. Aitken-Neville scheme; iterative method; numerical differentiation;

numerical integration.

c⃝2013 Authors retain the copyrights of their papers, and all
open access articles are distributed under the terms of the Creative Commons Attribution License.

104



NUMERICAL DIFFERENTIATION AND INTEGRATION 105

Definition 2.1. Define N
(r)
j [x] =

1
(xj−x)r

and

(2.1) N
(r)
j,j+1,...,j+i[x] =

1

xi+j − xj

∣∣∣∣∣ N
(r)
j,j+1,...,j+i−1[x] xj − x

N
(r)
j+1,j+2,...,j+i[x] xi+j − x

∣∣∣∣∣ .
Also, define Ñ

(r)
j [x] =

f(xj)
(xj−x)r

and

(2.2) Ñ
(r)
j,j+1,...,j+i[x] =

1

xi+j − xj

∣∣∣∣∣ Ñ
(r)
j,j+1,...,j+i−1[x] xj − x

Ñ
(r)
j+1,j+2,...,j+i[x] xi+j − x

∣∣∣∣∣ .
Where N

(r)
j,j+1,j+2,...,j+i[x] and Ñ

(r)
j,j+1,j+2,...,j+i[x] are constructed by Neville scheme

of interpolation, i = 1, 2, . . . , n and j = 0, 1, 2, 3 . . . n − i. At the nth iteration,
assume that N

(r)
0,1,2,...,n[x] = N

(r)(x), r = 0, 1, 2, . . . k.

Definition 2.2. Define A
(r)
0,j[x] =

1
xj−x0

∣∣∣∣∣
1

(x0−x)r
x0 − x

1
(xj−x)r

xj − x

∣∣∣∣∣ and
(2.3) A

(r)
0,1,2,...,i−1,i,j[x] =

1

xj − x0

∣∣∣∣∣ A
(r)
0,1,...,i−1,i[x] xi − x

A
(r)
0,1,...,i−1,j[x] xj − x

∣∣∣∣∣ .
Also, define Ã

(r)
0,j[x] =

1
xj−x0

∣∣∣∣∣
f(x0)

(x0−x)r
x0 − x

f(xj)
(xj−x)r

xj − x

∣∣∣∣∣ and
(2.4) Ã

(r)
0,1,2,...,i,j[x] =

1

xj − xi

∣∣∣∣∣ Ã
(r)
0,1,...,i−1,i[x] xi − x

Ã
(r)
0,1,...,i−1,j[x] xj − x

∣∣∣∣∣ .
Where A

(r)
0,1,2,...,i−1,i,j[x] and Ã

(r)
0,1,2,...,i−1,i,j[x] are constructed by Aitken’s scheme

of interpolation, i = 1, 2, . . . , n and j = i + 1, i + 2, . . . , n. At the nth iteration,

assume that A
(r)
0,1,2,...,n[x] = A

(r)(x), r = 1, 2, . . . k.

The following theorem gives recursive formulas for numerical differentiation through
Aitken-Neville schemes.

Theorem 2.3. Let x, x0, x1, x2, . . . xn are n + 1 distinct numbers on the interval
[a, b], k ∈ W and f ∈ Cn+k+1[a, b]. Then

(2.5)
k∑

i=0

f(i)(x)

i!
N(k−i)(x) = Ñ

(k)
0,1,...,n[x] + E(x).

and

(2.6)
k∑

i=0

f(i)(x)

i!
A(k−i)(x) = Ã

(k)
0,1,...,n[x] + E(x).

Where E(x) =
f(n+k+1)(ξ)
(n+k+1)!

∏n
i=0(x−xi) and for min{x, x0, . . . , xn} < ξ < max{x, x0, . . . , xn}.

Proof. Let Pn(x) is a polynomial of degree ≤ n in x that approximates the function
f. Using equation (2.7), gives

f[ x, . . . , x︸ ︷︷ ︸
k + 1 times

] = D0,1,2,...,n[x] + l(x)
f(n+k+1)(ξ)

(n + k + 1)!
.



106 RAMESH KUMAR MUTHUMALAI

Where Pn[ x, . . . , x︸ ︷︷ ︸
k + 1 times

/x0, x1, . . . xn] = D0,1,2,...,n[x] and l(x) =
∏n

i=0(x − xi) and

min{x, x0, . . . , xn} < ξ < max{x, x0, . . . , xn}. Expand Pn[x, . . . , x︸ ︷︷ ︸
k times

, xj] repeatedly,

using the recursive formula of divided difference [2, p-41] to get

Pn[x, . . . , x︸ ︷︷ ︸
k times

, xj] = −
k−1∑
r=0

P
(r)
n (x)

r!

1

(xj − x)k−r
+

Pn(xj)

0!(xj − x)k
.(2.7)

Applying equation (2.7) for two consecutive data xj, xj+1

Pn[ x, . . . , x︸ ︷︷ ︸
k + 1 times

/xj, xj+1] =
1

xj+1 − xj

∣∣∣∣∣∣∣∣
Pn[x, . . . , x︸ ︷︷ ︸

k times

/xj] xj − x

Pn[x, . . . , x︸ ︷︷ ︸
k times

/xj+1] xj+1 − x

∣∣∣∣∣∣∣∣ .
Using (2.7) and properties of determinants, finds that

Pn[ x, . . . , x︸ ︷︷ ︸
k + 1 times

/xj, xj+1] = −
k−1∑
r=0

P
(r)
n (x)

r!
N

(k−r)
j,j+1 [x] + Ñ

(k)
j,j+1[x].

Proceeding this for some i,

Pn[ x, . . . , x︸ ︷︷ ︸
k + 1 times

/xj, xj+1, . . . , xj+i]

= −
k−1∑
r=0

P
(r)
n (x)

r!
N

(k−r)
j,j+1,...,j+i[x] + Ñ

(k)
j,j+1,...,j+i[x].

Putting i = n and j = 0,(i.e at nth iteration)

Pn[ x, . . . , x︸ ︷︷ ︸
k + 1 times

/x0, x1, . . . , xn] = −
k−1∑
r=0

P
(r)
n (x)

r!
N

(k−r)
0,1,...,n[x] + Ñ

(k)
0,1,...,n[x].

Since Pn approximates f and N
(r)
0,1,...n[x] = N

(r)(x). Then

f[ x, . . . , x︸ ︷︷ ︸
k + 1 times

] − l(x)
f(n+k+1)(ξ)

(n + k + 1)!
= −

k−1∑
r=0

f(r)(x)

r!
N(k−r)(x) + Ñ

(k)
0,1,...,n[x].

Since f[ x, . . . , x︸ ︷︷ ︸
k + 1 times

] =
f(k)(x)

r!
and after simplification gives (2.5). Similar manner

using Definition 2.2, yields (2.6). �

Theorem 2.4. Let x, x0, x1, x2, . . . xn are n + 1 distinct numbers on the interval
[a, b], k ∈ W and f ∈ Cn+k+1[a, b]. Then direct formula for numerical differentia-
tion through Neville scheme as follows

(2.8)
f(t)(x)

t!
= atf(x)χ +

t∑
k=χ

at−kÑ
(k)
0,1,2,...n[x] + ED(a; x).



NUMERICAL DIFFERENTIATION AND INTEGRATION 107

and

a0 = 1, at = −
t−1∑
k=0

akN
(t−k)(x), t = 1, 2, 3, . . . .(2.9)

The direct formula for numerical differentiation through Aitken scheme as follows

(2.10)
f(t)(x)

t!
= btf(x)χ +

t∑
k=χ

bt−kÃ
(k)
0,1,2,...n[x] + ED(b; x).

and

b0 = 1, bt = −
t−1∑
k=0

bkA
(t−k)(x), t = 1, 2, 3, . . . .(2.11)

Where

χ =

{
1 if x = xi
0 if x ̸= xi,

ED(y; x) = l(x)
t∑

m=χ

yt−m
f(n+m+1)(ξm)

(n + m + 1)!
.

and y = [y0 y1 . . . yt] is an one dimensional array of length t + 1.

Proof. If one use equation (2.5) recursively to derive direct formulas for kth order
differentiation, then we obtain the following form, for some t = 0, 1, 2, 3 . . .

f(t)(x)

t!
= atf(x)χ +

t∑
k=χ

at−kÑ
(k)
0,1,2,...n[x] + ED(a; x).

To evaluate these unknown a′s, set f(x) = 1 in the above equation, gives a0 = 1
when t = 0 and for t = 1, 2, 3, . . . gives

t−1∑
k=1

at−kN
(k)(x) + at = 0.

This gives (2.9). Similarly, from (2.6), it can be easily find (2.10). �

2.2. Numerical experiment. In this subsection, we present algorithm (coded in
MATLAB) and numerical result to illustrate the performance of the new formulas
given in (2.9) and (2.10) for arbitrary spaced grids.

Algorithm 2.5. Numerical differentiation through Neville scheme
function [d]=differentiation(x,y,k,s)
% Input parameters
% s location where approximations are to be accurate
% x(1:n) grid point locations, found in x(1:n)
% y(1:n) functional value locations, found in y(1:n)
% k highest derivatives are sought at ’s’
% Output parameters
% d(1:k+1,1:n) sequence of approximate derivatives of order 0:k
n=length(x);
% Check whether ’s’ coincide with any one of ’x’. (i.e) functional value
% of ’s’ is known or not.



108 RAMESH KUMAR MUTHUMALAI

xs=x-s*ones(1,n);
pos=find(˜(xs&xs));
chi=˜all(xs);
% If chi=1 (functional value at ’s’ is known), then swap respective
% position of s and its functional value with first element of x and y.
if chi==1;
x([1 pos])=x([pos 1]);
y([1 pos])=y([pos 1]);
nevy(1,1:n)=y(1);
nev(1,1:n)=1;
end
f=1+chi;

% Determine the value of N(r)(s) andÑ
(r)
0,1,2,...,n[s] by Neville scheme and store %it in

’nev’ and ’nevy’ respectively
for ij=f:k+1

yy(f:n)=1./(x(f:n)-s).ˆ(ij-1);
yy1=y.*yy;
nev(ij,f)=yy(f);
nevy(ij,f)=yy1(f);
for i= f:n

yy(f:n+f-i-1)=(yy(f+1:n+f-i).*(s-x(f:n+f-i-1))-yy (f:n+f-i-1)
.*(s-x(i+1:n)))./(x(i+1:n)-x(f:n+f-i-1));

yy1(f:n+f-i-1)=(yy1(f+1:n+f-i).*(s-x(f:n+f-i-1))- yy1(f:n+f-i-1)
.*(s-x(i+1:n)))./(x(i+1:n)-x(f:n+f-i-1));

nev(ij,i)=yy(f);
nevy(ij,i)=yy1(f);

end
end
% Find the value of a’s and store them in a two dimensional array ’b’
for m=1:n

b(1,m)=1;
mn=min(m,k+1);
a(1,1)=1;
FACT(mn)=1;
for i=2:mn

b(i,m)=0;
for j=1:i-1
b(i,m)=b(i,m)-b(j,m)*nev(i-j+1,m);
end
a(i,i:-1:1)=b(1:i,m)’;
FACT(mn+1-i)=factorial(i-1);

end
% Evaluate sequence of approximate derivatives up to order ’k’

Diff=(a(mn:-1:1,:)*nevy(1:mn,m)).*FACT’;
d(1:mn,m)=Diff(mn:-1:1);

end



NUMERICAL DIFFERENTIATION AND INTEGRATION 109

The algorithm for numerical differentiation through Aitken scheme is similar to
the algorithm 2.5. The following are some notes regarding the implementation of
the algorithm:

• A call to differentiation to obtain value of mth derivative returns also value
of kth derivative, k = 0, 1, 2, . . . , m with no additional costs.

• The code returns all the data above also for points which extend only over
x0, x1, . . . , xj, j = 0, 1, 2, . . . , n, still no additional costs.

• It requires O
(
k2n

)
operations. Also, the maximum size of array used is

(k + 1) × n.
In the following example, we have taken 20 Chebyshev first kind of points on [-1,1]
of function f(x) = ex. We evaluate the first order derivatives at 100 equally spaced
points on [−1 + h, 1 − h], where h = 2/101. Figure 1 plots the errors for numerical

−1 −0.5 0 0.5 1
10

−20

10
−15

10
−10

10
−5

Neville
Aitken

Figure 1. Relative errors in computed f′(x) for 20 Chebyshev
first kind of points of f(x) = ex.

differentiation through Aitken-Neville schemes. We see that the Neville scheme
performs stably and Aitken scheme becomes very unstable as x approaches one end
of the interval.

3. Numerical Integration.

3.1. Formulae for Numerical integration. In this section, we derive numeri-
cal integration formulas for arbitrary spaced grids to any order of accuracy. Let
x, x0, x1, x2, . . . , xn are distinct numbers on the closed interval [a, b] and f ∈ C(n+1)[a, b].
To derive an expression for the definite integral

∫ x+h
x

f(x)dx , (h ̸= 0) expanding
through Taylor series, yields

(3.1)

∫ x+h
x

f(x)dx =
n∑

k=0

f(k)(x)

(k + 1)!
hk+1 + O

(
hn+2

)
.



110 RAMESH KUMAR MUTHUMALAI

Using (2.8), gives∫ x+h
x

f(x)dx =
n∑

k=0

hk+1

k + 1

(
akf(x)χ +

k∑
m=χ

ak−mÑ
(m)
0,1,2,...n[x]

+l(x)

k∑
m=χ

ak−m
f(n+m+1)(ξm)

(n + m + 1)!

)
+ O

(
hn+2

)
.

Setting γk =
∑n−k

m=0 am
hm+k+1

m+k+1
, k = 0, 1, 2, . . . n and rearranging above equation,

yields

(3.2)

∫ x+h
x

f(x)dx = γ0f(x)χ +

n∑
k=χ

γkÑ
(k)
0,1,2,...,n[x] + EI(γ; x).

Where

(3.3) EI(γ; x) = l(x)
n∑

m=χ

γm
f(n+1+m)(ξm)

(n + 1 + m)!
+ O

(
hn+2

)
.

Similarly, numerical integration formula through Aitken scheme can be easily found
from (2.1) as follows

(3.4)

∫ x+h
x

f(x)dx = γ′0f(x)χ +
n∑

k=χ

γ′kÃ
(k)
0,1,2,...,n[x] + EI(γ

′; x).

Where γ′k =
∑n−k

m=0 bm
hm+k+1

m+k+1
, k = 0, 1, 2, . . . n

3.2. Numerical Experiment. We evaluate the integral
∫ 1
−1 e

xdx on Chebyshev

and other point distributions. The exact value of the integral is e − e−1. Table

It Equally spaced points Chebyshev first kind Chebyshev second kind
Aitken Neville Aitken Neville Aitken Neville

1 6.8696e-01 6.8696e-01 1.2962e+00 1.2962e+00 6.8696e-01 6.8696e-01
2 3.4633e-01 3.4633e-01 9.1875e-01 9.1875e-01 2.9095e-01 2.9095e-01
3 1.3011e-01 1.3011e-01 4.2473e-01 4.2473e-01 9.8765e-02 9.8765e-02
4 3.6744e-02 3.6744e-02 1.3578e-01 1.3578e-01 2.4775e-02 2.4775e-02
5 7.8565e-03 7.8565e-03 3.0515e-02 3.0515e-02 4.4655e-03 4.4655e-03
6 1.2724e-03 1.2724e-03 4.7856e-03 4.7856e-03 5.5864e-04 5.5864e-04
7 1.5482e-04 1.5484e-04 5.1025e-04 5.1025e-04 4.6114e-05 4.6114e-05
8 1.3890e-05 1.3899e-05 3.5295e-05 3.5295e-05 2.3236e-06 2.3236e-06
9 1.2211e-05 7.2372e-07 1.4669e-06 1.4670e-06 6.3110e-08 6.3209e-08
10 3.3453e-04 1.0280e-06 2.1951e-09 3.7429e-08 4.2184e-08 4.1447e-09
11 2.2040e-02 2.9403e-05 1.9287e-05 4.6798e-06 2.5779e-05 1.1667e-04
12 1.0882e+00 5.4514e-04 1.9782e+00 1.0000e+00 1.5238e+03 2.1884e+02
13 2.9428e+01 6.3104e-03 1.8712e+12 4.6780e+12 1.1153e+05 8.6310e+08

Table 1. Comparisons between relative errors of computation of
∫
1
−1 e

xdx on various

point distributions through Aitken-Neville schemes

1 shows the comparisons between relative errors of computation of the required
integral through (3.2) and (3.4) on 13 grids of various point distributions (equally
spaced and Chebyshev grids). We observe that the relative errors of both schemes
converge up to 9th iteration on equally spaced grids and converge up to 10th iteration
for Chebyshev grids, after that, diverge rapidly in all point distributions. Also, it
gives an idea of the accuracy of the result per addition of a new grid. It is observed
that the new formulas give greater accuracy on Chebyshev point distributions than
evenly spaced grids.



NUMERICAL DIFFERENTIATION AND INTEGRATION 111

4. Conclusion.

The formulas for approximating derivatives and integrals for arbitrary spaced
grids through Aitken-Neville schemes have been developed in this article. We have
provided an efficient algorithm (in MATLAB) to iterate numerical differentiation
with cost of O(k2n) operations. Also, the numerical result shows that the numerical
differentiation through Neville scheme is stable than Aitken scheme. All available
numerical integration formulas do not iterate integral values per addition of a new
grid to any degree of accuracy. But the new formulas for integration described in
section 3, will do this with cost of O(n3) operations. Also, an interesting advantage
is that they are also applicable to unevenly grids.

References

[1] K. E. Atkinson, An introdction to numerical Analysis, 2 Ed, John Wiley & Sons, Newyork
(1989).

[2] S.D. Conte, Carl de boor, Elementary numerical Analysis, 3 Ed, McGraw-Hill, Newyork,
USA (1980).

[3] M. Dvornikov, Formulae for Numerical differentiation, JCAAM, 5 (2007), 77-88[e-print arx-
iv:math.NA/0306092].

[4] B. Fornberg, Calculation of weights in finite difference formulas, SIAM Rev, Vol 40, No 3,
685-691 (1998).

[5] F.B. Hildebrand, Introduction to Numerical analysis, 2 Ed, McGraw-Hill, Newyork (1974).
[6] J.Li, General Explicit difference formulas for Numerical differentiation, J.Comp & Appl.

Math, 183, 29-52 (2005).
[7] G. M. Phillips, Interpolation and approximation by polynomials, Springer-Verlag, Newyork

(2003).
[8] S.S. Sastry, Introdutory methods of Numerical Analysis, 4 Ed, Prentice hall of India, New

Delhi (2005).
[9] E. Süli & D. Mayers, An introduction to Numerical Analysis, Cambridge University Press,

UK (2003).

Department of Mathematics, D.G. Vaishnav College, Arumbakkam, Chennai-106, Tamil

Nadu, India