Dual Algorithms
F. Peñuñuri∗a , O. Carventea , M. A. Zambrano-Arjonaa , and Carlos A. Cruz-Villarb

arXiv:1411.2684v3 [cs.CE] 10 Jan 2017

a

Facultad de Ingenierı́a, Universidad Autónoma de Yucatán, A.P. 150, Cordemex, Mérida,
Yucatán, México.
b
Cinvestav-IPN, Departamento de Ingenierı́a Eléctrica, Av. IPN 2508, A. P. 14-740,
07300, México D.F., México.

Abstract
The cubic spline interpolation method, the Runge–Kutta method, and the Newton–Raphson
method are extended to dual versions (developed in the context of dual numbers). This extension allows the calculation of the derivatives of complicated compositions of functions which
are not necessarily defined by a closed form expression. The code for the algorithms has
been written in Fortran and some examples are presented. Among them, we use the dual
Newton–Raphson method to obtain the derivatives of the output angle in the RRRCR spatial
mechanism; we use the dual normal cubic spline interpolation algorithm to obtain the thermal
diffusivity using photothermal techniques; and we use the dual Runge–Kutta method to obtain
the derivatives of functions depending on the solution of the Duffing equation.
Keywords— Dual numbers, Derivatives, Runge–Kutta algorithm, Newton–Raphson algorithm Cubic spline interpolation.

1

Introduction

Analogous to a complex number z = a + i b where a and b are real numbers and i2 = −1, a dual
number is defined as r̂ = a + ǫ b with a and b real numbers and ǫ2 = 0. Such numbers were
introduced by Clifford who also developed their algebra in the late nineteenth century [1].
The fact that ǫ2 = 0 suggests that the dual numbers can be used to differentiate functions,
since in analogy to an infinitesimal dx, quantities of order dxn with n an integer greater than or
equal to two are usually neglected. It turns out that this reasoning is correct. This can be easily
proved for polynomials and then via the Taylor Series, generalized for any analytic function. So,
extending a real function to a dual function one can numerically obtain its derivatives. Nonetheless, most of the applications of dual numbers are in the area of mechanics (see for example [2]
where some contributions to mechanics based on dual numbers are presented) and only relatively
recently have they been used to obtain the derivatives of functions [3–6]. Moreover, there are
no papers addressing the dualization (from now on the term dualization will refer to extending a
function or algorithm to the context of dual numbers) of algorithms such as the Newton–Raphson
algorithm, the Runge–Kutta algorithm, or the cubic spline interpolation method. The present
∗

Corresponding Author. Tel: + 52 999 9300550x1051, 1052
E-mail addresses: francisco.pa@correo.uady.mx (F. Peñuñuri), osvaldo.carvente@correo.uady.mx (O. Carvente),
miguel.zambrano@correo.uady.mx (M. A. Zambrano-Arjona), cacruz@cinvestav.mx (Carlos A. Cruz-Villar).

1

paper shows that a dualization of these algorithms allows the numerical calculation of the derivatives of complicated compositions of functions which frequently arise in science and engineering
applications.
The problem of numerical differentiation has been addressed by many researchers. A complete
review of the literature on this topic is beyond the scope of this paper. Nevertheless we want
to cite some works which represent most of the techniques used to obtain derivatives numerically
[5–22]. It is worthwhile to mention that if there is a closed form expression of the function to be
differentiated, the automatic differentiation methods (AD) [6, 9, 19] are (from the point of view of
precision, accuracy, and efficiency) the best choice to calculate the derivatives. In the case when
the derivatives of a set of data are required, two excellent approaches are presented in [13, 22].
However, in some applications, one deals with derivatives of the composition of functions defined
by an expression rather than an explicit function or set of data. For example, one could be
interested in calculating the derivatives of the composition of a function with another that is the
solution of some differential equation, or in calculating the derivatives of a function which depends
on another function which is the spline interpolation of some data. In such cases, the use of dual
numbers is especially suited for calculating the derivatives. As illustrative examples, the thermal
diffusivity for solids is obtained by making a dual cubic spline interpolation of the amplitude of
the photothermal radiometry signal [23]. The derivatives of functions depending on the solution
of the Duffing equation [24–26] are calculated using the dual Runge–Kutta method. And the dual
Newton–Raphson algorithm is used to obtain the derivatives of the output angle in the RRRCR
spatial mechanism [27]. The Fortran code of the elemental dual functions as well as the dual
version for the mentioned algorithms are provided as additional material to this article and can be
downloaded from http://dual-algorithms.frp707.esy.es/. From these codes, the dualization
of many other functions and algorithms can be obtained.

2

Derivatives using dual numbers

How dual numbers can be used for calculating numerical derivatives can be seen in [28, 29]. Here
we briefly review the essential ideas, bearing in mind a numerical implementation.
Let f : R → R be an analytic function. Expanding this function in a Taylor Series and
evaluating at x̂ = x + ǫ (note that we have taken the dual part to be 1) we get
✯0
✟
✟

