Numeric-Symbolic Composite Derivative Calculations

I. INTRODUCTION

Composite function calculations arise in all areas of scientific and engineering computing. Low-order applications present no technical problems for manually derived calculations. High-order applications, however, rapidly become cumbersome to manage the bookkeeping involved. Alternatively, symbolic tools are available, such as Maple, Mathematica, Octave (used herein) and others, for generating arbitrary order derivative calculations. Unfortunately, run-time penalties make it impractical to embed a algebraic symbolic program within an application program. Since 1857 the work of Faa Di Bruno [1],[2],[3] has served as the gold standard for computing composite derivatives. This work combines the power of symbolic computing with the simple elegance of George Scott's 1861 formulation [4]. For five requested composite derivative calculations, the symbolic software generates a $5\times 5$ lower triangular array, that is embedded in the application software for computing numerical derivative values. A numerical example is provided to demonstrate the effectiveness of the proposed algorithm.

1.1 Faa Di Bruno Formulation

Since 1857 scalar composite function calculations have been theoretically handled by invoking the celebrated work of Faà Di Bruno [1],[2],[3]. His formula for computing the $n$th derivative of the function $u=F(G(x))$ is given by

where the summation process is subject to the following integer constraint conditions

No additional derivative calculations are required. Not surprisingly the calculation represents a daunting bookkeeping challenge.

1.2 Scott's Composite Derivative Formulation

In 1861 George Scott [4] proposed the following alternative composite function derivative formulation

which is analytically much simpler than Faà Di Bruno's identity. Unfortunately, the requirement for generating hand-derived derivatives for ${\frac{d^n}{d{x}^n}{G}^k(x)|}_{G=0}$ , where $n$ denotes the $n$th derivative and $k$ denotes derivative summation index, has limited its broader application in the computational research community.

II. NUMERIC-SYMBOLIC COMPOSITE DERIVATIVE FORMULATION

This paper presents a composite derivative algorithm that combines: (1) the simplicity of Scott's formulation [4], and, (2) the power of algebraic symbolic manipulation. Only two function routines are required (see Figure 1 and Appendix A). Symbolic function derivatives are computed for a generic $G(x)$ function. The use of a generic function (see Figure 1), eliminates the need for embedding a symbolic program in the application software. The algorithm consists of two parts. First, a symbolic preprocessor (see Figure 1) computes the derivatives for ${\frac{d^n}{d{x}^n}{G}^k(x)|}_{G=0}$ and enforces a $G=0$ boundary condition on each derivative. The derivative calculations become very complicated as the requested number of composite derivatives increase. Nevertheless, after the boundary condition $G=0$ is enforced, the elements of the output derivative array are very simple. The symbolically generated derivative array consists of a lower triangular $n\times n$ array; which is copied into the application software for numerical derivative calculations. The elements of the array consist of components of the vector array $\{V_1,V_2,\dots ,V_n\}$ , where $V_i$ denotes $\frac{d^i}{d{x}^i}G$ . The array is only recomputed if the number of requested derivative increases.

Second, the application program embeds the computational function routine appearing in Figure 1, for numerically computing the composite derivatives. An elegantly simple double summation algorithm generates the numerical results for the composite derivatives.

For $nd=five$ requested composite derivatives, the symbolic preprocessor, presented in Appendix A, generates the following $dGdx$ lower triangular $nd\times nd$ array:

where ${V_i=\frac{d^{i}G}{d{x}^i}|}_{G=0}$ . The $dGdx$ array of Eq. 2 enables 1 through 5 composite derivatives to be evaluated. Elements of the array are nonlinear, for example, the element $dGdx(4,2)=8V_1{V}_3+6V_2^2$ . During the preprocessing operations, the following change of variables eliminates potential problems with processing the nonlinear terms:

The change of variables substitution starts with the highest order derivative and continues until the lowest order derivative variable has been replaced (see Appendix A).

An additional benefit of this approach is that for applications only requiring three composite derivatives, one only uses the sub portion of the array dGdx(1:3,1:3) for all numerical derivative calculations.

III. COMPOSITE DERIVATIVE CALCULATIONS

Taking the symbolically generated $dGdx$ array of Eq. 2 as input for the application software; composite derivatives are computed with following function routine: Function Composite $=$ Derivative(nd,x,DF,DG,Rfac) %% This function computes composite derivatives for nd derivatives. %% The Maximum number of derivatives is five(5), otherwise the %% symbolic preprocessor is recomputed to generate a new $dGdx$ array. %% INPUT: %% nd number of composite derivatives requested %% x evaluation point for the composite derivative

where $x$ denotes the evaluation point for the derivatives

%% DF, DG are user supplied nd x 1 arrays of derivatives numerically

%% evaluated at the application points

%% $Rfac=$ reciprocal factorial array nd x 1

%% OUTPUT:

%% Composite denotes the nd x 1 array of requested composite derivatives

%% %% Copyright(2025) James D. Turner, All Rights Reserved

%%

%% Introduce Current $\mathrm{G}(\mathrm{x})$ Derivative values in the following $dGdx$ array

%% dGdx is copied from the symbolic preprocessor of Appendix A

$dGdx(1,1)=V_1,dGdx(2,1)=V_2,dGdx(3,1)=V_3,dGdx(4,1)=V_4,dGdx(5,1)=V_5$

$dGdx(2,2)=2V_1^2$

$dGdx(3,2)=6V_1{V}_2$

$dGdx(4,2)=8V_1{V}_3+6V_2^2$

$dGdx(5,2)=10V_1{V}_4+20V_2{V}_3$

$dGdx(3,3)=6V_1^3$

$dGdx(4,3)=36V_1^2{V}_2$

$dGdx(5,3)=60V_1^2{V}_3+90V_1{V}_2^2$

