How dual numbers can be used for calculating numer- ical derivatives can be seen in ^{[30, 31]}. 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 $\hat{x}$=x+ϵ (note that we have taken the dual part to be 1) we get

(1)

As we can see, f (x) + f′(x)ϵ is a dual number, so we associate the dual function

(2)

to it. As in the case of the complex numbers, there is an isomorphism between the dual numbers and R², thus a dual number can be written as

(3)

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

(4)

The next step is to dualize the composition of f (x) with another function g(x). From the chain rule, we get

(5)

For instance, if h(x) = f (g(u(x))), the dual component of $\hat{h}\left(\hat{x}\right)=\hat{f}\left(\hat{g}\left(\hat{u}\left(\hat{x}\right)\right)\right)$ 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 (forward mode of AD).

The generalization to second derivatives is straight- forward. To this end we define a new dual number

(6)

with a, b and c being real numbers and ϵ₁ and ϵ₂ having the following multiplication table

(7)

Now, evaluating the Taylor expansion for f (x) in $\tilde{x}$ = x + ϵ₁ we have,

(8)

(9)

These expressions are of the form of Eq. (6) so we can write

(10)

The fact that the third component is a half of the second derivative does not matter. We can choose:

(11)

as our definition for this new dual function (in fact, by choosing ϵ² = 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 ϵ₀ = 1, ϵ₁, ϵ₂, . . . , ϵn having the multiplication table

(12)

It is also interesting to prove Eqs. (8, 9, 10) from a matrix algebra point of view. To this end we define the following matrices having the multiplication table (7),

(13)

(14)

is in its Jordan canonical form, we will have ^[32]

(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 ^[33]).

A generalization of Eqs. (15, 16) to obtain higher order derivatives, can be done by calculating f(X) with

(17)

From Eq. (10), the analog of $\hat{x}$={x, 1} will be = {x, 1, 0} and the analog of Eq. (4) will be

(18)

where f₀ = f (x), f₁ = f′(x) and f₂ = f′′(x). Similarly, we will have

(19)

for the composition of two dual functions.

Eq. (19), is of central importance in dualizing a function or algorithm. Notice that, unlike the fnite diference methods, the use of Eq. (19) to obtain derivatives does not have the problems of truncation or cancellation errors.

In order to facilitate the use of the dual number given by Eq. (6), it is convenient to introduce a class for this kind of number. In the Matlab programming language, this can be done as:

Now, the dual number of Eq. (6) can be written as $\tilde{r}$ = dual(a,b,c). This class as well as many other overloaded functions are included in the additional material to this article.

Suppose that u(x) is defined by the equation f (u, x) = 0 where

Let us consider the functions g₁(x) = sin(u(x)) + x and g₂(x) = u(sin x + x2). The problem is to find u(x₀), u′(x₀), u′′(x₀), g₁(x₀), g₁′ (x₀), g₁′′(x₀), g₂(x₀), g₂′ (x₀), g₂′′(x₀), with x₀ = 0.7. This can be accom- plished 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, u₀ 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$\tilde{f}\left(\tilde{u},\tilde{x}\right)$, the dual version of f (u, x); gdual is the dual point where we want to evaluate the functions. The following script computes the required values.

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

The example here presented concerns the RRRCR spa- tial mechanism ^[28]. Following Eq. (2) of ^[28] and the definitions given in the aforementioned reference, the output angle $\oint$ as a function of the input angle $∲$ is given by the equation

(44)

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 sin₂ x and the point x₀ = 2.0. The problem is to calculate f (φ(x₀)), f′(φ(x₀)), f′′(φ(x₀)), φ(f (x₀)), φ′(f (x₀)) and φ′′(f (x₀)) for the same set of parameters given in ^[28], and reproduced in Table 2 for clarity.

Table 2
Parameters used for the Newton{Raphson example 2

Analogously to the example in Section 3.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 φ = 2.1351.

Table 3
Results for the Newton{Raphson example 2.

Let y(x) be the normal cubic spline interpolation for the data shown in Table 4.

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

The values y(x₀), y′(x₀), f (x₀), f′(x₀), g(x₀) and g′(x₀) with f (x) = x sin₂(y(x)), g(x) = y(x sin² x) and x₀ = 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=dual (1.75,1,0). For instance the first (sec- ond) component of NCSplinedual(A,xd) will be y(x₀) (y′(x₀)).

The derivatives f′(x₀) and g′(x₀), respectively, are obtained by taking component f1 of

and

The results are given in Table 5.

Table 5
Results for the dual cubic spline interpolation example 1

This example concerns the determination of the thermal properties of a solid using photothermal radiometry ^[24]. In particular, we are interested in determining the thermal di_usivity from the experimental data of table. 4 from ^[24]. The experimental data can be extracted using EasyNData ^[41]. Assuming that such data are stored in the matrix pts, the dual cubic spline interpolation is done by NCSplinedual (pts,x). Now, according to ^[24], the thermal diffusivity can be calculated from

(45)

where L_s = 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 is

Taking x₀ = 10.0 as the initial point where the Newton–Raphson algorithm will start to look for a solution, fq can be found by

After this, the obtained thermal diffusivity was αs =6.00 10⁻⁶ m 2s⁻¹ which is in good agreement with the value reported in ^{[24, 42]}.

Consider the Du_ng equation studied in ^[27]

(46)

with initial conditions

(47)

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 compo- sitions are defined) can be calculated by dualizing the Runge–Kutta algorithm. Such an algorithm is dualized in the function rk4dual(u₁,u₂,t₀,x₁₀,x₂₀,np,td). In this function, x₁₀ and x₂₀ are the initial conditions for f (t₀) and f′(t₀) respectively (see Section 2.3.3 for details and also for the definitions of u₁ and u₂); td is the dual point where we want to evaluate the solution; and np is the number of steps between t₀ and and the real component of td. The values for f (t), (f g)(t), (g f )(t), as well as their first and second derivatives at t = 1.0 are shown in Table 6.

Table 6
Results for the Runge-Kutta example. The functions are evaluated at t = 1:0.