✟2 ).
O(ǫ
f (x + ǫ) = f (x) + f ′ (x)ǫ + ✟

(1)

As we can see, f (x) + f ′ (x)ǫ is a dual number, so we associate the dual function
fˆ(x̂) = f (x) + f ′ (x)ǫ

(2)

to it. As in the case of the complex numbers, there is an isomorphism between the dual numbers
and R2 , thus a dual number can be written as
r̂ = {a, b}.

(3)

From now on we will use the notation given by Eq. (3), thus Eq. (2) will be written as

where f0 = f (x) and f1 = f ′ (x).

fˆ(x̂) = {f0 , f1 },

2

(4)

The next step is to dualize the composition of f (x) with another function g(x). From the chain
rule, we get
fˆ(ĝ) = {f0 (g0 ), f1 (g0 )g1 }.
(5)
For instance, if h(x) = f (g(u(x))), the dual component of ĥ(x̂) = fˆ(ĝ(û(x̂))) will be h′ (x). That
is the power of the dual method of obtaining derivatives. We only need to implement the chain
rule once. This makes the dual number method of obtaining derivatives an AD method.
The generalization to second derivatives is straightforward. To this end we define a new dual
number
r̃ = a + b ǫ1 + c ǫ2 ,
r̃ = {a, b, c},

(6)

with a, b and c being real numbers and ǫ1 and ǫ2 having the following multiplication table
1
ǫ1
ǫ2

1

ǫ1

ǫ2

1
ǫ1
ǫ2

ǫ1
ǫ2
0

ǫ2
0
0

.

(7)

Now, evaluating the Taylor expansion for f (x) in x̃ = x + ǫ1 we have,
1 ′′
f (x)ǫ21 ,
2!
1
f (x̃) = f (x) + f ′ (x)ǫ1 + f ′′ (x)ǫ2 .
2!
These expressions are of the form of Eq. (6) so we can write
f (x̃) = f (x) + f ′ (x)ǫ1 +

(8)
(9)

1
f˜(x̃) = {f (x), f ′ (x), f ′′ (x)}.
(10)
2!
The fact that the third component is a half of the second derivative does not matter. We can
choose
f˜(x̃) = {f (x), f ′ (x), f ′′ (x)}
(11)

as our definition for this new dual function (in fact, by choosing ǫ21 = 2ǫ2 the factor 1/2 disappears).
Also the inappropriate name, new dual function, can be changed to dual function; leaving to the
context of the problem (number of components of f ) if we are talking about a simple dual function
or of an extended dual function of 3 components. Clearly Eq. (10) can be generalized to obtain
higher order derivatives. For instance, to obtain derivatives until order n, define ǫ0 = 1, ǫ1 , ǫ2 , . . . ,
ǫn having the multiplication table

0
if i + j > n,
ǫi · ǫj =
(12)

ǫi+j otherwise,
with i, j = 0, 1, . . . , n.
It is also interesting to prove Eqs. (8, 9, 10) from a matrix algebra
we define the following matrices having the multiplication table (7),





0
0 1 0
1 0 0





1 = 0 1 0 , ǫ1 = 0 0 1 , ǫ2 = 0
0
0 0 0
0 0 1
3

point of view. To this end

0 1
0 0 .
0 0

(13)

Now, since

X = x 1 + ǫ1

is in its Jordan canonical form, we will have [30]


f (x) f ′ (x) f ′′ (x)/2!
f (X) =  0
f (x)
f ′ (x)  ,
0
0
f (x)
1
f (X) = f (x) 1 + f ′ (x) ǫ1 + f ′′ (x) ǫ2 .
2!

(14)

(15)
(16)

Thus, the derivatives of f are determined by calculating the matrix function f (X). This matrix
function can be calculated by using the Cauchy integral formula for operator-value functions (see
for example Sec. 16.8 of [31]).
A generalization of Eqs. (15, 16) to obtain higher order derivatives, can be done by calculating
f (X) with
(17)
X = x 1(n+1)×(n+1) + ǫ1;(n+1)×(n+1) ,


ǫα;(n+1)×(n+1) i,j = δi,j−α ; being δi,j the Kronecker’s delta function, and α = 1, 2, . . . , n.
From Eq. (10), the analog of x̂ = {x, 1} will be x̃ = {x, 1, 0} and the analog of Eq. (4) will be
f˜(x̃) = {f0 , f1 , f2 },

(18)

where f0 = f (x), f1 = f ′ (x) and f2 = f ′′ (x). Similarly, we will have
f˜(g̃) = {f0 (g0 ), f1 (g0 )g1 , f2 (g0 )g12 + f1 (g0 )g2 },

(19)

for the composition of two dual functions.
Eq. (19), is of central importance in dualizing a function or algorithm. Notice that, unlike the
finite difference methods, the use of Eq. (19) to obtain derivatives does not have the problems of
truncation or cancellation errors.

2.1

Fortran implementation

In order to facilitate the use of the dual number giving by Eq. (6), it is convenient to introduce a
derived data type. In the Fortran programming language, this can be done as:
type, public :: dual2
real(8) :: f0, f1, f2
end type dual2
Now, the dual number of Eq. (6) can be written in Fortran as dual2(a,b,c). This derived data
type is defined in the module called dual2 mod. Also many Fortran functions and operators are
overloaded to deal with this kind of number. Such a module is included in the additional material
to this article.

4

3

Dualization of algorithms

3.1

Dual Newton–Raphson algorithm

Let us suppose that we are given equation
F (u(x), x) = 0

(20)

and that u′ (x0 ) is required. In the case when a closed form expression for u(x) can be obtained,
the derivatives can be calculated by writing u(x) in its dual form ũ(x̃). If there is no closed form
expression for u(x), applying the chain rule yields
u′ (x0 ) = −

1
∂F
(u(x0 ), x0 ).
∂F/∂u(u(x0 ), x0 ) ∂x

(21)

Assuming that there are closed form expressions for ∂F/∂x and ∂F/∂u, the derivative u′ (x0 ) can
be calculated by solving numerically Eq. (20) and substituting u(x0 ) into Eq. (21). If ∂F/∂x
and ∂F/∂u cannot be obtained in closed form, then there is no conventional, easy, accurate, and
precise way to calculate u′ (x0 ). However if we code the dual version of some numerical solution
method to solve Eq. (20), then we will obtain u(x0 ) and automatically its derivatives. In the
present paper, we dualize the Newton–Raphson method [32, 33].
Let f (q) be a real function of a real variable. The Newton–Raphson method allows finding a
solution of the equation f (q) = 0. Starting to look for a solution in q0 , the algorithm determines
an approximate solution as1
q1 = q0 −

f (q0 )
f ′ (q0 )

..
.
qn+1 = qn −

f (qn )
.
f ′ (qn )

(22)

In order to obtain u′ (x0 ) from Eq. (20) we proceed as follows: Write F̃ for the dual version
of the function F . Define the dual starting point u˜0 = {u0 , 0, 0}. Since we need F ′ (u), we
define x˜0 = {x0 , 0, 0} and ũ = {u0 , 1, 0}, then from f˜1 = F̃ (ũ, x˜0 ) we obtain the first and second
derivatives of F with respect to u. This is an extra bonus of the method: we do not need to worry
about the derivative of F . In practice, this is an issue, and the derivative must be provided by
hand. Finally, write Eq. (22) in its dual form, keeping in mind that since we want u′ (x) we need
to evaluate F̃ in u˜0 and x̃ = {x0 , 1, 0}.
Putting u˜0 = u0d, x˜0 = x0d, and F̃ = fd; the essential part of the method in Fortran is
x0d = x0
u0d = u0
do k = 1,
fx =
fu =
u0d =
end do

n
fd(u0d, dual2(x0,1d0,0d0))
fd(dual2(u0d%f0,1d0,0d0), x0d)
u0d - fx/(fu%f1)

1

The convergence of the Newton–Raphson method is not a trivial subject. In some cases, a small variation of q0
causes a convergence/divergence of the method. For a study of the convergence of the method, we refer the reader
to [34].

5

This algorithm will produce ũ(x). Using Eq. (19), we can construct the general dual function for
u: ũ(g̃). The complete algorithm is coded in the module NR dual mod of the additional material.
We also coded Halley’s method [35], so the user can choose between the simple Newton–Raphson
method or Halley’s method. Some examples are presented in Section 4.

3.2

Dual cubic spline interpolation

Let P : A ⊂ R → B ⊂ R be the function representing a cubic spline interpolation and let f be a
real function such that f (P (x), x) is defined. The problem is to find f ′ (P (x), x) with x ∈ A. We
solve this problem by writing the dual version of the natural cubic spline interpolation.
Natural cubic splines
Suppose that we have n data points {(x1 , y1 ), . . . , (xn , yn )}, n > 1. A cubic spline is a spline
constructed of piecewise third-order polynomials
Yi (t) = ai + bi t + ci t2 + di t3 ; t ∈ [0, 1], i = 1, . . . , n − 1,

(23)

which pass through this set of points.
In determining the coefficient of the ith piece of the spline, a linear system of 4(n − 1) − 2 equations and 4(n − 1) unknowns will appear [36]. Two more equation can be obtained by demanding
that Y ′′ (0) = 0 and Y ′′ (1) = 0. This yields the so called natural cubic spline interpolation. With
this, the coefficients of Yi are given by
ai = y i
bi = Di
ci = 3(yi+1 − yi ) − 2Di − Di+1

di = 2(yi − yi+1 ) + Di + Di+1

(24)

where the Ds numbers are determined by solving the symmetric tridiagonal system:
TD = R
with



2 1
 1 4 1


1 4


..
T=
.





1
..
.
1








D=





D1
D2
D3
..
.
Dn−2
Dn−1
Dn
6

(25)







..

.


4 1

1 4 1 
1 2












(26)

(27)



y2 − y1
y3 − y1
y4 − y2
..
.






R = 3

 yn−1 − yn−3

 yn − yn−2
yn − yn−1








.





(28)

The inversion of an n × n tridiagonal matrix can be done by an O(n) algorithm [32], and in general
it is this kind of algorithm which is used in order to find the inverse of the matrix T in Eq. (26).
However, as we will show, an analytical formula for the inverse of T can be obtained.
Let us consider the n × n nonsingular tridiagonal matrix


a 1 b1
 c1 a2 b2





.
.
.
.
.
.
.
M=
c2




.
.
.
.

.
. bn−1 
cn−1 an
The inverse of such a matrix can be written as [37]:

(−1)i+j bi . . . bj−1 θi−1 φj+1 /θn
−1
Mij =

(−1)i+j cj . . . ci−1 θj−1 φi+1 /θn

if i 6 j,
(29)
if i > j,

where θi is obtained by solving the recurrence equation

θi = ai θi−1 − bi−1 ci−1 θi−2 , for i = 2, . . . , n,

(30)

with θ0 = 1 and θ1 = a1 . Then φi is obtained by solving the recurrence equation
φi = ai φi+1 − bi ci φi+2 , for i = n − 1, . . . , 1,

(31)

with φn+1 = 1 and φn = an .
The use of Eq. (29) is not an efficient way to find the inverse of a tridiagonal matrix. However
it can be used to deduce an analytical formula for the inverse of the matrix T.
Applying Eqs. (29, 30 and 31) to Eq. (26), we get
 √ s
√ s
2− 3) +(2+ 3)