$dGdx(4,4)=24V_1^4$

$dGdx(5,4)=240V_1^3{V}_2$

$dGdx(5,5)=120V_1^5$

%% Begin Main Derivative Summation Loop

for $n=1$ :nd %% Number of requested composite derivatives

[ $\mathrm{s}$ = 0; \quad % $\text{initialize each derivative summation}$ ]

for $k=1:n$ Number of summation terms for each derivative

s = s + DF(k)*dGdx(n,k)*Rfac(k);

endfor %% Sum current Derivative Series Expansion

Composite(n) = s; %% Save nth value of the composite derivative endfor %% double summation algorithm endfunction Figure 1: Numerical Composite Derivative Function Routine This double summation algorithm is extremely simple when compared with Faà di Bruno's formulation [1],[2],[3].

IV. NUMERICAL APPLICATION

This Section presents a numerical composite derivative example, where $f^{(5)}$ composite derivatives are requested. The assumed composite derivative function is given by:

There are two numerical evaluation points for the function appearing in Eq. (3). Namely, Evalpt = 0.5 for G(Evalpt), and $\mathrm{FEvalpt}=\cosh \left(\mathrm{Evalpt}\right)$ The user is assumed to have provided numerical values for the following analytic derivative arrays for $F$ and $G$ :

These numerical arrays are processed in the double summation algorithm presented in Figure 1. Assuming that the evaluation point for the composite derivative is $x=0.5$ , the above derivative arrays have the following truncated numerical values

DG = \left[ \begin{array}{lllll} 1.5e-01 & 1.01e+00 & 1.5e-01 & 1.01e+00 & 1.5e-01 \end{array} \right]

Using the double summation algorithm presented in Figure 1, the five requested composite derivatives are given by:

V. CONCLUSION

This work links (1) the power of symbolic processing for computing derivative values for $\frac{d^n}{d{x}^n}{G}^k(x){|}_{G=0}$ and the simplicity of Scott's algorithm presented in Eq. (1). Numerical composite derivative calculations are performed using the elegantly simple double summation algorithm presented in Figure 1. Up to five composite derivatives are handled by the algorithm. If more composite derivatives are required; the symbolic preprocessor of Appendix A is recomputed to generate an updated derivative array for $dGdx$ . The symbolic preprocessor performs derivative calculations and enforces a boundary condition on each derivative. The output of the symbolic preprocessor consists of a lower triangular $5\times 5$ matrix. The elements $dGdx$ functionally depend on the nonlinear derivatives of generic function $G(x)$ . The derivative array $dGdx$ is embedded in function routine presented in Figure 1 for computing the numerical composite derivatives. Assuming that the upper limit for the requested number of composite derivatives does not change; the preprocessor is only executed one time. The algorithm scales for tens to hundreds of composite derivatives. A numerical example is presented that demonstrates the effectiveness of the proposed methodology.

APPENDIX A

Symbolic Constrained Derivative Preprocessor

$dx={\int }_{G=0}^{x}dx$ The use of a generic function $G(x)$ eliminates the need for embedding a symbolic manipulation tool in the application software. The generation of symbolic derivatives is straightforward. Enforcing the $G=0$ boundary condition of the derivative expressions requires an indirect approach to avoid a potential loss of nonlinear terms. This is handled by introducing a change of variables to eliminate the explicit derivative expressions in the symbolic derivative results. The calculation is started by eliminating the highest derivatives first. The calculation generates a lower triangular array $dGdx$ , that is inserted in the application software as the function routine appearing in Figure 1. This array is computed only one time, and is only recomputed if the number of requested composite derivative increases.

function $dgk=$ Sym_dGdk_preprocessor ( nd )

%% Symbolic PreProcessor for composite derivatives

$\mathrm{\%}\mathrm{\%}$ Math Model: $\mathrm{phi}=\mathrm{f}(\mathrm{g}(\mathrm{x}))$ , $\mathrm{x}=$ evaluation point

%% INPUT: nd number of implicit derivatives required on output

$\mathrm{\%}\mathrm{\%}$ INPUT: nd Number of implicit derivatives required

%% OUTPUT: dgk nxn lower triangular array of derivatives @ G = 0

$\mathrm{\%}\mathrm{\%}$ Copyright (2025) James D. Turner All rights reserved

syms G(x); $\mathrm{\%}\mathrm{\%}$ Declare variables to be symbolic

syms k;

$\mathrm{\%}$ Create change of variables array to eliminate analytic derivs:

symbols = sym( zeros(1, nd)); % Preallocate symbolic array

for $\mathrm{i}=1:\mathrm{nd}$

symbols(i) $=$ sym(sprintf('v%d',i));

% Create symbols = {v1, v2, ..., vn}

end (G' = V1, G'' = V2, ..., G(n) = Vn)

old $=$ sym (zeros(nd,1));

new = sym( zeros(nd, 1)); %% allocate memory for derivative transformations

for $\mathrm{i}=1:\mathrm{nd}$

old(i,1) = diff(G(x), x, i); %% Load G derivatives

new(i,1) = symbols(i); %% Load {v1, v2, ..., vn}

endfor $\mathrm{\%}\mathrm{\%}$ Store derivative transformation equations

oldr = old(end:-1:1);

newr = new( end:-1:1 ); %% Reverse order of derivative & V arrays

$\mathrm{dgk}=\mathrm{sym}(\text{zeros(nd, nd)})$ ; $\mathrm{\%}\mathrm{\%}$ allocate array for enforcing G = o BC

%%= = = = = = = = = = = = = = = = = = =

%% Symbolic Preprocessor Loop

%%=--------

for $\mathrm{n}=1:\mathrm{nd}$

for $\mathrm{k}=1:\mathrm{n}$ %% summation terms for each derivative