1 Introduction

Analogous to a complex number z = a + i b where a and b are real numbers and i² = 1, a dual number is defined as $\hat{r}$ = a + ϵ b with a and b real numbers and ϵ² = 0. Such numbers were introduced by Clifford who also developed their algebra in the late nineteenth century ^[1].

The fact that ϵ² = 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, by 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,4,5,6]. Moreover, there are no papers addressing the dualization¹ of algorithms such as the Newton–Raphson algorithm, the Runge–Kutta algorithm, or the cubic spline interpolation method. The present 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 deriva- tives numerically ^{[5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22]}.

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 automatic differentiation methods (AD) ^[6,9,19,23] and in particular AD by using dual numbers are especially suited for calculating such 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 ^[24]. The derivatives of functions depending on the solution of the Duffing equation ^[25,26,27] 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 ^[28]. The Matlab code of the elemental dual functions as well as the dual version for the mentioned algorithms are provided as additional material to this article, which is available at ^[29].

2 Methodology

2.1 Derivatives using dual numbers

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.

2.2 Matlab implementation

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.

2.3 Dualization of algorithms

2.3.1 Dual Newton-Raphson algorithm

Let us suppose that we are given equation

(20)

and that u′(x₀) 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 $\tilde{u}\left(\tilde{x}\right)$.

If there is no closed form expression for u(x), applying the chain rule yields

(21)

Assuming that there are closed form expressions for ∂F/∂x and ∂F/∂u, the derivative u′(x₀) can be calculated by solving numerically Eq. (20) and substituting u(x₀) into Eq. (21). If ∂F/∂x and ∂F/∂u cannot be obtained in closed form, its numerical values can be obtained using AD. Nevertheless, by coding the dual version of some numerical solution method to solve Eq. (20), we will obtain u(x₀) and automatically its derivatives. In the present paper, we dualize the Newton– Raphson method ^{[34, 35]}.

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 q₀, the algorithm determines an approximate solution as²

(22)

In order to obtain u′(x₀) from Eq. (20) we proceed as follows: Write for the dual version of the function F . Define the dual starting point ${\tilde{u}}_0$ = {u₀, 0, 0}. Since we need F′(u), we define ${\tilde{x}}_0$ = {x₀, 0, 0} and $\tilde{u}$= {u₀, 1, 0}, then from $\tilde{f}$ ₁ = $\tilde{F}\left(\tilde{u},{\tilde{x}}_0\right)$ 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 $\tilde{F}$ in ${\tilde{u}}_0$ and $\tilde{x}$ = {x₀, 1, 0}.

Putting ${\tilde{u}}_0$ = u₀d, ${\tilde{x}}_0$ = x₀d, and $\tilde{F}$ = fd; the essential part of the method is

This algorithm will produce$\tilde{u}\left(x\right)$. Using Eq. (19), we can construct the general dual function for u:$\tilde{u}\left(g\right)$. The complete algorithm is coded in the function NR -dual of the additional material. We also coded Halley’s method ^[37], so the user can choose between the simple Newton–Raphson method or Halley’s method. Some examples are presented in Section 3.

2.3.2 Dual cubic spline interpolation

Let P: A ⊂ R → B ⊂ R be the function representing a cubic spline interpolation and let f b 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 (x₁, y₁), . . . , (x_n, y_n) , n > 1. A cubic spline is a spline constructed of piecewise third-order polynomials which oass thtough this set of points.

(23)

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 ^[38]. 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

(24)

where the D_s numbers are determined by solving the symmetric tridiagonal system:

(25)

with

(26)

(27)

(28)

The inversion of an n x n tridiagonal matrix can be done by an O(n) algorithm ^[34], 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 x n nonsingular tridiagonal matrix

The inverse of such a matrix can be written as ^[39]:

(29)

where θ_i is obtained by solving the recurrence equation

(30)

with θ₀ = 1 and θ₁ = a₁. Then φ_i is obtained by solving the recurrence equation

(31)

The use of Eq. (29) is not an eficient 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

(32)

(33)

(34)

writing explicitly φ_k+1/θ_n and carrying out some elementary algebra, we obtain

(35)

with

(36)

From this, the inverse of the matrix T is

(37)

Thus the coefficients D_s are given by

(38)

Now that all the coeficients have been determined, the parametric equation for the interpolated points will be

(39)

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

(40)

The interpolated points in the whole interval [x₁, x_n] are

(41)

Once P (x) has been promoted to a dual function (see the function NCSplinedual of the aditional material), the derivatives of an arbitrary function f = f (P (x), x) for any x ∈ [x₁, x_n],as well as P′(f (x), x), are calculated by writing f in its dual form. Some examples are presented in Section 3.