(
if s 6= n,
2
(32)
θs =
√ n
√
√ n
√

 (2+ 3) (−3+2 3) − (2− 3) (3+2 3)
if s = n,
2
 √ 1+n
√ s
√ s
√ 1+n
2− 3)
(2+ 3) + (2− 3) (2+ 3)

(
if s 6= 1,
2
φs =
(33)
√
√
√
√
n
n

 (2+ 3) (−3+2 3)−(2− 3) (3+2 3)
if s = 1.
2
Defining inv as

inv(s, k) = (−1)s+k θs−1 φk+1 /θn ,
7

(34)

writing explicitly φk+1 /θn and carrying out some elementary algebra, we obtain
i
(−1)s+k 1 + (α† /α)s−1 h (k+1) ln α† + (s−1) ln α
† (s+k) ln α + n ln(α† /α)
inv(s, k) =
α
e
+
α
e
,
2
β † − (α† /α)n β

with

√
3
√
α† = 2 − 3
√
β =2 3+3
√
β † = 2 3 − 3.

(35)

α=2 +

From this, the in verse of the matrix T is
(
inv(s, k) if s 6 k,
−1
Tsk =
inv(k, s) if s > k.

(36)

(37)

Thus the coefficients Ds are given by
Ds =

n
X

T−1
sk Rk , s = 1, . . . , n.

(38)

k=1

Now that all the coefficients have been determined, the parametric equation for the interpolated
points will be
ri (t) = {xi + (xi+1 − xi ) t, Yi (t)} .
(39)

Eliminating the parameter t, the equation for the ith polynomial of the normal cubic spline will
be


x − xi
, x ∈ [xi , xi+1 ].
(40)
Pi (x) = Yi
xi+1 − xi
The interpolated points in the whole interval [x1 , xn ] are
P (x) =

n−1
[

Pi (x).

(41)

i=1

Once P (x) has been promoted to a dual function (see the module cubic spline dual mod of
the aditional material), the derivatives of an arbitrary function f = f (P (x), x) for any x ∈ [x1 , xn ]2 ,
as well as P ′ (f (x), x), are calculated by writing f in its dual form. Some examples are presented
in Section 4.

3.3

Dual Runge–Kutta algorithm

Let us consider the following ordinary differential equation (ODE)
f ′′ (t) = F (t, f, f ′ ),

(42)

with initial conditions f (t0 ) = f0 and f ′ (t0 ) = v0 . Let g(t) be a function defined in such a way
that f (g(t)) and g(f (t)) are defined. The problem we want to address is to find the derivatives of
the composition of the functions f and g.
2

Thus, dualizing an interpolation method we can calculate the derivatives of experimental data. Nevertheless we
recommend using the methodology presented in [13, 22] instead, unless the derivatives between the nodes [xi , xi+1 ]
are required.

8

Dual Runge–Kutta 4th order method
One of the most often used methods for numerically solving ODEs is the Runge–Kutta Method
[32, 33, 38]. In order to calculate f ′ (g(t)), f ′′ (g(t)), g′ (f (t)) and g′′ (f (t)), we will dualize the 4th
order Runge–Kutta method (RK4). Putting x2 = f ′ , x1 = f , u1 (t, x1 , x2 ) = x2 , u2 (t, x1 , x2 ) =
F (t, x1 , x2 ), the RK4 method produces x2 and x1 and hence f ′′ (t). From this, the dual version of
f as a function of the real variable t (see the rk4 function in the module runge kutta dual mod
of the additional material) is

f˜(t) = f (t), f ′ (t), f ′′ (t) .
(43)
The dual version of f for any dual variable ũ, namely f˜(ũ), can be constructed following Eq.
(19). (See the rk4dual function in module runge kutta dual mod of the additional material).
Once we have f˜(ũ), the derivatives are calculated from the compositions f˜(g̃(t̃)) and g̃(f˜(t̃))—
actually, we can calculate the derivatives of any function G(f (t), t). Some examples are presented
in Section 4.

4

Worked examples

All the used functions are coded in the additional material to this article.

4.1

Dual Newton–Raphson example 1

Suppose that u(x) is defined by the equation f (u, x) = 0 where
f (u, x) = cos(u x) − u3 + x + sin(u2 x).
Let us consider the functions g1 (x) = sin(u(x)) + x and g2 (x) = u(sin x + x2 ). The problem is to
find u(x0 ), u′ (x0 ), u′′ (x0 ), g1 (x0 ), g1′ (x0 ), g1′′ (x0 ), g2 (x0 ), g2′ (x0 ), g2′′ (x0 ), with x0 = 0.7. This can
be accomplished by writing the Newton–Raphson algorithm in its dual form. Such an algorithm
is coded in the function NR_dual(NRkind,u0,n,fd,gdual). In this function, NRkind is the kind
of method used. NRkind = NR1 is for the simple Newton–Raphson method and NRkind = NR2 is
for Halley’s method, u0 is the dual point where the method starts to look for a solution, n is the
number of iterations, fd is the equation to solve; in this example it will be f˜(ũ, x̃), the dual version
of f (u, x); gdual is the dual point where we want to evaluate the functions. Writing the module
equation_mod.f90, where the function f (u, x) is coded in its dual version, we have
1
2
3
4
5
6
7
8
9
10
11
12
13

module equation_mod
use dual2_mod
implicit none
contains
function f(u,x) result(f_result)
type(dual2), intent(in) :: u, x
type(dual2)
:: f_result
f_result = cos(u*x) - u**3 + x + sin(u**2*x)
end function f
end module equation_mod
9

—a driver program could be as:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30

include "dual2_mod.f90"
include "equation_mod.f90"
include "NR_dual_mod.f90"
program nr1
use dual2_mod
use NR_dual_mod
use equation_mod
implicit none
real(8), parameter :: x = 0.7d0, u0 = 1.6d0
type(dual2)
:: xd
integer
:: n
xd = dual2(x,1.d0,0.d0)
n = 10
print*,"The function u and its derivatives in x = 0.7 are"
print*,NR_dual("NR1", u0, n, f, xd)
print*," "
print*,"The function g1 and its derivatives in x = 0.7 are"
print*, sin(NR_dual("NR1", u0, n, f, xd)) + xd
print*," "
print*,"The function g2 and its derivatives in x = 0.7 are"
print*,NR_dual("NR1", u0, n, f, sin(xd) + xd*xd)
end program nr1

Notice the advantages of this method: we can calculate very complicated derivatives involving
u(x) by using the dual functions. The results are shown in Table 1.
Table 1: Results for the Newton–Raphson example 1.

4.2

u(x0 )

u′ (x

1.3085

0.1163

0)

′′ u(x

0)

−0.9337

g1 (x0 )

g1′ (x0 )

g1′′ (x0 )

g2 (x0 )

g2′ (x0 )

g2′′ (x0 )

1.6658

1.0301

−0.2551

1.2963

−0.2556

−1.1425

Dual Newton–Raphson example 2

The example here presented concerns the RRRCR spatial mechanism [27]. Following Eq. (2)
of [27] and the definitions given in the aforementioned reference, the output angle φ as a function

10

of the input angle θ is given by the equation
a2 c21 c22 − 2ac1 c22 s1 (b − e) − 2ac21 c22 L cos θ + 2ac1 c22 R cos φ − c21 c22 (b − e)2 +

2c1 c22 Ls1 (b − e) cos θ + 2c1 c2 Ls2 (b − e) sin θ − 2c1 Rs2 (b − e) sin φ−

2Rs1 (b − e) cos φ + (b − e)2 − c21 c22 l2 + c21 c22 L2 + c21 c22 R2 cos2 φ−

2c21 c2 LR sin θ sin φ − c21 R2 1 − 2 sin2 φ − 2c1 c22 LR cos θ cos φ−
2c1 c2 LRs1 s2 sin θ cos φ + 2c1 R2 s1 s2 sin φ cos φ + R2 cos2 φ = 0.

(44)

Table 2: Parameters used for the Newton–Raphson example 2.
Parameter

Value

L
l
a
R
s1
s2
b
e

0.3933578023
0.4174323687
0.9526245468
0.4484604992
0.6298138891
0.2506389576
2.0
1.0

Now, if one is interested in calculating the velocities and accelerations, complicated functions
involving φ and φ′ (θ) will appear. As an example, let us consider the function f (x) = 2 sin2 x and
the point x0 = 2.0. The problem is to calculate f (φ(x0 )), f ′ (φ(x0 )), f ′′ (φ(x0 )), φ(f (x0 )), φ′ (f (x0 ))
and φ′′ (f (x0 )) for the same set of parameters given in [27], and reproduced in Table 2 for clarity.
Analogously to the example in Section 4.1, we obtain the results shown in Table 3. Notice that
the solution of Eq. (44) for θ = 2 is not unique. The values in Table 3 are for φ = 20.9846.
Table 3: Results for the Newton–Raphson example 2.

4.3

f (φ(x0 ))

f ′ (φ(x0 ))

f ′′ (φ(x0 ))

φ(f (x0 ))

φ′ (f (x0 ))

φ′′ (f (x0 ))

1.4279

-1.7693

-1.2856

1.7817

-1.6171

-3.5137

Dual cubic spline interpolation example 1

Let y(x) be the normal cubic spline interpolation for the data shown in Table 4. The values y(x0 ),
y ′ (x0 ), f (x0 ), f ′ (x0 ), g(x0 ) and g′ (x0 ) with f (x) = x sin2 (y(x)), g(x) = y(x sin2 x) and x0 = 1.75
can be calculated by dualizing the normal cubic spline interpolation method. This has been done
in the function NCSplinedual(A,xd). In such a function, A is the matrix containing the points for
which the spline will be constructed (in this case, they would be the points of Table 4) and xd is
the dual point where we want to evaluate, in this case xd=[1.75_dp,1_dp,0._dp]. For instance
the first (second) component of NCSplinedual(A,xd) will be y(x0 ) (y ′ (x0 )).
The derivatives f ′ (x0 ) and g ′ (x0 ), respectively, are obtained by taking component f1 of
xd*sin(NCSplinedual(A,xd))**2
and
11

Table 4: Points used for the dual cubic spline interpolation example 1.
x

y

1.
1.25
1.5
1.75
2.
2.25
2.5
2.75
3.

0
0.22314355
0.40546511
0.55961579
0.69314718
0.81093022
0.91629073
1.0116009
1.0986123

NCSplinedual(A, xd*sin(xd)**2)
The results are given in Table 5.
Table 5: Results for the dual cubic spline interpolation example 1.

4.4

y(x0 )

y ′ (x0 )

f (x0 )

f ′ (x0 )

g(x0 )

g′ (x0 )

0.5596

0.5727

0.4931

1.1836

0.5272

0.2097

Dual cubic spline interpolation example 2

This example concerns the determination of the thermal properties of a solid using photothermal
radiometry [23]. In particular, we are interested in determining the thermal diffusivity from the
experimental data of Fig. 4 from [23]. The experimental data can be extracted using EasyNData
[39]. Assuming that such data are stored in the matrix pts, the dual cubic spline interpolation is
done by NCSplinedual(pts,x). Now, according to [23], the thermal diffusivity can be calculated
from
64 L2s fq
αs =
,
(45)
9π
where Ls = 522 µm is the thickness of the studied sample and fq is the frequency for which the
derivative of the amplitude of the radiometry signal is zero. The frequency fq can be calculated
using the Newton–Raphson algorithm. To this end we define the dual function fej2(x) containing
the first and second derivatives of the spline. The essential part of the code in Fortran is
function fej2(x) result(f_result)
...
auxresult = NCSplinedual(pts,x)
f_result = dual2(auxresult%f1,auxresult%f2,0.d0)
...
Taking x0 = 10.d0 as the initial point where the Newton–Raphson algorithm will start to look
for a solution, fq can be found by
...
12

fq = x0
do k=1,50
fauxd = fej2(dual2(fq,1.d0,0.d0))
fq
= fq - fauxd%f0/fauxd%f1
end do
...
After this, the obtained thermal diffusivity was αs = 6.00×10−6 m 2 s−1 which is in good agreement
with the value reported in [23, 40].

4.5

Dual Runge–Kutta example

Consider the Duffing equation studied in [26]
f ′′ (t) + 0.4 f ′ (t) + 1.1 f (t) + f 3 (t) = 2.1 cos(1.8 t)

(46)

f (0) = 0.3, f ′ (0) = −2.3.

(47)

with initial conditions
The values f (t), (f ◦ g)(t), (g ◦ f )(t), as well as their first and second derivatives at t = 1.0 (or
at some other t where the functions are defined) for the function g(t) = sin t (or some other function where the above compositions are defined) can be calculated by dualizing the Runge–Kutta
algorithm. Such an algorithm is dualized in the function rk4dual(u1,u2,t0,x10,x20,np,td). In
this function, x10 and x20 are the initial conditions for f (t0 ) and f ′ (t0 ) respectively (see Section
3.3 for details and also for the definitions of u1 and u2 ); td is the dual point where we want to
evaluate the solution; and np is the number of steps between t0 and and the real component of td.
In order to apply the function rk4dual(u1,u2,t0,x10,x20,np,td), it is convenient to write
a module containing the functions u1 (t, x1 , x2 ) and u2 (t, x1 , x2 ). This can be done as follows:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

module equation_S45_mod
implicit none
contains
function u1(t,x1,x2) result(f_result)
real(8), intent(in) :: t, x1, x2
real(8)
:: f_result
f_result = x2 + 0d0*t + 0d0*x1
end function u1
function u2(t,x1,x2) result(f_result)
real(8), intent(in) :: t, x1, x2
real(8)
:: f_result
f_result = 2.1d0*cos(1.8d0*t) - 0.4d0*x2 - 1.1d0*x1 - x1**3

13

end function u2
end module equation_S45_mod

19
20

From this, f (t), (f ◦ g)(t), (g ◦ f )(t), as well as their first and second derivatives at t = 1.0 can be
calculated with the following driver program:
include "dual_modules/dual2_mod.f90"
include "dual_modules/runge_kutta_dual_mod.f90"
include "eqs_modules/equation_S45_mod.f90"

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

program
use
use
use

test_S45_dual
dual2_mod
equation_S45_mod
runge_kutta_dual_mod

implicit none
integer, parameter
real(8), parameter
type(dual2)
type(dual2)

::
::
::
::

np = 100
x10 = 0.3d0, x20 = -2.3d0, t0 = 0.d0
solution0, solution1, solution2
tvard

tvard = dual2(1.d0,1.d0, 0.d0)
solution0 = rk4dual(u1,u2,t0,x10,x20,np,tvard)
solution1 = sin(rk4dual(u1,u2,t0,x10,x20,np,tvard))
solution2 = rk4dual(u1,u2,t0,x10,x20,np,sin(tvard))
print*,solution0
print*,solution1
print*,solution2
end program test_S45_dual

Assuming that this program is saved in the file test_S45_dual.f90, we can compile it by executing
the Fotran compiler in a terminal window3 . The results are shown in Table 6.
Table 6: Results for the Runge–Kutta example. The functions are evaluated at t = 1.0.

5

f

f′

f ′′

(f ◦ g)

(f ◦ g)′

(f ◦ g)′′

(g ◦ f )

(g ◦ f )′

(g ◦ f )′′

−0.7474

−0.1282

0.8140

−0.6797

−0.0940

0.6081

−0.7144

−0.1638

0.6608

Conclusions

After a dualization of the normal cubic spline interpolation algorithm, the Runge–Kutta algorithm, and the Newton–Raphson algorithm, it is possible to calculate the derivatives of functions
3

For instance, if the intel R Fortran compiler is used, ifort test S45 dual.f90 will produce the executable file.

14

efficiently, precisely, and accurately. So we can calculate the derivatives of functions depending
on the solution of algebraic or differential equations as well as functions resulting from the spline
interpolation of experimental data. As an added value to the normal cubic spline interpolation
algorithm, it is shown that a closed form expression for its coefficients can be obtained. Interesting
applications in science and engineering were studied. Those examples can be used as a guide to
dualize many other functions and algorithms. For example, it would be interesting to dualize the
trapezium rule although this would be only for academic purposes since there is not much to gain
because its components would be the integral, the first derivative (which is actually the function
to integrate), and the second derivative (which is actually the first derivative of the function to
integrate). Nevertheless, the dualization of an integration method could be necessary for dualizing
some other algorithm.

References
[1] W. Clifford, Preliminary sketch of biquaternions, Proc. London Mathematical Society 1 (1-4)
(1873) 381–395.
[2] E. Pennestrı̀, R. Stefanelli, Linear algebra and numerical algorithms using dual numbers,
Multibody System Dynamics 18 (2007) 323–344.
[3] H. Leuck, H.-H. Nagel, Automatic differentiation facilitates of-integration into steering-anglebased road vehicle tracking, IEEE Computer Society Conference on Computer Vision and
Pattern Recognition 2 (5) (1999) 2360.
[4] D. Piponi, Automatic differentiation, C++ templates, and photogrammetry, Journal of
Graphics, GPU, and Game Tools 9 (4) (2004) 41–55.
[5] J. A. Fike, J. J. Alonso (Eds.), The Development of Hyper-Dual Numbers for Exact SecondDerivative Calculations, Proceedings of the 49th AIAA Aerospace Sciences Meeting, Orlando
FL, USA, 2011.
[6] W. Yu, M. Blair, DNAD, a simple tool for automatic differentiation of Fortran codes using
dual numbers, Computer Physics Communications 184 (2013) 1446–1452.
[7] J. N. Lyness, C. B. Moler, Numerical differentiation of analytic functions, SIAM Journal on
Numerical Analysis 4 (1967) 202–210.
[8] R. E. Rowlands, T. Liber, I. M. Daniel, P. G. Rose, Higher-order numerical differentiation
of experimental information, Experimental Mechanics 13 (1973) 105–112.
[9] A. Griewank, On automatic differentiation, in Mathematical Programming: Recent Developments and Applications, M. Iri and K. Tanabe, eds., Kluwer, Dordrecht, The Netherlands,
1998.
[10] I. Knowles, R. Wallace, New finite difference formulas for numerical differentiation, Numerische Mathematik 70 (1995) 91–110.
[11] I. Khan, R. Ohba, Closed-form expressions for the finite difference approximations of first
and higher derivatives based on Taylor series, Journal of Computational and Applied Mathematics 107 (1999) 179–193.
15

[12] I. Khan, R. Ohba, New finite difference formulas for numerical differentiation, Journal of
Computational and Applied Mathematics 126 (2001) 269–276.
[13] J. Li, General explicit difference formulas for numerical differentiation, Journal of Computational and Applied Mathematics 183 (2005) 29–52.
[14] Z. Qian, C.-L. Fu, X.-T. Xiong, T. Wei, Fourier truncation method for high order numerical
derivatives, Applied Mathematics and Computation 181 (2006) 940–948.
[15] K. Ahnert, M. Abel, Numerical differentiation of experimental data: Local versus global
methods, Computer Physics Communications 177 (2007) 764–774.
[16] A. K. Singh, B. S. Bhadauria, Finite difference formulae for unequal sub-intervals using
Lagranges interpolation formula, International Journal of Mathematical Analysis 3 (2009)
815–827.
[17] Z. Wang, R. Wen, Numerical differentiation for high orders by an integration method, Journal
of Computational and Applied Mathematics 234 (2010) 941–948.
[18] F.-F. Dou, C.-L. Fu, Y.-J. Ma, A wavelet-Galerkin method for high order numerical differentiation, Applied Mathematics and Computation 215 (2010) 3702–3712.
[19] R. D. Neidinger, Introduction to automatic differentiation and MATLAB object-oriented
programming, SIAM Review 52 (3) (2010) 545–563.
[20] M. Ramachandran, Fast derivative computation using smooth X-splines, Applied Numerical
Mathematics 62 (2012) 1654–1662.
[21] V. I. Dmitriev, Z. G. Ingtem, Numerical differentiation using spline functions, Computational
Mathematics and Modeling 23 (2012) 179–193.
[22] H. Hassan, A. Mohamad, G. Atteia, An algorithm for the finite difference approximation
of derivatives with arbitrary degree and order of accuracy, Journal of Computational and
Applied Mathematics 236 (2012) 2622–2631.
[23] M. Depriester, P. Hus, S. Delenclos, A. H. Sahraoui, New methodology for thermal parameter
measurements in solids using photothermal radiometry, Review of Scientific Instruments 76
(2005) 074902–1–074902–6.
[24] S. Wiggins, Introduction to Applied Nonlinear Dynamical Systems and Chaos, SpringerVerlag, New York, 1990.
[25] H. Nijmeijer, H. Berghuis, On Lyapunov control of the Duffing equation, IEEE Transactions
on Circuits and Systems–I (42) (1995) 473–477.
[26] Elçin Yusufoǧlu, Numerical solution of Duffing equation by the Laplace decomposition algorithm, Applied Mathematics and Computation (177) (2006) 572–580.
[27] J. J. Cervantes-Sánchez, J. M. Rico-Martı́nez, V. H. Pérez-Muñoz, A. Bitangilagy, Function
generation with the RRRCR spatial linkage, Mechanism and Machine Theory 74 (2014)
58–81.
16

[28] Y.-L. Gu, J. Y. S. Luh, Dual-number transformation and its applications to robotics, IEEE
Journal of Robotics and Automation 3 (6) (1987) 615–623.
[29] H. H. Cheng, Programming with dual numbers and its applications in mechanisms design,
Engineering with Computers 10 (4) (1994) 212–229.
[30] N. J. Higham, Functions of Matrices: Theory and Computation, Society for Industrial and
Applied Mathematics, Philadelphia, PA, USA, 2008.
[31] S. Hassani, Mathematical Physics. A modern Introduction to its Foundations, Springer, New
York, 2000.
[32] W. H. Press, S. A. Teukolsky, W. T. Vetterling, B. P. Flannery, Numerical Recipes in Fortran
77 The Art of Scientific Computing, 2nd Edition, Vol. 1, Cambridge Univerity Press, 1995.
[33] H. Cohen, Numerical Approximation Methods, Springer-Verlag, New York, 2010.
[34] H. Susanto, N. Karjanto, Newtons methods basins of attraction revisited, Applied Mathematics and Computation 215 (2009) 1084–1090.
[35] W. Gander, On Halley’s iteration method, The American Mathematical Monthly 92 (1985)
131–134.
[36] R. H. Bartels, J. C. Beatty, B. A. Barsky, An Introduction to Splines for Use in Computer
Graphics and Geometric Modeling, Morgan Kaufmann, Los Altos, CA, 1987.
[37] R. A. Usmani, Inversion of a tridiagonal jacobi matrix, Linear Algebra and its Applications
212213 (0) (1994) 413 – 414.
[38] D. F. Griffiths, D. J. Higham, Numerical Method for Ordinary Differential Equations,
Springer, London, 2010.
[39] P. Uwer, Easyndata: A simple tool to extract numerical values from published plots, arXiv:
0710.2896v1 [physics.comp-ph].
[40] A. L. Edwards, A compilation of thermal property data for computer heat-conduction calculations, UCRL-50589, University of California Lawrence Radiation Laboratory, 1969.

17