2.3.3 Dual Runge-Kutta algorithm

Let us consider the following ordinary diferential equation (ODE)

(42)

with initial conditions f (t₀) = f₀ and f′(t₀) = v₀. 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.

Dual Runge{Kutta 4th order method

One of the most often used methods for numerically solv- ing ODEs is the Runge–Kutta Method ^[34,35,40]. In or- der 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 x₂ = f′, x₁ = f , u₁(t, x₁, x₂) = x₂, u₂(t, x₁, x₂) = F (t, x₁, x₂), the RK4 method produces x₂ and x₁ and hence f′′(t). From this, the dual ver- sion of f as a function of the real variable t (see the rk4 function of the additional material) is

(43)

The dual version of f for any dual variable$\tilde{u}$, namely$\tilde{f}\left(\tilde{u}\right)$, can be constructed following Eq. (19). (See the rk4dual function of the additional material). Once we have$\tilde{f}\left(\tilde{u}\right)$, the derivatives are calculated from the com- positions $\tilde{f}\left(\tilde{g},\left(\tilde{t}\right)\right)$ and$\tilde{g}\left(\tilde{f},\left(\tilde{t}\right)\right)$—actually, we can calculate the derivatives of any function G(f (t), t). Some exam- ples are presented in Section 3.

3 Results

All the examples and required folders are included in the additional material to this article.

3.1 Dual Newton-Raphson example 1

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

3.2 Dual Newton-Raphson example 2

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.

3.3 Dual cubic spline interpolation example 1

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

3.4 Dual cubic spline interpolation example 2

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]}.

3.5 Dual Runge-Kutta example

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.

4 Conclusions

After a dualization of the normal cubic spline interpo- lation algorithm, the Runge–Kutta algorithm, and the Newton–Raphson algorithm, it is possible to calculate the derivatives of functions 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 aca- demic 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.

Referencias

[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-angle-based 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 Second-Derivative Calcu- lations, 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 Appli- cations, 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.

[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] A. Griewank, A. Walther, Evaluating Derivatives: Principles and Techniques of Algorithmic Differentiation, Society for Industrial Mathematics, 2008.

[24] 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.

[25] S. Wiggins, Introduction to Applied Nonlinear Dynamical Systems and Chaos, Springer-Verlag, New York, 1990.

[26] H. Nijmeijer, H. Berghuis, On Lyapunov control of the Duffing equation, IEEE Transactions on Circuits and Systems–I (42) (1995) 473–477.

[27] El¸cin Yusufoˇglu, Numerical solution of Duffing equation by the Laplace decomposition algorithm, Applied Mathematics and Computation (177) (2006) 572–580.

[28] J. J. Cervantes-S´anchez, J. M. Rico-Mart´ınez, V. H. P´erez-Mun˜oz, A. Bitangilagy, Function generation with the RRRCR spatial linkage, Mechanism and Machine Theory 74 (2014) 58–81.

[29] F. Pen˜un˜uri, O. Carvente, M. Zambrano-Arjona, R. Pe´on, C. A. Cruz-Villar, Matlab class for dual numbers, http://algorithmic-differentiation.frp707.esy.es/.

[30] 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.

[31] H. H. Cheng, Programming with dual numbers and its applications in mechanisms design, Engineering with Computers 10 (4) (1994) 212–229.

[32] N. J. Higham, Functions of Matrices: Theory and Computation, Society for Industrial and Applied Mathe- matics, Philadelphia, PA, USA, 2008.

[33] S. Hassani, Mathematical Physics. A modern Introduction to its Foundations, Springer, New York, 2000.

[34] 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.

[35] H. Cohen, Numerical Approximation Methods, Springer-Verlag, New York, 2010.

[36] H. Susanto, N. Karjanto, Newtons methods basins of attraction revisited, Applied Mathematics and Compu- tation 215 (2009) 1084–1090.

[37] W. Gander, On Halley’s iteration method, The American Mathematical Monthly 92 (1985) 131–134.

[38] 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.

[39] R. A. Usmani, Inversion of a tridiagonal jacobi matrix, Linear Algebra and its Applications 212213 (0) (1994) 413 – 414.

[40] D. F. Griffiths, D. J. Higham, Numerical Method for Ordinary Differential Equations, Springer, London, 2010.

[41] P. Uwer, Easyndata: A simple tool to extract numerical values from published plots, arXiv: 0710.2896v1 [physics.comp-ph].

[42] A. L. Edwards, A compilation of thermal property data for computer heat-conduction calculations, UCRL- 50589, University of California Lawrence Radiation Laboratory, 1969.

Notes

Author notes