The Complex Taylor Series Expansion Method

The Complex Taylor Series Expansion Method (CTSE)

The starting point to learning HYPAD (HYPercomplex Automatic Differentiation) is the Complex Taylor Series Expansion method (CTSE), also called “complex step”. This method is convenient because it uses only complex variables, complex algebra, and complex functions contained in any modern computer language. Therefore, one can immediately compute numerical derivatives using whatever language is preferred.

CTSE is an ingenious numerical method that can be used to compute highly accurate numerical derivatives. The essence of CTSE can be seen in the simple one dimensional Taylor series expansion of a complex-valued function

Eq. 1 $\qquad \qquad f(z)=\sum _{n=0}^{\infty }{\frac {f^{(n)}(a)}{n!}}(z-a)^{n}=f(a)+{\frac {f'(a)(z-a)}{1!}}+{\frac {f''(a)(z-a)^{2}}{2!}}+\dots$

where $f'$ denotes the first derivative of $f$ with respect to $x$, $f''$ the second derivative, and $f^{(n)}$ the $n$'th derivative. Evaluating the series with the complex number $z=x+ih$, and the real number $a=x\in \mathbb {R}$, then $z-a=ih$, where $h$ is the step size. As a result,

Eq. 2 $\qquad \qquad \qquad {\begin{aligned}f(x+ih)&=f(x)+f'(x)(ih)+{\frac {f''(x)}{2}}(ih)^{2}+{\frac {f'''(x)}{6}}(ih)^{3}+\cdots \\&={f(x)}+{ihf'(x)}-{\frac {h^{2}}{2}}f''(x)-{\frac {ih^{3}}{6}}f'''(x)+\cdots \end{aligned}}$

We see that the terms within the complex Taylor series alternate between real and imaginary. We can sort this equation into real terms and imaginary terms as

Eq. 3 $\qquad \qquad \qquad \qquad \qquad \quad \Re (f(x+ih))=f(x)-{\frac {h^{2}}{2}}f''(x)+\operatorname {H.O.T.}$

Eq. 4 $\qquad \qquad \qquad \qquad \qquad ~~\Im (f(x+ih))=hf'(x)-{\frac {ih^{3}}{6}}f'''(x)+\operatorname {H.O.T.}$

where $\operatorname {H.O.T.}$ denotes “Higher Order Terms”. Equation 3 shows that the imaginary terms divided by $h$ can be written as

Eq. 5 $\qquad \qquad \qquad \qquad \qquad \qquad \qquad f'(x)={\frac {\Im (f(x+ih))}{h}}+O(h^{2})$

As a result, the first imaginary term contains the first derivative ${df \over dx}$ times $x$, and the first derivative can be approximated using the operation

Eq. 6 $\qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad f'(x)\approx {\frac {\Im (f(x+ih))}{h}}$

Equation 6 is not exact; however, as shown below, the terms of $O(h^{2})$ can be made to be numerically insignificant.

The key to CTSE is that it does not involve subtractive error, that is, it does not involve the subtraction of terms. As a result, the step size $h$ can be made very small. A small step size is valuable as it reduces truncation error by driving the higher order terms, represented by $O(h^{2})$, below machine precision. The end result is that, although Eq. 3 is not exact, it can be used to compute highly accurate derivatives.

A small $h$ is advantageous but what constitutes a sufficiently small $h$? A numerical study shown below suggests that one should use a value of $h<10^{-8}$ or smaller. However the size of the parameter being perturbed should also be considered. In other words, for example, if the units of the parameter of interest is $10^{-10}$, CTSE will not provide accurate results using a value of $h=10^{-8}$. Hence as a practical solution, $h$ should be at most $10^{-10}$ times the parameter of interest, but it can be smaller. Since the smallest number that can be represented with double precision arithmetic is $10^{-308}$, it follows that the perturbation applied must be greater than this value to prevent underflow. More practically, the variable of interest and its perturbation must be greater than $10^{-308}$. To summarize, there is a huge swath of numbers that can be used for the step size $h$ to obtain machine precision accuracy. A practical solution $h$ should be at least $10^{-10}$ times the parameter of interest. The suggested limits for $h$ are

$\qquad \qquad \qquad \qquad \qquad \qquad \qquad \quad {\begin{array}{cc}10^{-308}<h<&10^{-10}*parameter\end{array}}$

An additional aspect to CTSE is its effects on the real portion of $f(x+ih)$. As shown in equation 4, the real results get “polluted” by terms of order $h^{2}$. However, as stated previously, these terms are made numerically insignificant through the use of a small $h$. As a result, we numerically obtain the same real-valued result for $f(x+ih)$ as for $f(x)$.

In summary, if we evaluate $f(x+ih)$ using a sufficiently small step size $h$ instead of only $f(x)$, we obtain real result $f(x)$ and $f'(x)$.

1 Analogy with the finite difference method

CTSE has a direct analogy with the well-known finite difference (FD) method for estimating a derivative using a numerical approximation. Figure 1 shows a graphical representation comparing FD and CTSE. In this figure, the first order derivative of a function $f$ with respect to parameter $x$ is determined numerically by both approaches, FD by perturbing a parameter of interest along the real axis ($\Re$) by a step size ${\hat {h}}$, and CTSE by perturbing along the imaginary axis ($\Im$) by a step size $h$. An estimate for the derivative using a forward finite difference approach is

Eq. 7 $\qquad \qquad \qquad \qquad \qquad \qquad \qquad \quad {\frac {df(x)}{dx}}\approx {\frac {f(x+{\hat {h}})-f(x)}{\hat {h}}}$

In an analogous manner, CTSE perturbs the parameter of interest by a step size $h$ but along the imaginary axis. That is, real-variable parameter of interest is made complex-variable and its imaginary component is given a value of $h$. The derivative is then estimated as shown in equation 6

Eq. 8 $\qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad f'(x)\approx {\frac {\Im (f(x+ih))}{h}}$

Both approaches are numerical estimates, however, as will be shown, CTSE provides far superior results.

As we shall see, $h<<{\hat {h}}$, that is, the step size used to perturb the function for CTSE ($h$) is many orders of magnitude smaller than the step size ${\hat {h}}$ used in the finite difference method.

There is one major algorithmic difference between FD and CTSE. Comparing equations 7 and 8, notice that the estimate of the derivative using FD is the difference between two results, whereas CTSE has no differencing operation. As a result, FD has subtractive error if ${\hat {h}}$ is too small. Conversely, truncation error results if ${\hat {h}}$ is too large. Consequently, there is a sweet spot for the step size, i.e. ${\hat {h}}$ cannot be too large nor can it be too small - the bounds on the step size are two sided. In addition, a satisfactory step size is parameter and problem dependent - what works in one situation will not necessarily work in another. Finally, there usually is no way to know when a good step size has been reached.

CTSE also suffers from truncation error but not subtractive error - this lack of subtractive error is the significant difference between the two methods and the reason why CTSE is numerically superior. Hence, the bounds on the step size for CTSE are one sided - $h$ can be too large but it cannot be too small (within the representation of a real number by the computer). As a result, the step size $h$ can easily be made sufficiently small such that truncation error is eliminated without suffering from subtractive error. In summary, both FD and CTSE suffer from truncation error - a large step size leads to errors; however, only FD suffers from subtractive error - a small step size leads to errors.

The one drawback to CTSE is that the function $f$ must be complex valued and the intrinsic real variables within the computer model must become defined as complex variable. Hence CTSE is an intrusive method - the source code must be changed; however, the change is simple in concept. Often, this leads to a conceptual challenge - what is the physical meaning of an imaginary displacement, force, stress, strain, or pressure, as examples? The answer is that the imaginary component is merely a mathematical convenience for computing the derivative and of no physical significance.

Figure 1: FD & CTSE perturbation

As an aside, you may be wondering if a negative step is useful in CTSE? A negative step in finite differencing results in the “backwards difference” method, which sometimes has advantages. That is, when the step sizes are of significant size, as in the finite difference method, the a backwards step will sometimes yield an improved estimate over the forward step. However, in the case of CTSE, although a negative step is certainly possible, there no value in adopting a negative step since the step size is so small, that is, forwards and backwards step results yield no difference. Hence, a forward step will be used for all examples. However, the mathematical equations for a backwards step are listed below.

The Taylor series for a negative step is

Eq. 9 $\qquad \qquad \qquad \qquad f(x-i\,h)=f(x)-ihf'(x)+{\frac {h^{2}}{2}}f''(x)+{\frac {ih^{3}}{6}}f'''(x)+\cdots$

Therefore, the derivative can be determined from the equation

Eq. 10 $\qquad \qquad \qquad \qquad \qquad \qquad \qquad ~~f'(x)\approx {\frac {\Im (f(x-ih))}{(-h)}}$

Although the real power of CTSE lies in its application to sophisticated numerical models such as large finite element analyses, it is instructive to observe how it operates on several common functions. Therefore, analytical examples are provided followed by numerical examples. In all cases, the numerical result is within machine precision. Of course, since all computer languages contain complex numbers, the user can easily verify CTSE on these functions and many others.

2 Closed-form examples of elementary functions

It is instructive to apply CTSE analytically to elementary functions to assess their behavior as $h\rightarrow 0$. In these examples, it is straightforward to determine the behavior as $h\rightarrow 0$; however, many computer languages have a “Limit” operator such that the limit as $h\rightarrow 0$ can be implemented in software. With each example is a figure that is plot of the real and imaginary portions of the function. It is also shown that CTSE automatically accounts for the product and chain rules from calculus. These rules do not have to be explicitly implemented.

Example 2.1:

$\qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad f(x)=x^{2}$

Consider the simple function $f(x)=x^{2}$. Then

$\qquad \qquad \qquad \qquad \qquad \qquad \quad {\begin{aligned}f(x+ih)&=(x+ih)^{2}=x^{2}+2\,ih\,x-h^{2}\\{\frac {df}{dx}}&={\frac {\Im \left(f(x+ih)\right)}{h}}\\&={\frac {2h\,x}{h}}=2x\end{aligned}}$

We see in that for this example CTSE is exact as there are no higher order terms with respect to $h$.

Figure 2: Graphical derivative of $x^{2}$ and $dx^{2}/dx$

Example 2.2:

$\qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad f(x)=x^{3}$

Consider the function $f(x)=x^{3}$. Then

$\qquad \qquad \qquad \qquad \qquad ~{\begin{aligned}f(x+ih)&=(x+ih)^{3}=(x^{3}-3h^{2}x)+i(3hx^{2}-h^{3})\\{\frac {df}{dx}}&\approx {\frac {\Im \left(f(x+ih)\right)}{h}}={\frac {3hx^{2}-h^{3}}{h}}=3x^{2}-h^{2}\end{aligned}}$

Here the derivative is exact to within $O(h^{2})$. However, the exact derivative is recovered as $h\to 0$. In effect,

$\qquad \qquad \qquad \qquad \qquad \qquad \qquad \quad {\frac {df}{dx}}=\lim _{h\to 0}{\frac {\Im ((x+ih)^{3})}{h}}=3x^{2}$

Figure 3: Graphical derivative of $x^{3}$ and $dx^{3}/dx$

Example 2.3:

$\qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad f(x)=e^{x}$

Consider $f(x)=e^{x}$. Then

Eq. 11 $\qquad \qquad \qquad {\begin{aligned}f(x+ih)&=\exp(x+ih)=\exp(x)\exp(ih)=\exp(x)\left(\cos(h)+i\sin(h)\right)\\\\{\frac {df}{dx}}&\approx {\frac {\Im \left(f(x+ih)\right)}{h}}={\frac {\exp(x)\sin(h)}{h}}\end{aligned}}$

As shown in Eq. 11, CTSE is not exact but as $h\to 0$, $\sin(h)/h\rightarrow 1$ and the numerical derivative is recovered as $h$ approaches zero. In effect,

$\qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \quad {\frac {df}{dx}}=\lim _{h\to 0}{\frac {\Im (e^{x+ih})}{h}}=e^{x}$

Figure 4: Graphical derivative of $e^{x}$ and $de^{x}/dx$

Example 2.4:

$\qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \quad ~f(x)=\sin(x)$

Consider $f(x)=\sin(x)$. Using the identity for the sum of two angles

$\qquad \qquad \qquad \qquad \qquad \qquad \quad \sin(\alpha +\beta )=\sin(\alpha )\cos(\beta )+\cos(\alpha )\sin(\beta )$

we can express the sine of a complex number $x+ih$ as

$\qquad \qquad \qquad \qquad \qquad \qquad ~~\sin(x+ih)=\sin(x)\cos(ih)+\cos(x)\sin(ih)$

using $\sin(ih)=i\sinh(h)$, and $\cos(ih)=\cosh(h)$

$\qquad \qquad \qquad \qquad \qquad \qquad \sin(x+ih)=\sin(x)\cosh(h)+i\cos(x)\sinh(h)$

Then

$\qquad \qquad \qquad \qquad \qquad \qquad \quad ~{\frac {df}{dx}}\approx {\frac {\Im (\sin(x+ih))}{h}}={\frac {\cos(x)\sinh(h)}{h}}$

As $h\to 0$, $\sinh(h)/h\to 1$ and the numerical derivative is recovered as $h$ approaches zero. In effect,

$\qquad \qquad \qquad \qquad \qquad \qquad \qquad {\frac {df}{dx}}=\lim _{h\to 0}{\frac {\Im (\sin(x+ih))}{h}}=\cos(x)$

Figure 5: Graphical derivative of $\sin(x)$ and $d\sin(x)/dx$

Example 2.5:

$\qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \quad ~f(x)=\cos(x)$

Consider $f(x)=\cos(x)$. Using the identity for the sum of two angles

$\qquad \qquad \qquad \qquad \qquad \qquad \cos(\alpha +\beta )=\cos(\alpha )\cos(\beta )-\sin(\alpha )\sin(\beta )$

we can express the cosine of a complex number $x+ih$ as

$\qquad \qquad \qquad \qquad \qquad \quad {\begin{aligned}\cos(x+ih)&=\cos(x)\cos(ih)-\sin(x)\sin(ih)\\&=\cos(x)\cosh(h)-i\sin(x)\sinh(h)\end{aligned}}$

Then

$\qquad \qquad \qquad \qquad \qquad \qquad \qquad \quad {\frac {\Im (\cos(x+ih)}{h}}={\frac {-\sin(x)\sinh(h)}{h}}$

As $h\to 0,\,\sinh(h)/h\to 1$. In effect,

$\qquad \qquad \qquad \qquad \qquad \qquad \quad {\frac {df}{dx}}=\lim _{h\to 0}{\frac {\Im (\cos(x+ih))}{h}}=-\sin(x)$

Figure 6: Graphical derivative of $\cos(x)$ and $d\cos(x)/dx$

Example 2.6:

$\qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \quad ~f(x)=\tan(x)$

Consider $f(x)=\tan(x)$. Using the identity for the sum of two angles

$\qquad \qquad \qquad \qquad \quad \tan(\alpha +\beta )={\frac {\sin(\alpha +\beta )}{\cos(\alpha +\beta )}}={\frac {\sin(\alpha )\cos(\beta )+\cos(\alpha )\sin(\beta )}{\cos(\alpha )\cos(\beta )-\sin(\alpha )\sin(\beta )}}$

we can express the tangent of the complex number $x+ih$ as

$\qquad \qquad \qquad \qquad \qquad \quad {\begin{aligned}\tan(x+ih)&={\frac {\sin(x)\cos(ih)+\cos(x)\sin(ih)}{\cos(x)\cos(ih)-\sin(x)\sin(ih)}}\\&={\frac {\sin(x)\cosh(h)+i\cos(x)\sinh(h)}{\cos(x)\cosh(h)-i\sin(x)\sinh(h)}}\end{aligned}}$

Multiplying the numerator and the denominator by the conjugate of the denominator gives:

$\quad {\begin{aligned}\tan(x+ih)&={\frac {\sin(x)\cos(x)(\cosh ^{2}(h)-\sinh ^{2}(h))+i\cosh(h)\sinh(h)(\cos ^{2}(x)+\sin ^{2}(x))}{\left(\cos(x)\cosh(h)\right)^{2}+\left(\sin(x)\sinh(h)\right)^{2}}}\\&={\frac {\sin(x)\cos(x)+i\cosh(h)\sinh(h)}{\left(\cos(x)\cosh(h)\right)^{2}+\left(\sin(x)\sinh(h)\right)^{2}}}\end{aligned}}$

Then

$\qquad \qquad \qquad \quad {\frac {df}{dx}}\approx {\frac {\Im (\tan(x+ih))}{h}}={\frac {\sinh(h)\cosh(h)}{h((\cos(x)\cosh(h))^{2}+(\sin(x)\sinh(h))^{2})}}$

As $h\to 0$, $\sinh(h)/h\to 1$. In effect,

$\qquad \qquad \qquad \qquad \quad ~~~~~{\begin{aligned}{\frac {df}{dx}}=\lim _{h\to 0}{\frac {\Im (\tan(x+ih))}{h}}&={\frac {1}{((1)\cos(x))^{2}+((0)\sin(x))^{2}}}\\&={\frac {1}{\cos ^{2}(x)}}\\&=\sec ^{2}(x)\end{aligned}}$

Figure 7: Graphical derivative of $\tan(x)$ and $d\tan(x)/dx$

Example 2.7:

$\qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \quad ~f(x)=\sinh(x)$

Consider $f(x)=\sinh(x)$. Then using the complex variable version of the hyperbolic sine

$\qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \quad ~~\sinh(z)={\frac {1}{2}}(e^{z}-e^{-z})$

$\qquad \qquad \qquad \qquad {\begin{aligned}\sinh(x+ih)&={\frac {1}{2}}\left(e^{x+ih}-e^{-x-ih}\right)={\frac {1}{2}}\left(e^{x}e^{ih}-e^{-x}e^{-ih}\right)\\&={\frac {1}{2}}e^{x}\left(\cos(h)+i\sin(h)\right)-{\frac {1}{2}}e^{-x}\left(\cos(h)-i\sin(h)\right)\\&={\frac {1}{2}}\cos(h)\left(e^{x}-e^{-x}\right)+i{\frac {1}{2}}\sin(h)\left(e^{x}+e^{-x}\right)\\&=\cos(h)\sinh(x)+i\sin(h)\cosh(x)\end{aligned}}$

Then,

$\qquad \qquad \qquad \qquad \qquad \qquad \qquad \quad {\frac {\Im (\sinh(x+ih))}{h}}={\frac {\cosh(x)\sin(h)}{h}}$

As $h\to 0$, $\sinh(h)/h\to 1$ and the numerical derivative is recovered as $h$ approaches zero. In effect,

$\qquad \qquad \qquad \qquad \qquad \qquad \qquad ~~~{\frac {df}{dx}}={\frac {\Im (\sinh(x+ih))}{h}}=\cosh(x)$

Figure 8: Plot of $\sinh(x)$ and $d\sinh(x)/dx$

Example 2.8:

$\qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad f(x)=\cosh(x)$

Consider $f(x)=\cosh(x)$. Then using the complex variable version of the hyperbolic cosine $\qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \quad ~~\cosh(z)={\frac {1}{2}}(e^{z}+e^{-z})$

$\qquad \qquad \qquad {\begin{aligned}\cosh(x+ih)&={\frac {1}{2}}\left(e^{x+ih}+e^{-x-ih}\right)={\frac {1}{2}}\left(e^{x}e^{ih}+e^{-x}e^{-ih}\right)\\&={\frac {1}{2}}e^{x}\left(\cos(h)+i\sin(h)\right)+{\frac {1}{2}}e^{-x}\left(\cos(h)-i\sin(h)\right)\\&={\frac {1}{2}}\cos(h)\left(e^{x}+e^{-x}\right)+i{\frac {1}{2}}\sin(h)\left(e^{x}-e^{-x}\right)\\&=\cos(h)\cosh(x)+i\sin(h)\sinh(x)\end{aligned}}$

Taking

$\qquad \qquad \qquad \qquad \qquad \qquad \qquad ~{\frac {\Im (\cosh(x+ih))}{h}}={\frac {\sinh(x)\sin(h)}{h}}$

As $h\to 0$, $\sinh(h)/h\to 1$ and the numerical derivative is recovered as $h$ approaches zero. In effect,

$\qquad \qquad \qquad \qquad \qquad \qquad \qquad \quad {\frac {df}{dx}}={\frac {\Im (\cosh x+ih)}{h}}=\sinh(x)$

Figure 9: Plot of $\cosh(x)$ and $d\cosh(x)/dx$

Example 2.9:

$\qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad f(x)=ln(x)$

Consider the function $f(x)=\log(x)$ with $x>0$. Expressing the complex number $x+ih$ in its polar form.

$\qquad \qquad \qquad \qquad \qquad \qquad \qquad \quad \log(x+ih)=\log(re^{i\phi })=\log(r)+i\phi$

Taking the imaginary part and dividing by $h$

$\qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad {\frac {\Im (\log(x+ih))}{h}}={\frac {\phi }{h}}$

As $h\to 0$

$\qquad \qquad \qquad \qquad \qquad \quad ~\lim _{h\to 0}{\frac {\Im (\log(x+ih))}{h}}=\lim _{h\to 0}{\frac {\phi }{h}}=\lim _{h\to 0}{\frac {1}{h}}\tan ^{-1}({\frac {h}{x}})$

This form is an indeterminate form of the type $0/0$. Applying L'Hôpital's rule we obtain,

$\qquad \qquad \qquad \qquad \quad {\begin{aligned}\lim _{h\to 0}{\frac {\Im (\log(x+ih))}{h}}&=\lim _{h\to 0}{\frac {{\frac {d}{dh}}\tan ^{-1}({\frac {h}{x}})}{{\frac {d}{dh}}h}}=\lim _{h\to 0}{\frac {{\frac {d}{dh}}({\frac {h}{x}})}{1+({\frac {h}{x}})^{2}}}\\{\frac {df}{dx}}&=\lim _{h\to 0}{\frac {x}{x^{2}+h^{2}}}={\frac {1}{x}}\end{aligned}}$

Figure 10: Plot of $\ln(x)$ and $d\ \ln(x)/dx$

Example 2.10

$\qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad f(x)={\sqrt {x}}$

Consider the function $f(x)={\sqrt {x}}$ with $x>0$. Expressing $z=x+ih$ in its polar form,

$\qquad \qquad \qquad \qquad \qquad \qquad ~~{\sqrt {z}}=r^{1/2}e^{i\theta /2}=r^{1/2}(\cos(\theta /2)+i\sin(\theta /2))$

Using the fact that $r={\sqrt {x^{2}+h^{2}}}$ and for small $h$ that $\sin(\theta /2)={\frac {h/2}{r}}$,

$\qquad \qquad \qquad \qquad {\frac {df}{dx}}\approx \Im (f(x+ih)/h=r^{1/2}{\frac {h/2}{rh}}={\frac {r^{1/2}}{2r}}={\frac {1}{2r^{1/2}}}={\frac {1}{2{\sqrt {x^{2}+h^{2}}}}}$

As $h\rightarrow 0$,

$\qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad {\frac {df}{dx}}=\lim _{h\to 0}{\frac {\Im ({\sqrt {x+ih}})}{h}}={\frac {1}{2{\sqrt {x}}}}$

Figure 11: Plot of ${\sqrt {x}}$ and $d{\sqrt {x}}/dx$

Example 2.11

$\qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \quad ~~f(x)=xe^{x}$

Consider the function $f(x)=xe^{x}$. Replacing $x$ by $x+ih$, (notice that we need to replace all occurrences of $x$ with $x+ih$)

$\qquad \qquad \qquad \qquad {\begin{aligned}(x+ih)e^{(x+ih)}&=(x+ih)e^{x}e^{ih}=(x+ih)e^{x}(cos(h)+isin(h))\\&=e^{x}(x\cos(x)-h\sin(h))+e^{x}(x\sin(h)+h\cos(h))i\\{\frac {df}{dx}}&\approx {\frac {{\text{Im}}\left(f(x+ih)\right)}{h}}=e^{x}(xsin(h)/h+cos(h))\end{aligned}}$

As $h\to 0$, $sin(h)/h\to 1$, and $cos(h)\to 1$, hence,

$\qquad \qquad \qquad \qquad \qquad \qquad \qquad \quad {\frac {df}{dx}}=\lim _{h\to 0}(x+ih)e^{x+ih}=e^{x}(x+1)$

and the exact derivative is obtained.

Figure 12: Plot of $f(x)=xe^{x}$ and $d\ xe^{x}/dx$

Example 2.12

$\qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \quad f(x)={\frac {sin(x)}{x}}$

Consider the function $f(x)={\frac {sin(x)}{x}}$. Replacing $x$ by $x+ih$,

$\qquad \qquad \qquad \qquad \qquad \qquad {\frac {sin(x+ih)}{x+ih}}={\frac {cosh(h)sin(x)+icos(x)sinh(h)}{x+ih}}$

Multiplying the numerator and denominator by $x-ih$ yields

$\qquad \qquad \qquad \qquad \quad {\begin{aligned}{\frac {sin(x+ih)(x-ih)}{x+h^{2}}}&={\frac {(xcosh(h)sin(x)+hcos(x)sinh(h))}{x+h^{2}}}\\&+{\frac {i(-cosh(h)sin(x)h+xcos(x)sinh(h))}{x+h^{2}}}\end{aligned}}$

Extracting the imaginary part and dividing by $h$ yields

$\qquad \qquad \qquad \qquad \qquad {\frac {df}{dx}}\approx \Im (f(x+ih)/h={\frac {-cosh(h)sin(x)+xcos(x)sinh(h)/h}{x+h^{2}}}$

In the limit as $h\rightarrow 0$, $cosh(h)\rightarrow 1$, $sinh(h)/h\rightarrow 1$, and $h^{2}\rightarrow 0$. As a result,

$~~~\Im (f(x+ih)/h=f'(x)={\frac {-sin(x)+xcos(x)}{x}}.{\frac {df}{dx}}=\lim _{h\to 0}{\frac {sin(x+ih)}{x+ih}}={\frac {-sin(x)+xcos(x)}{x}}$

Figure 13: Plot of ${\frac {sin(x)}{x}}$ and its derivative

Example 2.13

$\qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad f(x)=e^{x^{2}}$

Consider the function $f(x)=e^{x^{2}}$. Replacing $x$ by $x+ih$,

$\qquad \qquad \qquad e^{(x+ih)^{2}}=e^{(x^{2}+2hxi-h^{2})}=e^{(x^{2}-h^{2})}e^{2hxi}=e^{(x^{2}-h^{2})}(cos(2hx)+isin(2hx))$

Extracting the imaginary part and dividing by $h$ yields

$\qquad \qquad \qquad \qquad \qquad \qquad \qquad {\frac {df}{dx}}\approx \Im (f(x+ih)/h=e^{(x^{2}-h^{2})}{\frac {sin(2hx)}{h}}$

In the limit as $h\rightarrow 0$, $\sin(2xh)/h\rightarrow 2x$, and $h^{2}\rightarrow 0$.

Note, the first few terms of the Taylor series for

$\qquad \qquad \qquad \qquad \qquad \qquad \qquad \quad \sin(2hx)=2hx-4h^{3}x^{3}/3+O(x^{5})$,

hence,

$\qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad {\underset {h\rightarrow 0}{\lim }}\sin(2hx)/h=2x$.

As a result,

$\qquad \qquad \qquad \qquad \qquad \qquad \qquad \quad {\frac {df}{dx}}={\underset {h\to 0}{\lim }}\ e^{(x^{2}-h^{2})}{\frac {sin(2hx)}{h}}=2xe^{x^{2}}$

Figure 14: Plot of $e^{x^{2}}$ and its derivative

3 Cauchy-Riemann (CR) form

Although it is almost always possible to use complex algebra within any given situation, it is valuable to understand the Cauchy-Riemann form of complex numbers and how to use this form for computing derivatives. One example would be where one needs to solve a complex system of equations but only has access to a real-variable solver. In this case, one can expand the complex system into a larger but real-only system of equations, then apply the real-variable solver.

Any complex number can be written as a matrix of all real numbers using the Cauchy-Riemann matrix form, i.e.

Eq. 12$\qquad \qquad \qquad \qquad \qquad \qquad \qquad a+bi=\left[{\begin{array}{cc}a&-b\\b&a\end{array}}\right]$

In the CR matrix, the real part is located along the diagonal and the imaginary part is located on the off diagonals. Note, that the matrix is non-symmetric.

Thus, to differentiate a function $f$ using CTSE, one can compute $f(x+ih)=f\left[{\begin{array}{cc}x&-h\\h&x\end{array}}\right]$, where $h$ is the small step size. After evaluation of $f\left[{\begin{array}{cc}x&-h\\h&x\end{array}}\right]$, the result will be of the form $\left[{\begin{array}{cc}f(x)&-f'(x)h\\f'(x)h&f(x)\end{array}}\right]$. Hence, the derivative will be located in the 2nd row, 1st column, times the step size $h$. Thus, the equivalent form of the equation $f'(x)=\Im (f(x+ih)/h$ is to extract the value from the 2nd row, 1st column, and divide by $h$. Several examples that demonstrate how to compute the derivative using the CR matrix form of a complex number are provided below. Note, however, that functions of matrices are often required. For example, $\sin \left[{\begin{array}{cc}x&-h\\h&x\end{array}}\right]$ must use the matrix sine function.

Example 3.1:

$\qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad f(x)=x^{2}$

Consider the simple function $f(x)=x^{2}$.

$\qquad \qquad \qquad \qquad \qquad \qquad \qquad f(x)=\left[{\begin{array}{cc}x&-h\\h&x\end{array}}\right]^{2}=\left[{\begin{array}{cc}x^{2}-h^{2}&-2hx\\2hx&x^{2}-h^{2}\end{array}}\right]$

Hence, $f'(x)={\frac {2hx}{h}}=2x$.

Example 3.2:

$\qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad f(x)=x^{3}$

Consider the simple function $f(x)=x^{3}$.

$\qquad \qquad \qquad \qquad f(x)=\left[{\begin{array}{cc}x&-h\\h&x\end{array}}\right]^{3}=\left[{\begin{array}{cc}x\left(x^{2}-h^{2}\right)-2h^{2}x&-h\left(x^{2}-h^{2}\right)-2hx^{2}\\h\left(x^{2}-h^{2}\right)+2hx^{2}&x\left(x^{2}-h^{2}\right)-2h^{2}x\end{array}}\right]$

Hence, $f'(x)={\frac {h\left(x^{2}-h^{2}\right)+2hx^{2}}{h}}=\left(x^{2}-h^{2}\right)+2x^{2}$. The exact derivative is recovered as $h\rightarrow 0$.

Example 3.3:

$\qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad f(x)=e^{x}$

Consider the simple function $f(x)=e^{x}.\,\ f(x)=expM\left[{\begin{array}{cc}e^{x}&-h\\h&e^{x}\end{array}}\right]=\left[{\begin{array}{cc}e^{x}\cos(h)&-e^{x}\sin(h)\\e^{x}\sin(h)&e^{x}\cos(h)\end{array}}\right]$. Here we use $expM$ to indicate the matrix exponential function. The derivative is determined as, $f'(x)={\frac {1}{h}}e^{x}\sin(h)$. The exact derivative is recovered as $h\rightarrow 0$.

4 Numerical step size study

The comparison of the effect of the step size on the accuracy of the derivative can be shown through the following example. Consider the single variable function shown below and its analytical derivative evaluated at $x=1.5$, taken from Reference [ref Fike]. The exact values given up to 15 significant figures are $f(1.5)=4.49778005394616$ and $4.05342789389862.$

Eq. 13 $\qquad \qquad \qquad \qquad \qquad \qquad f(x)={\frac {e^{x}}{\sqrt {\sin ^{3}(x)+\cos ^{3}(x)}}}$

$\qquad \qquad \qquad \qquad \quad ~~f'(x)={\frac {e^{x}(9\sin(x)+\sin(3x)+3\cos(x)+5\cos(3x))}{8\left(\sin ^{3}(x)+\cos ^{3}(x)\right)^{3/2}}}$

The derivative of this function is computed using the forward difference (FD) formula,$f'(x)\approx {\frac {f(x+h)-f(x)}{h}}$, the central differencing (CD) formula, $f'(x)\approx {\frac {f(x+h)-f(x-h)}{2h}}$, and CTSE, each as a function of step size. The results are plotted in below with all computations computed as double precision with 15 significant digits. The relative error is defined as $Relerr={\frac {f'_{Numerical}(1.5)-f'_{Exact}(1/5)}{f'_{Exact}(1.5)}}$, where $f'_{Numerical}(1.5)$ is the derivative computed by FD, CD, or CTSE. As shown, the accuracy of the forward differencing method increases as the step size is initially reduced due to a reduction in truncation error; however, for step sizes below $10^{-8}$, the error increases due to subtraction error dominating the truncation error. Similarly, the central differencing method shows improved accuracy until a step size of $10^{-5}$ has been reached, beyond which the error increases. As a result, one can see that neither the forward nor central differencing methods can achieve machine precision for any step size. In addition, this “best” step size is problem and parameter dependent, that is, there is no straightforward way to know a priori the appropriate step size, (although an optimum step size has been proposed that requires an approximation to the second order derivative [Haftka]). These methods behave as they do due to the conflict between truncation error (smaller step size is better), and subtractive error (larger step size is better). CTSE only suffers from truncation error, hence, the smaller the step size the better the accuracy, with no limit until machine precision is

Figure 15: Accuracy of forward difference (blue), central difference (purple), and CTSE (red) as a function of step size for the function $f(x)=e^{x}{sin^{3}(x)+cos^{3}(x)}$ evaluated at $x=1.5$

reached. In fact, a step size of $10^{-8}$ the relative error, defined as

$RelErr={\frac {(f'_{CTSE}-f'_{Exact})}{f'{}^{Exact}}}$, is less than $10^{-15}$, which is machine precision. (Note, the relative error was computed as zero for some steps sizes. As zero cannot be plotted in a logarithmic plot, a value of 1E-16 was used.) In section 1, we discussed why the step size $h$ must be greater than $10^{-308}$. In fact, one can continue the step size study all the way until $h=10^{-308}$ and verify that the relative error will be approximately $10^{-15}$. Note, there is nothing exceptional about the function chosen for the step size study. Almost identical results will be obtained for any other function.

5 Numerical evaluation of special functions

It is quite convenient to use CTSE to compute the derivative for any number of numerical functions such as the Error and complementary Error functions (Erf, Erfc), the Gamma function, the Beta function, the Mathieu function, inverse trigonometric functions, elliptic integrals, Airy functions and many others, as these functions have complex variable versions.

One can see CTSE in action by plotting $\Re (f(x+ih))$, and $\Im (f(x+ih)/h)$ as a function of $x$ for several special functions as shown below. As expected, the function and its derivative are contained in the plot. Several numerical examples are provided.

Figure 16: Plot of $f(x)=\sin ^{-1}(x)$ and $f(x)=d\sin ^{-1}(x)/dx$	Figure 17: Plot of $f(x)=Erf(x)$ and dErf(x)/dx
Figure 18: Plot of the complementary error function and its derivative	Figure 19: Plot of Gamma function and its derivative
Figure 20: Plot of the Bessel function of the first kind $J_{0}$ and its derivative	Figure 21: Plot of the elliptic integral of the first kind and its derivative
Figure 22: Plot of Beta function and its derivative	Figure 23: Plot of the Airy function and its derivative

6 Finite element example using a beam element

A simple finite element example using an Euler-Bernoulli beam element is shown in Figure 23. It is straightforward to apply CTSE to compute derivatives with respect the model parameters and obtain exact results as $h\to 0$.

The beam is under load $P$ and has displacement, $\delta$, and rotation, $\theta$, degrees of freedom at each node. The stiffness matrix for a beam element of elastic modulus $E$, cross sectional moment of inertia $I$, and length $L$, with displacement, $\delta$, and rotation, $\theta$, degrees of freedom at each node is

$\qquad \qquad \qquad \qquad \qquad \qquad \quad ~~~K={\frac {EI}{L^{3}}}\left[{\begin{array}{cccc}12&6L&-12&6L\\6L&4L^{2}&-6L&2L^{2}\\-12&-6L&12&-6L\\6L&2L^{2}&-6L&4L^{2}\end{array}}\right]$

After applying the constraints at the right end $(\delta =\theta =0)$, the system of equations becomes

$\qquad \qquad \qquad \qquad \qquad \qquad \qquad ~~~\left\{{\begin{array}{c}P\\0\end{array}}\right\}={\frac {EI}{L^{3}}}\left[{\begin{array}{cc}12&6L\\6L&4L^{2}\end{array}}\right]\left\{{\begin{array}{c}\delta \\\theta \end{array}}\right\}$

The solution to this system of equations can be obtained as

$\qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \quad ~~\left\{{\begin{array}{c}\delta \\\theta \end{array}}\right\}=K^{-1}\left\{{\begin{array}{c}P\\0\end{array}}\right\}$

yielding

$\qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \quad \delta ={\frac {PL^{3}}{3EI}},\,\theta =-{\frac {PL^{2}}{2EI}}$

6.1 Sensitivity with respect to loading P

Consider the computation of the derivative of the displacement and rotation with respect to the applied load $P$. The exact solutions are

$\qquad \qquad \qquad \qquad \qquad \qquad \qquad \quad \quad ~{\frac {d\delta }{dP}}={\frac {L^{3}}{3EI}}\operatorname {;} {\frac {d\theta }{dP}}=-{\frac {L^{2}}{2EI}}$

Based on the CTSE method, we perturb the load $P$ by $ih$ and solve the system of equations

$\qquad \qquad \qquad \qquad \qquad \qquad \qquad \quad \left\{{\begin{array}{c}P+ih\\0\end{array}}\right\}={\frac {EI}{L^{3}}}\left[{\begin{array}{cc}12&6L\\6L&4L^{2}\end{array}}\right]\left\{{\begin{array}{c}\delta \\\theta \end{array}}\right\}$

The solution is

$\qquad \qquad \qquad \qquad \qquad \qquad \quad \delta ={\frac {PL^{3}}{3EI}}+ih{\frac {L^{3}}{3EI}}\ ;\ \theta =-{\frac {PL^{2}}{2EI}}-ih{\frac {L^{2}}{2EI}}$

The derivatives of the displacement and rotation with respect to the applied load $P$ are obtained exactly as

$\qquad \qquad \qquad \qquad ~{\frac {\partial \delta }{\partial P}}=\Im (\delta (P+ih))/h={\frac {L^{3}}{3EI}}\ ;\ {\frac {\partial \theta }{\partial P}}=\Im (\theta (P+ih))/h=-{\frac {L^{2}}{2EI}}$

6.2 Sensitivity with respect to elastic modulus E

Consider the computation of the derivative of the displacement and rotation with respect to the elastic modulus. The exact solutions are

$\qquad \qquad \qquad \qquad \qquad \qquad \qquad \quad {\frac {d\delta }{dE}}=-{\frac {PL^{3}}{3E^{2}I}}\ ;\ {\frac {d\theta }{dE}}={\frac {PL^{2}}{2E^{2}I}}$

Based on the CTSE method, we perturb the load elastic modulus by $ih$ and solve the system of equations

$\qquad \qquad \qquad \qquad \qquad \qquad \qquad \left\{{\begin{array}{c}P\\0\end{array}}\right\}={\frac {(E+ih)I}{L^{3}}}\left[{\begin{array}{cc}12&6L\\6L&4L^{2}\end{array}}\right]\left\{{\begin{array}{c}\delta \\\theta \end{array}}\right\}$

The solutions are

$\qquad \qquad \qquad \qquad \qquad \qquad \quad ~~\delta ={\frac {PL^{3}}{3EI+i(h3I)}}\ ;\ \theta =-{\frac {PL^{2}}{2EI+i(h2I)}}$

Applying the complex conjugate to $\delta$ and $\theta$, the derivatives with respect to the elastic modulus $E$ are

$\qquad \quad {\frac {\partial \delta }{\partial E}}\approx Im(\delta (E+ih))/h=-{\frac {PL^{3}}{3I(E^{2}+h^{2})}}\ ;\ {\frac {\partial \theta }{\partial E}}\approx Im(\theta (E+ih))/h={\frac {PL^{2}}{2I(E^{2}+h^{2})}}$

We see that as $h\to 0$, the exact results are obtained.

6.3 Sensitivity with respect to length L

Consider the computation of the derivative of the displacement and rotation with respect to the length of the beam. The exact solutions are

$\qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad ~{\frac {d\delta }{dL}}={\frac {PL^{2}}{EI}};{\frac {d\theta }{dL}}=-{\frac {PL}{EI}}$

Based on the CTSE method, to compute the sensitivity with respect to the element length, we perturb the length by $ih$ and solve the system of equations

$\qquad \qquad \qquad \qquad \qquad \quad \left\{{\begin{array}{c}P\\0\end{array}}\right\}={\frac {EI}{(L+ih)^{3}}}\left[{\begin{array}{cc}12&6(L+ih)\\6(L+ih)&4(L+ih)^{2}\end{array}}\right]\left\{{\begin{array}{c}\delta \\\theta \end{array}}\right\}$

The solution to this system of equations is

$\qquad \qquad \qquad \delta =\left({\frac {PL^{3}}{3EI}}-{\frac {Ph^{2}}{EI}}\right)+i\left({\frac {PL^{2}h}{EI}}-{\frac {Ph^{3}}{3EI}}\right)\ ;\ \theta =\left({\frac {PL^{2}}{2EI}}+{\frac {Ph^{2}}{2EI}}\right)-i\left({\frac {PLh}{EI}}\right)$

The derivatives of the displacement and rotation with respect to the beam length are approximated as

$\qquad \qquad \qquad {\frac {\partial \delta }{\partial L}}\approx \Im (\delta (L+ih))/h={\frac {PL^{2}}{EI}}-{\frac {Ph^{2}}{3EI}}\ ;\ {\frac {\partial \theta }{\partial L}}\approx \Im (\theta (L+ih))/h=-{\frac {PL}{EI}}$

We see that the displacement derivative is exact up to $O(h^{2})$ and the derivative of the rotation is exact. We see that as $h\to 0$, the exact results are obtained.

Figure 24: Schematic of cantilever beam under load P

7 Numerical integration

CTSE can be applied within numerical integration algorithms to compute the derivatives with respect to parameters within the integrand and the integration limits. The following examples demonstrate this using two numerical integration methods, Simpson's rule and Gauss-Legendre quadrature; however, other numerical integration algorithms can be used in an equally straightforward way. The important point to realize is that no changes are required to the numerical algorithm. In effect, CTSE can be implemented unaltered with the exception of perturbing one of the parameters along the imaginary axis by a small step size $h$.

7.1 Simpson's rule

The particular Simpson's rule used for the demonstration is Simpson's 1/3 rule, defined as

Eq. 14 $\qquad \qquad I(a,b)=\int \limits _{1}^{3}f(x)dx\approx {\frac {b-a}{8}}[f(a)+3f({\frac {2a+b}{3}})+3f({\frac {a+2b}{3}})+f(b)]$

where $f$ is the function to be integrated.

To compute derivatives with respect to parameters $a,b$, we perturb the parameter of interest, compute the integral, extract the imaginary component, then divide by the step size. To be more explicit, to compute $dI/da$, we define $a^{*}=a+ih$ and use this value in the formula for Simpson's rule, that is

Eq. 15 $dI(a,b)/da\approx {\frac {1}{h}}Im\left[{\frac {b-(a+ih)}{8}}[f(a+ih)+3f({\frac {2(a+ih)+b}{3}})+3f({\frac {(a+ih)+2b}{3}})+f(b)]\right]$

and similarly for $dI/db$. Each perturbation of $a$, and $b$ is conducted independently; hence, two analyses are required to compute derivatives with respect to $a$, and $b$. It should be clear that the integration is conducted using real variables and not integration within the complex plane.

Example 7.1.1:

$\qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad ~~~I(a,b)=\int _{a}^{b}2xdx$

Consider the integral $I(a,b)=\int _{a}^{b}2xdx$ with the analytical solutions $I(a,b)=(b^{2}-a^{2})$, $dI/db=2b,\ dI/da=-2a$. The integral $I$ and the derivatives with respect to $a$, and $b$ were computed using numerical values $a=1,\ b=2$ and perturbing each in turn by a step size of $h=10^{-10}$. For example, to compute the derivative with respect to $b$, we perform the numerical operations

Eq. 16 $dI/db\approx {\frac {1}{h}}\Im \left[{\frac {(2+ih)-1}{8}}[f(1)+3f({\frac {2(1)+(2+ih)}{3}})+3f({\frac {1+2(2+ih)}{3}})+f(2+ih)]\right]$

The results are shown in Table In this case, the evaluation of the integral and the derivatives with respect to $a$ and $b$ are accurate to within machine precision, i.e., 15 decimal places.

Notice that it is not necessary to have any special consideration when computing derivatives with respect to a parameter in the integration limits, that is, we do not need to explicitly apply Leibniz's rule for this case, that is, for example $dI/db=f(b)+\int _{a}^{b}(df(x)/db)dx=f(b)$.

$\qquad \qquad \qquad \qquad \qquad \qquad {\begin{array}{|c|c|c|}\hline \ &\operatorname {Exact} &\operatorname {Simpson's\ w\ CTSE} \\\hline Integral&3&3\\\hline dI/da&-2&-2\\\hline dI/db&4&4\\\hline \end{array}}$$\qquad {\begin{aligned}\operatorname {Table:1\ CTSE\ applied\ within\ the\ Simpson'} &\operatorname {s\ rule\ for\ numerical\ integration\ of} \int _{a}^{b}2xdx\\\operatorname {with\ values\ a=1,\ b=2.\ A\ } &\operatorname {step\ size\ of} h=10^{-10}\operatorname {was\ used.} \end{aligned}}$

Example 7.1.2:

$\qquad \quad \qquad \qquad \qquad \quad \qquad \qquad \qquad \qquad I(a,b,c)=\int _{a}^{b}xe^{-cx}dx$

Consider the integral $I(x,a,b,c)=\int _{a}^{b}xe^{-cx}dx$ where $a,\ b,$ and $c$ are constants. In this case, we will compute the derivates with respect to the limits $a$ and $b$, and a parameter $c$ contained within the integrand. The exact indefinite integral (using integration by parts) is $I(x,c)=-{\frac {e^{-cx}(cx+1)}{c^{2}}}$ and the exact definite integral is $I(a,b,c)={\frac {(ac+1)e^{-ac}-(bc+1)e^{-bc}}{c^{2}}}$. The exact derivatives of this integral with respect to the constants are: ${\frac {dI(a,b,c)}{da}}=-ae^{-ac},{\frac {dI(a,b,c)}{db}}=be^{-bc}$, and ${\frac {dI(a,b,c)}{dc}}={\frac {e^{c(-(a+b))}\left(e^{ac}\left(b^{2}c^{2}+2bc+2\right)-\left(a^{2}c^{2}+2ac+2\right)e^{bc}\right)}{c^{3}}}$.

The analytical derivatives are:

${\frac {d}{da}}\int _{a}^{b}f(x)dx=f(a)$

The computation of the derivatives $dI/da$ and $dI/db$ are computed as in the previous example, i.e., one can use equation 15 to compute $dI/db$. The derivative $dI/dc$ is conducted by applying Simpson's rule to the function $\int _{a}^{b}Im(f(x,c+ih))/hdx$ using CTSE. Simpson's rule becomes for this case

$\qquad \quad {\begin{aligned}dI/dc&=\int _{a}^{b}(df(x,c)/dc)dx\approx {\frac {1}{h}}Im(\int _{a}^{b}f(x,c+ih))dx)\\&\approx {\frac {b-a}{8h}}[f(a,c+ih)+3f({\frac {2a+b}{3}},c+ih)+3f({\frac {a+2b}{3}},c+ih)+f(b,c+ih)]\end{aligned}}$

The results are shown in Table 2 using the values $a=1,\ b=2,$ and $c=2$. These results indicate that neither the integral nor the derivatives are exact. This is an important point. Although we have shown that CTSE can compute derivatives of machine precision, this does not mean that the derivatives are always exact - it depends upon the application. In general, if the real-valued operation, Simpson's rule in this case, is not exact, neither will be the derivatives. While there is no formal proof, empirical results usually show that the derivative results are slightly less accurate than the real-valued results. For example, for this case, the relative error in the computation of the integral using Simpson's rule is $\approx 10^{-4}$, whereas the relative error in the derivatives are $\approx 10^{-3}$. Of course, improved results for the integral and the derivatives could be obtained by subdividing the integral and applying equation 14 repeatedly.

Accuracy analysis of dI/dc

A natural question arises as to the accuracy of the derivatives. It is instructive to compute the “best” estimate of $dI/dc$ that can be obtained using Simpson's rule in order to compare with the results computed using CTSE. To determine the “best” solution for $dI/dc$ we take the derivative inside the integral and analytically differentiation the integrand before integration. In this way, no error is incurred during differentiation and all error will be due to the discretization employed within Simpson's rule, for example $\qquad \qquad \qquad \quad dI/dc=\int _{a}^{b}(df(x,c)/dc)dx=-\int _{a}^{b}x^{2}e^{-cx}dx={\frac {1}{c^{3}}}e^{-cx}(c^{2}x^{2}+2cx+2)$

Applying Simpson's rule to this integrand with $a=1,\ b=2,\ c=2$ yields the result $dI/dc=-0.109557$, identical to the CTSE result shown in Table 2 - in fact, the two solutions agree to 15 decimal places. This example supports the claim that CTSE, applied with a sufficiently small step size, will provide the most accurate derivative possible and that all errors in the result are due to errors in the numerical method itself. In fact, the claim that CTSE will provide the “best” possible derivative is a general claim that is true regardless of the numerical algorithm. Figure 24 below shows the convergence of $dI/dc$ as a function of step size. The results converge at about a value of $h=10^{-5}$ and are constant thereafter. This behavior is typical when applying CTSE within numerical algorithms. While CTSE does not add any error, the accuracy of the derivatives will be dependent upon the numerical methods employed.

$\qquad \qquad \qquad \qquad \quad {\begin{array}{|c|c|c|c|}\hline &\operatorname {Exact} &\operatorname {Simpson's\ w\ CTSE} &\operatorname {Relative\ Error} \\\hline \operatorname {Integral} &0.0786069&0.0785338&9.3*10^{-4}\\\hline dI/da&-0.135335&-0.134807&3.9*10^{-3}\\\hline dI/db&0.0366313&0.0363799&6.9*10^{-3}\\\hline dI/dc&-0.109643&-0.109557&7.8*10^{-3}\\\hline \end{array}}$

${\begin{aligned}\operatorname {Table\ 2:CTSE\ applied\ to\ a\ parameter\ within\ the} &\operatorname {\ function\ using\ Simpson's\ rule\ for\ the} \\\operatorname {numerical\ integration\ of} \int _{a}^{b}xe^{-cx}dx\operatorname {with\ values\ a=1,\ b} &\operatorname {=2,\ c=2.\ A\ step\ size\ of} h=10^{-10}\operatorname {was\ used.} \end{aligned}}$

Figure 25: Step size study of the convergence of dI/dc as a function of step size using CTSE employed within Simpson's rule

7.2 Gauss-Legendre Quadrature

Gauss-Legendre quadrature is a popular numerical integration scheme. In this case the integral is computed using the formula

$\qquad \qquad \qquad \qquad \qquad I(a,b)=\int _{a}^{b}f(x)\approx \left({\frac {b-a}{2}}\right)\sum _{i=1}^{n}w_{i}f({\frac {b-a}{2}}\xi _{i}+{\frac {b+a}{2}})$

where $w_{i},\xi _{i}$ are the weights and evaluation points, respectively, and $n$ defines the number of integration points. Gauss integration with $n$ integration points can exactly integrate polynomials of order $2n-1$. Consider 3 point Gauss integration with evaluation points $\xi =(-{\sqrt {3/5}},\ 0,\ {\sqrt {3/5}})$ and weights $w=(5/9,8/9,5/9)$, respectively. The derivatives of the integral using Gauss integration can be obtained using CTSE by perturbing $a$ or $b$ as $a+ih$ or $b+ih$. For example, $dI/da$ can be computed as

Eq. 17 $\qquad \qquad \qquad \qquad \qquad \qquad \qquad dI/da={\frac {1}{h}}\Im (I(a+ih,b))$

Eq. 18 $\qquad \qquad \quad \approx {\frac {1}{h}}\Im \left(\left({\frac {b-(a+ih)}{2}}\right)\sum _{i=1}^{n}w_{i}f\left({\frac {b-(a+ih)}{2}}\xi _{i}+{\frac {(b+(a+ih)}{2}}\right)\right)$

and similarly for $dI/db.$

Example 7.1.3:

$\qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad ~~I(a,b)=\int _{a}^{b}2xdx$

Consider the integral $I(a,b)=\int _{a}^{b}2xdx$ with the analytical solutions $I(a,b)=(b^{2}-a^{2})$, $dI/db=2b,\ dI/da=-2a$. In this case, the Gauss integration with 3 integration points is exact. The integral and derivatives were computed using the values $a=1,\ b=2,\ c=2$. The derivatives with respect to $a$ and $b$ computed using CTSE are also exact and agree with the values shown in Table

Example 7.1.4:

$\qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad ~~~~~~~I(a,b,c)=\int _{a}^{b}xe^{-cx}dx$

Consider the integral $I(a,b,c)=\int _{a}^{b}xe^{-cx}dx$, where now the integrand contains a parameter $c$. Derivatives with respect to $a$, and $b$ were computed as shown in equation 17. Derivatives with respect to $c$ were computed using the formula

Eq. 19 $\qquad \qquad \qquad \qquad \qquad \qquad \qquad dI/dc=\Im (I(a,b,c+ih))/h$

Eq. 20 $\qquad \qquad \qquad \qquad \quad {\frac {1}{h}}\Im \left({\frac {b-a}{2}}\sum _{i=1}^{n}w_{i}f\left({\frac {b-a}{2}}\xi _{i}+{\frac {b+a}{2}},c+ih\right)\right)$

$\qquad \qquad \qquad \qquad \qquad \quad ~~{\begin{array}{|c|c|c|c|}\hline &\operatorname {Exact} &\operatorname {Gauss\ Int.} &\operatorname {Relative\ Error} \\\hline \operatorname {Integral} &0.0786069&0.0786094&3.17\times 10^{-5}\\\hline dI/da&-0.135335&-0.135356&1.55\times 10^{-4}\\\hline dI/db&0.0366313&0.0366457&3.93\times 10^{-4}\\\hline dI/dc&-0.109643&-0.109642&-1.27\times 10^{-5}\\\hline \end{array}}$

8 Newton-Raphson Method

The Newton-Raphson Method is an algorithm that can be used to approximate the root of a function. The estimate of the root of function $f$ at iteration $j+1$ can be estimated from the result at step $j$ as

Eq. 21 $\qquad \qquad \qquad \qquad \qquad \qquad \qquad x_{j+1}=x_{j}-{\frac {f(x_{j})}{f'(x_{j})}}$

Rewriting Equation 21 with the help of the CTSE where $x_{j}^{*}=x_{j}+ih$, then $f(x_{i})=\Re (f(x_{j}^{*}))$ and $f'(x_{i})=\Im (f(x_{j}^{*}))/h$ and one obtains

Eq. 22 $\qquad \qquad \qquad \qquad \qquad \qquad x_{j+1}=x_{j}-{\frac {\Re (f(x_{j}+ih))}{\Im (f(x_{j}+ih))/h}}$

Example 8.1:

$\qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad ~~f(x)=\cos(x)-x$

The goal is to find the value of $x$ such that $\cos(x)=x$, or in other words, find the root of the equation $\cos(x)-x=0$; therefore, we define $f(x)=\cos(x)-x$. The Newton-Raphson equation becomes

$\qquad \qquad \qquad \qquad \qquad \qquad ~x_{j+1}=x_{j}-{\frac {\Re (\cos(x+ih)-(x+ih))}{\Im (\cos(x+ih)-(x+ih))/h}}$

This function converges quickly to the numerical value x=0.739085133215161 as shown in the table and figure below starting from an initial value of $x=0.2$.

$\qquad \qquad \qquad \qquad \qquad {\begin{array}{|c|c|c|}\hline \operatorname {Iteration\ no.} &\operatorname {x\ value} &\operatorname {Relative\ Error} \\\hline 0&0.2&\\\hline 1&0.850777122431116&1.51*10^{-1}\\\hline 2&0.741530193469262&3.31*10^{-3}\\\hline 3&0.739086449877212&1.78*10^{-6}\\\hline 4&0.739085133215543&5.18*10^{-13}\\\hline 5&0.739085133215161&0\\\hline \end{array}}$

$\qquad \operatorname {Table\ 4:Convergence\ of\ the\ Newton-Raphson\ algorithm\ for\ the\ function} f(x)=\cos(x)-x.$

Figure 26: Convergence of the Newton-Raphson algorithm for the function f(x)=cos(x)-x starting from a value of x=0.2. The dashed red line indicates the exact solution within 15 decimal places.

9 ODE solver

CTSE can be applied within an ODE solver to compute derivatives of the ODE solution with respect to parameters in the differential equation. The beauty of CTSE is that no alterations to the numerical algorithms are required; merely perturb the parameter of interest using a small step size $h$ along the imaginary axis, then use complex algebra during the computations. The results will be complex with the derivatives times $h$ located in the imaginary part.

In this section, two simple ODE solution methods are employed to demonstrate the approach but CTSE can be applied within more sophisticated ODE solution methods in a similar manner. The examples are provided for first order ODEs but derivatives of parameters for higher order ODEs can be obtained similarly.

Consider the first order ODE

$\qquad \qquad \qquad \qquad \qquad \qquad \qquad ~~y'(x)=f(x,y)=ay(x)+bx-x^{2}$

with the initial condition $y(0)=g$ with the values $a=1,\ b=2,\ g=1$. The solution to the ODE, $y(x)$, will be determined along with the derivatives: $dy(x)/da,\ dy(x)/db,\ dy(x)dg$. The derivates will be determined using CTSE during the solution process to determine $y(x)$.

The analytical solutions are (using $a=1,\ b=2,\ g=1$):

$\qquad \qquad \qquad \qquad \qquad \qquad \quad ~{\begin{aligned}y(x)&=x^{2}+e^{x}\\{\frac {dy(x)}{da}}&=-2-2x-x^{2}+e^{x}(2+x)\\{\frac {dy(x)}{db}}&=-1-x+e^{x}\\{\frac {dy(x)}{dg}}&=e^{x}\end{aligned}}$

The ODE to be solved, using $a=1,\ b=2,\ g=1$, is $y'(x)=y(x)+2x-x^{2}$

9.1 Euler's method

Euler's method uses the knowledge of $y'(x)$ at an initial point $x_{n}$ to project the solution to $x_{n+1}$, that is

$\qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad y(x_{n+1})=h_{ode}f(x_{n},y_{n})$

with $(n,0,1,2,...).\ h_{ode}$ is the step size for the ode solver; this step size is much different than the CTSE step size. In effect, the slope $y'(x,\ y)=f(x,\ y),\ at\ x_{n},\ y_{n}$ is assume to hold for the entire step $h_{ode}$.

To compute derivates of $y$ with respect to parameters within the ODE, we solve the following ODEs using a small step size, e.g., h=10^{-10}, applied to the parameter of interest. In all cases, the derivative is determined as typical using equation 6.

Derivative $dy/da$

To compute $dy/da$, we perturb the parameter $a$ as $a=1+hi$ with $h=10^{-10}$. This results in the ODE

$\qquad \qquad \qquad \qquad \qquad \qquad \qquad \quad y'(x)=(1+ih)y(x)+2x-x^{2}$

with $y(x)=1$.

Table 5 below shows the first 4 steps taken using Euler's method using an ODE step size of $h_{ode}=0.1$, and $h=10^{-10}$ for CTSE. The derivative, $dy/da$, can be obtained as always using CTSE by taking the imaginary part and dividing by $h$. The results indicate that the relative error of $dy/da$ is slightly larger than for $y$ and that both grow larger with increasing $x$.

${\begin{array}{|c|c|c|c|c|c|}\hline x&\operatorname {Euler's\ method\ output} &y&y_{exact}&\operatorname {RelativeError} &dy/da&dy/da_{exact}&\operatorname {RelativeError} \\\hline 0.0&-&1.0&1.0&-&-&0.0&-\\\hline 0.1&1.1+1.*10^{-11}i&1.1&1.1152&1.36*10^{-2}&0.1&0.1108&9.80*10^{-2}\\\hline 0.2&1.229+2.2*10^{-11}i&1.229&1.2614&2.57*10^{-2}&0.22&0.2471&1.09*10^{-1}\\\hline 0.3&1.3879+3.649*10^{-11}i&1.388&1.4399&3.61*10^{-2}&0.3649&0.4147&1.20*10^{-1}\\\hline 0.4&1.578+5.402*10^{-11}i&1.578&1.6518&4.49*10^{-2}&0.5402&0.6204&1.29*10^{-1}\\\hline \end{array}}$

$\quad ~{\begin{aligned}\operatorname {Table:5\ Results\ using\ Euler's\ method\ to\ } &\operatorname {solve\ for\ y(x)\ and\ dy(x)/da\ using\ an\ ode\ stepsize\ of\ } \\h_{ode}=0.1\operatorname {and\ a} &\operatorname {\ CTSE\ step\ size\ of} h=10^{-10}\end{aligned}}$

Derivative $dy/db$

To compute $dy/db$, we perturb the parameter $b$ as $b=2+hi$ with $h=10^{-10}$. This results in the ODE

$\qquad \qquad \qquad \qquad \qquad \qquad \qquad \quad y'(x)=y(x)+(2+ih)x-x^{2}$

with $y(x)=1$.

Table 6 below shows the first 4 steps taken using Euler's method using an ODE step size of $h_{ode}=0.1$, and $h=10^{-10}$ for CTSE. The derivative, $dy/db$, can be obtained as always using CTSE by taking the imaginary part and dividing by $h$. The results indicate that the relative error of $dy/db$ is slightly larger than for $y$ and that both grow larger with increasing $x$. The results for $y(x)$ are identical to those in Table 5.

$\qquad \qquad \qquad \quad {\begin{array}{|c|c|c|c|c|}\hline x&\operatorname {Euler's\ method\ output} &dy/db&dy/db_{exact}&\operatorname {RelativeError} \\\hline 0.0&-&-&0.0&-\\\hline 0.1&1.1+0i&0.0&0.0052&1.0\\\hline 0.2&1.229+1.0*10^{-12}i&0.01&0.0214&5.33*10^{-1}\\\hline 0.3&1.3879+3.1*10^{-12}i&0.031&0.0499&3.78*10^{-1}\\\hline 0.4&1.578+6.41*10^{-12}i&0.0641&0.0918&3.02*10^{-1}\\\hline \end{array}}$

$\qquad {\begin{aligned}\operatorname {Table\ 6:Results\ using\ Euler's\ method\ to\ solve\ for} &\ dy(x)/db\operatorname {using\ an\ ode\ size\ of} h_{ode}=0.1\\\operatorname {and\ a\ CTSE\ step\ } &\operatorname {size\ of} h=10^{-10}\end{aligned}}$

Derivative $dy/dg$

To compute $dy/dg$, the differential equation is unchanged but we perturb the parameter of the initial condition $g$ as $g=1+hi$ with $h=10^{-10}$. This results in the ODE

$\qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad ~~~~y'(x)=y(x)+2x-x^{2}$

with $y(x)=1+hi$.

Table 7 below shows the first 4 steps taken using Euler's method using an ODE step size of $h_{ode}=0.1$, and $h=10^{-10}$ for CTSE. The derivative, $dy/dg$, can be obtained as always using CTSE by taking the imaginary part and dividing by $h$. The results indicate that the relative error of $dy/dg$ is slightly larger than for $y$ and that both grow larger with increasing $x$. The results for $y(x)$ are identical to those in Table 5.

$\qquad \qquad \qquad \quad {\begin{array}{|c|c|c|c|c|}\hline x&\operatorname {Euler's\ method\ output} &dy/dg&dy/dg_{exact}&\operatorname {RelativeError} \\\hline 0.0&-&-&0.0&-\\\hline 0.1&1.1+1.1*10^{-10}i&1.1&1.011&4.68*10^{-3}\\\hline 0.2&1.229+1.1*10^{-10}i&1.21&1.221&9.34*10^{-3}\\\hline 0.3&1.388+1.1*10^{-10}i&1.331&1.350&1.40*10^{-2}\\\hline 0.4&1.578+1.1*10^{-10}i&1.464&1.492&1.86*10^{-2}\\\hline \end{array}}$

$\qquad {\begin{aligned}\operatorname {Table\ 7:Results\ using\ Euler's\ method\ to\ solve\ for} &\ dy(x)/dg\operatorname {using\ an\ ode\ size\ of} h_{ode}=0.1\\\operatorname {and\ a\ CTSE\ step\ } &\operatorname {size\ of} h=10^{-10}\end{aligned}}$

Results for $y(x)$, $dy(x)/da$, $dy(x)/db$, and $dy(x)/dg$ from solving the ODE $y'(x)=ay(x)+bx-x^{2}$ using Euler's method with step size $0.1$ and CTSE step size $10^{-10}$ (red dots) compared with the analytical solutions (solid line).

10 Cauchy-Riemann (CR) representation of a complex number

CTSE requires support for complex variables. Although all modern languages have this support, it may be needed in some cases, or more convenient, to augment an existing code using the Cauchy-Riemann (CR) representation of a complex number rather than using complex variables. The CR matrix form of a complex number is

$\qquad \qquad \qquad \quad \qquad \qquad \qquad \quad \qquad \quad x^{*}=A^{*}=a+bi=\left[{\begin{array}{cc}a&-b\\b&a\end{array}}\right]$

where $a$ and $b$ are real numbers.

As shown, the real value is always located along the diagonal with the imaginary number located in the off-diagonal locations differing by a minus sign. Therefore, in a CR formulation of CTSE, the real number can be retrieved from the appropriate location of the CR matrix, e.g., $a=A_{11}^{*}or\ a=A_{22}^{*}$, and the imaginary number retrieved as $b=A_{21}^{*}or\ b=-A_{12}^{*}$.

The important aspect of the CR representation is that the matrix contains all real numbers. In effect, mathematical operations on complex numbers can be replaced by matrix operations of all real matrices. For example, consider$(a+bi)(c+di)=(ac-bd)+(ad+bc)i=\left[{\begin{array}{cc}(ac-bd)&-(ad+bc)\\(ad+bc)&(ac-bd)\end{array}}\right]$. This same result can be obtained from the matrix multiplication of the equivalent CR matrices, e.g., $\left[{\begin{array}{cc}a&-b\\b&a\end{array}}\right]\left[{\begin{array}{cc}c&-d\\d&c\end{array}}\right]=\left[{\begin{array}{cc}(ac-bd)&-(ad+bc)\\(ad+bc)&(ac-bd)\end{array}}\right]$. Therefore, we can implement CTSE without complex numbers by using the CR form of the complex number and matrix operations. In fact, all the closed-form examples listed above can be solved using CR matrices. One important point though, is that matrix functions must be used for intrinsic functions such as the $\exp ,\ \sin ,\ \cos ,\ \log$, and square-root functions.

Example 10.1:

$\qquad \qquad \qquad \quad \qquad \qquad \qquad \quad \qquad \qquad \qquad f(x)=x^{2}$

Consider the simple function $f(x)=x^{2}$. Then

$\qquad \qquad \qquad \qquad ~{\begin{aligned}f(x+ih)=A^{*}=\left[{\begin{array}{cc}x&-h\\h&x\end{array}}\right]\left[{\begin{array}{cc}x&-h\\h&x\end{array}}\right]&=\left[{\begin{array}{cc}x^{2}-h^{2}&-2hx\\2hx&x^{2}-h^{2}\end{array}}\right]\\{\frac {df}{dx}}=A_{21}^{*}={\frac {2h\,x}{h}}=&2x\end{aligned}}$

We see in that for this example CTSE is exact as there are no higher order terms with respect to $h$.

Example 10.2:

$\qquad \qquad \qquad \quad \qquad \qquad \qquad \quad \qquad \qquad \quad ~~f(x)=\sin(x)$

$\qquad \qquad \quad {\begin{aligned}f(x+ih)=A^{*}=\sin(\left[{\begin{array}{cc}x&-h\\h&x\end{array}}\right]=\left[{\begin{array}{cc}\cosh(h)\sin(x)&-\cos(x)\sinh(h)\\\cos(x)\sinh(h)&\cosh(h)\sin(x)\end{array}}\right],&\operatorname {then} \\{\frac {df}{dx}}=A_{21}^{*}={\frac {\cos(x)\sinh(x)}{h}}&=cos(x)\end{aligned}}$

as $h\rightarrow 0$. Note, that we must use the matrix function for sine.

Example: Finite element example using a beam element

The finite element example presented in Section 6 is repeated here and solved using the Cauchy-Riemann form. Along the way, the conversion from complex to CR form is demonstrated. This example is instructive because it shows how the complex-variable system of equations can be solved using real numbers. In other words, if a complex-variable solver is not available, then the CR form can be employed in conjunction with a real-variable solver.

Formulation of the Cauchy-Riemann system of equations

The complex-variable system of equations, $K^{*}u^{*}=f^{*}$ is reformulated to consist of only real numbers. To convert the system of equations to the CR form, each complex variable must be expanded to a $2\times 2$ matrix of all real numbers.

For example, $k_{11}^{*}=\left[{\begin{array}{cc}Re[k_{11}^{*}]&-Im[k_{11}^{*}]\\Im[k_{11}^{*}]&Re[k_{11}^{*}]\end{array}}\right]$, etc. The right hand side must be expanded for compatibility, that is, $\left\{{\begin{array}{c}P\\0\end{array}}\right\}=\left[{\begin{array}{cc}Re[P]&-Im[P]\\Im[P]&Re[P]\\0&0\\0&0\end{array}}\right]$. The complete system of equations for this example becomes

$\left[{\begin{array}{cccc}\Re [k_{11}^{*}]&-\Im [k_{11}^{*}]&\Re [k_{12}^{*}]&-\Im [k_{12}^{*}]\\\Im [k_{11}^{*}]&\Re [k_{11}^{*}]&\Im [k_{12}^{*}]&\Re [k_{12}^{*}]\\\Re [k_{21}^{*}]&-\Im [k_{21}^{*}]&\Re [k_{22}^{*}]&-\Im [k_{22}^{*}]\\\Im [k_{21}^{*}]&\Re [k_{21}^{*}]&\Im [k_{22}^{*}]&\Re [k_{22}^{*}]\end{array}}\right]\left\{{\begin{array}{cc}\delta &-d\delta /dL\cdot h\\d\delta /dL\cdot h&\delta \\\theta &-d\theta /dL\cdot h\\d\theta /dL\cdot h&\theta \end{array}}\right\}=\left\{{\begin{array}{cc}\Re [P]&-\Im [P]\\\Im [P]&\Re [P]\\0&0\\0&0\end{array}}\right\}$ .

The deflection and rotation can be extracted from rows 1 and 3 from column 1, and their derivatives can be obtained from rows 2 and 4 of column 1 after dividing by h.

Numerical Example

A numerical example will be presented for clarity using the numerical values: $L=0.5m,\ E=20010^{9}Pa,\ I=1/12*(0.02m)^{4}$, i.e., rectangular cross section of 2 cm in each direction, load $P=10N$, and $h=10^{-10}$. As a result, the stiffness matrix and right hand side become

$k=\left[{\begin{array}{cc}256,000&64,000\\64,000&256,000\end{array}}\right]$, and $\left\{{\begin{array}{c}10\\0\end{array}}\right\}$ . Solving this system yields $\left\{{\begin{array}{c}\delta \\\theta \end{array}}\right\}=\left\{{\begin{array}{c}0.00015625\\-0.000468\end{array}}\right\}$.

Sensitivity with respect to load P

Consider the case where the load $P$ is perturbed as $P+ih$. The stiffness matrix is real only (no imaginary terms required) as

$k^{*}=\left[{\begin{array}{cc}256,000&64,000\\64,000&256,000\end{array}}\right]$, whereas the right hand side becomes $\left\{{\begin{array}{c}P\\0\end{array}}\right\}=\left\{{\begin{array}{c}10+10^{-10}i\\0\end{array}}\right\}$Solving for the complex-variable degrees-of-freedom as $\left\{{\begin{array}{c}\delta ^{*}\\\theta ^{*}\end{array}}\right\}=K^{-1}\left\{{\begin{array}{c}P\\0\end{array}}\right\}$ yields

$\qquad \qquad \qquad \qquad \qquad \qquad ~\left\{{\begin{array}{c}\delta ^{*}\\\theta ^{*}\end{array}}\right\}=\left\{{\begin{array}{c}0.00015625-1.5625*10^{-15}i\\-0.00046875-4.6875*10^{-15}i\end{array}}\right\}$.

Extracting the real and imaginary parts and dividing the imaginary parts by $h$ yields the results $\delta =0.00015625,\ d\delta /dP=0.000015625,\ \theta =-0.00046875,\ d\theta /dP=-0.00046875$.

The CR system of equations becomes

$\left[{\begin{array}{cccc}256,000&0&64,000&0\\0&256,000&0&64,000\\64,000&0&21,333.3&0\\0&64,000&0&21,333.3\end{array}}\right]\left\{{\begin{array}{cc}\delta &-d\delta /dL\cdot h\\d\delta /dL\cdot h&\delta \\\theta &-d\theta /dL\cdot h\\d\theta /dL\cdot h&\theta \end{array}}\right\}=\left\{{\begin{array}{cc}10&-10^{-10}\\10^{-10}&10\\0&0\\0&0\end{array}}\right\}$ .

Solving the system of equations yields $\left\{{\begin{array}{cc}\delta &-d\delta /dL\cdot h\\d\delta /dL\cdot h&\delta \\\theta &-d\theta /dL\cdot h\\d\theta /dL\cdot h&\theta \end{array}}\right\}=\left\{{\begin{array}{cc}0.00015625&-1.5625*10^{-15}\\1.5625*10^{-15}&0.00015625\\-0.00046875&4.6875*10^{-15}\\-4.6875*10^{-15}&-0.00046875\end{array}}\right\}$.

From these results we obtain values that agree with the complex-variable solution, namely $\delta =0.00015625,\ d\delta /dP=0.000015625,\ \theta =-0.00046875,\ d\theta /dP=-0.00046875$.

Sensitivity with respect to elastic modulus E

Consider the case where $E$ is perturbed as $E+ih$. Then the complex stiffness matrix becomes

$\qquad \qquad \qquad \qquad \qquad k^{*}=\left[{\begin{array}{cc}256,000+9.6*10^{-9}i&64,000+2.4*10^{-9}i\\64,000+2.4*10^{-9}i&256,000+8.0*10^{-9}i\end{array}}\right]$.

Solving for the complex-variable degrees-of-freedom as ${\left\{{\begin{array}{c}\delta ^{*}\\\theta ^{*}\end{array}}\right\}}={K^{*^{-1}}\left\{{\begin{array}{c}P\\0\end{array}}\right\}}$ yields

$\qquad \qquad \qquad \qquad \qquad \quad \left\{{\begin{array}{c}\delta ^{*}\\\theta ^{*}\end{array}}\right\}=\left\{{\begin{array}{c}0.00015625-5.85938*10^{-18}i\\-0.00046875+1.75781*10^{-17}i\end{array}}\right\}$

Extracting the real and imaginary parts and dividing the imaginary parts by $h$ yields the results $\qquad \qquad \delta =0.00015625,\ d\delta /dE=-5.85938*10^{-8},\ \theta =-0.00046875,\ d\theta /dE=1.75781*10^{-7}$.

The CR system of equations becomes

$\left[{\begin{array}{cccc}256000&-9.6*10^{-9}&64000&-2.4*10^{-9}\\9.6*10^{-9}&256000&2.4*10^{-9}&64000\\64000&-2.4*10^{-9}&256000&-8.0*10^{-9}\\2.4*10^{-9}&64000&8.0*10^{-9}&256000\end{array}}\right]\left\{{\begin{array}{cc}\delta &-d\delta /dL\cdot h\\d\delta /dL\cdot h&\delta \\\theta &-d\theta /dL\cdot h\\d\theta /dL\cdot h&\theta \end{array}}\right\}=\left\{{\begin{array}{cc}10&0\\0&10\\0&0\\0&0\end{array}}\right\}$ .

Solving the system of equations yields $\qquad \qquad \qquad \left\{{\begin{array}{cc}\delta &-d\delta /dL\cdot h\\d\delta /dL\cdot h&\delta \\\theta &-d\theta /dL\cdot h\\d\theta /dL\cdot h&\theta \end{array}}\right\}=\left\{{\begin{array}{cc}0.00015625&5.85938*10^{-18}\\-5.85938*10^{-14}&0.00015625\\-0.00046875&-1.75781*10^{-17}\\1.75781*10^{-17}&-0.00046875\end{array}}\right\}$.

From these results we obtain values that agree with the complex-variable solution, namely $\qquad \qquad \delta =0.00015625,\ d\delta /dE=-5.85938*10^{-8},\ \theta =-0.00046875,\ d\theta /dE=1.75781*10^{-7}$.

Sensitivity with respect to length L

Consider the case where the length $L$ is perturbed as $L+ih$. Then the complex stiffness matrix becomes

$\qquad \qquad \qquad \qquad \qquad k^{*}=\left[{\begin{array}{cc}256,000-0.0001536i&64,000-0.0000256i\\64,000-0.0000256i&256,000-0.0001536i\end{array}}\right].$

Solving for the complex-variable degrees-of-freedom as ${\left\{{\begin{array}{c}\delta ^{*}\\\theta ^{*}\end{array}}\right\}}={K^{*^{-1}}\left\{{\begin{array}{c}P\\0\end{array}}\right\}}$ yields

$\qquad \qquad \qquad \qquad \qquad \qquad \qquad \left\{{\begin{array}{c}\delta ^{*}\\\theta ^{*}\end{array}}\right\}=\left\{{\begin{array}{c}0.00015625+9.375*10^{-14}i\\-0.00046875-1.875*10^{-13}i\end{array}}\right\}$.

Extracting the real and imaginary parts and dividing the imaginary parts by $h$ yields the results $\qquad \qquad \qquad \delta =0.00015625,\ d\delta /dL=0.00009375,\ \theta =-0.00046875,\ d\theta /dL=-0.0001875$.

The CR system of equations becomes

$\left[{\begin{array}{cccc}256,000&0.0001536&64,000&0.0000256\\-0.0001536&256,000&-0.0000256&64,000\\64,000&0.0000256&21,333.3&4.267*10^{-6}\\-0.0000256&64,000&-4.267*10^{-6}&21,333.3\end{array}}\right]\left\{{\begin{array}{cc}\delta &-d\delta /dL\cdot h\\d\delta /dL\cdot h&\delta \\\theta &-d\theta /dL\cdot h\\d\theta /dL\cdot h&\theta \end{array}}\right\}=\left\{{\begin{array}{cc}10&0\\0&10\\0&0\\0&0\end{array}}\right\}$ . Note that the stiffness matrix is no longer symmetric.

Solving the system of equations yields $\left\{{\begin{array}{cc}\delta &-d\delta /dL\cdot h\\d\delta /dL\cdot h&\delta \\\theta &-d\theta /dL\cdot h\\d\theta /dL\cdot h&\theta \end{array}}\right\}=\left\{{\begin{array}{cc}0.00015625&-9.375*10^{-14}\\9.375*10^{-14}&0.00015625\\-0.00046875&1.875*10^{-13}\\-1.875*10^{-13}&-0.00046875\end{array}}\right\}$

From these results we obtain values that agree with the complex-variable solution, namely $\qquad \qquad \qquad \delta =0.00015625,\ d\delta /dL=0.00009375,\ \theta =-0.00046875,\ d\theta /dL=-0.0001875$.

11 Overloading functions to provide CTSE support

There are some functions for which CTSE will not work. For example, a complex-variable version may not exist, or the complex-variable version may return a only a real result, or the complex-variable implementation may not provide machine precision derivatives. For example, the inverse error function, $InverseErf$, does not have a complex-variable version. Another situation is where a functions has a complex-variable form but the result is real only. Such a case is the Absolute function, that is $Abs(a+bi)={\sqrt {a^{2}+b^{2}}}$. Third, in some implementations of the inverse trigonometric functions, the relative error can be reduced to $10^{-11}$ but no further. However, as shown below, an artificial “pseudo-complex” version can be created such that CTSE can be used in its standard form for all these cases. In other words, when the CTSE equation 6 is applied, the pseudo-complex version is called which returns the correct result in the imaginary part.

The key is to develop a first order complex-variable Taylor series of the form

Eq. 23 $\qquad \qquad \qquad \qquad \qquad \qquad \qquad f(x^{*})=f(x)+f^{'}(x)(x^{*}-x_{0})$

for the function of interest. Here $x^{*}=x+ih$, $f'(x)$ represents the derivative of the function evaluated at $x$, as a result$(x^{*}-x)i=hi$ and equation $x$ can be written

Eq. 24 $\qquad \qquad \qquad \qquad \qquad \qquad \qquad \quad f(x^{*})=f(x)+f^{'}(x)hi$

Thus, the only requirement is to define the imaginary portion of the function as $f'(x)h$ so that when equation 6 is applied $\qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad ~~~~~~{\frac {{\text{Im}}\left[f(x+ih)\right]}{h}}=f'(x)$

11.1 Example:

$\cos(x)$ function

To fix ideas, let's compute an overloaded function for cosine that implements equation 24. Clearly we don't need to do this since there exists complex versions of the cosine function but it instructive to do so. Applying equation 24, we have

Eq. 25$\qquad \qquad \qquad \qquad \qquad \qquad \quad \cos(x^{*})=\cos(x)-\sin(x)hi$

Hence, when this version of cosine is called with a complex argument, the imaginary result will be $-\sin(x)h$. Therefore equation 6 will return -$\sin(x)$ as the derivate. It should be clear that this version of the complex-variable cosine function is not equivalent to the standard cosine function. It it only useful for computing derivates using CTSE.

11.2 Example:

$Abs$ function

Consider the absolute operator $Abs(x)$. CTSE does not work on this operator since the $Abs$ operator applied to a complex number yields a real result, i.e., $Abs(a+bi)={\sqrt {a^{2}+b^{2}}}$. However, we can derive a function that uses the real function as $\Re (Abs(x))={\sqrt {x^{2}}}$ and its derivative as the imaginary part, $\Im (hf^{(1)}(x_{0}))={\frac {dAbs(x)}{dx}}h={\frac {x}{\sqrt {(x^{2})}}}h$. Therefore, we can define the overloaded $Abs$ for a complex argument of $x^{*}=x_{0}+hi$ as

Eq. 26 $\qquad \qquad \qquad \qquad \qquad \qquad \qquad Abs(x^{*})={\sqrt {x^{2}}}+{\frac {x}{\sqrt {x^{2}}}}hi$

Then if we evaluate, for example, $\Im (Abs(-1+hi)/h=-1$ for any $h$ and $\Im (Abs(1+hi)/h=1$ for any $h$. A plot of the overloaded $Abs$ function is shown in Figure 26.

Figure 27: Plot of the Absolute function and its derivative

11.3 Example:

$InverseErf$ function

As another example, consider the inverse error function, $InverseErf(x)$. In the current programming languages$InverseErf$ only accepts a real variable. We need will overload this definition to support complex arguments.

This function accepts a real argument between $-1<x<+1$.

The derivative is $dInverseErf(x)/dx={\frac {\sqrt {\pi }}{2}}e^{InverseErf(x)^{2}}$. Hence, the overloaded function for a complex argument $x^{*}=x+hi$ is

Eq. 27 $\qquad \qquad \qquad \quad InverseErf(x+hi)=InverseErf(x)+{\frac {\sqrt {\pi }}{2}}e^{InverseErf(x)^{2}}hi$

A plot of the this $InverseErf$ function and its derivative is shown in Figure 27.

Figure 28: Plot of the InverseErf function and its derivative

12 Miscellaneous topics

A number of miscellaneous topics are discussed below that may be of interest to the reader.

12.1 Accuracy of the numerical derivatives computed by CTSE

The CTSE method itself can obtain machine precision accuracy, typically 15 decimal places, on simple functions; however, the accuracy of the derivative depends upon the numerical algorithm to which it is being applied and the output that is being computed. That is, the accuracy of the derivative will only be as good as the accuracy of the numerical method itself. This result was clearly shown in the application of CTSE within numerical integration using Simpson's rule, see Table Experience to date indicates that the derivatives computed using CTSE are usually of the same or one order less accurate than the real-variable results. For example, if the relative error of the real-variable result was of the order $O(10^{-4})$, the relative error of the derivative will be of the order $O(10^{-4}-10^{-3})$. Results to date are entirely empirical, no formal mathematical investigation has been conducted.

12.2 Computing mixed and higher order derivatives using CTSE combined with the finite difference method

A natural question is whether CTSE can be used to calculate higher order derivatives such as second order. The short answer is no, CTSE cannot be used to calculate higher order derivates without subtractive error. For example, the application of CTSE to compute derivatives of second order, say by taking 2 steps up the imaginary axis, leads to a subtraction of terms, and, hence, the method becomes step size dependent and it loses significant accuracy. To compute higher order derivatives without subtractive error requires multicomplex or multidual variables and algebra. These topics are covered in subsequent chapters.

However, one can combine CTSE for first order derivatives with the classical finite difference estimation for first order derivatives to obtain second order derivatives that are numerically superior to using the classical finite difference method for second order derivatives. In a nutshell, one can apply a first order finite difference formula to the derivatives estimated by CTSE to approximate a second order derivative, i.e.,

Eq. 28 $\qquad \qquad \qquad \qquad \qquad \qquad \qquad f''\approx {\frac {f'(x+h_{Re})-f'(x)}{h_{Re}}}$

Contrast this equation with the classical second order finite difference formula

Eq. 29 $\qquad \qquad \qquad \qquad \qquad ~~f''(x)\approx {\frac {f(x+2h_{Re})-2f(x+h_{Re})+f(x)}{h_{Re}^{2}}}$

Using CTSE to compute $f'$, yields the equation

Eq. 30 $\qquad \qquad \qquad \qquad \quad f''(x)\approx {\frac {\Im (f((x+ih)+h_{\Re })h-\Im (f(x+ih))/h}{h_{\Re }}}$

Eq. 31 $\qquad \qquad \qquad \qquad \qquad \quad ~~={\frac {\Im (f(x+ih+h_{\Re })-\Im (f(x+ih))}{hh_{\Re }}}$

In contrast to Figure 1 where the perturbation is conducted along the real OR imaginary axes, this method perturbs the variable of interest along BOTH axes simultaneously. As shown in figure 28, one perturbs the variable $x$ by $h_{\Re }$ along the real axis and $h$ along the imaginary axis. However, as we know, $h<<h_{\Re }$. As an example, consider the second order derivative estimate of the cosine function computed using equation 31,

$\qquad \qquad \qquad \qquad \qquad \qquad f''(x)\approx {\frac {\Im (\cos(x+h_{\Re }+ih)-\Im (\cos(x+ih))}{hh_{\Re }}}$

Figure 29: Strategy for computing second order derivatives using CTSE and forward differencing

Figure 29 below shows the relative error of the second order derivative of $\cos(x)$ evaluated at $x=\pi /6$ as a function of the real-variable finite difference step size $h_{\Re }$. For all runs, the CTSE step size was set equal to $h=10^{-10}$ for all analyses so that the first order derivative computed using CTSE was accurate to machine precision. Hence, all inaccuracy is due to the truncation and subtractive error of the finite difference formulas. As shown in the figure, three additional orders of magnitude reduction in the relative error can be obtained using a combination of CTSE for the first order derivative combined with forward differencing compared to the classical second order forward difference formula. Also, it should be noted that for value of $h_{\Re }$ between $10^{-9}$ and $10^{-15}$ the 2nd order finite difference produced grossly inaccurate estimates of the 2nd order derivatives. These values are not shown in the figure.

Figure 30: Comparison of second order derivative of $\cos(x)$ at $x=\pi /6$. Second order finite difference (blue triangles) vs. combined CTSE and first order finite difference (red circles); $h=10^{-10}$ for all CTSE analyses.

12.3 Can CTSE be applied to a functions that are already complex

Can one apply CTSE to a function that is already complex? In short, no. One cannot apply equation 6 to a function that is already complex. One must use a different imaginary axis other than $i$ (since that axis is being used in the complex function) to apply the perturbations and extract the derivative. One can, for example, use dual numbers that have an imaginary number $\epsilon$  thereby combining complex and dual numbers.

13 The use of Quaternions to compute first order derivates

Readers familiar with quaternions may be interested to know that they can be used to compute first order derivatives without subtractive error in the same manner as using CTSE.

Quaternions are constructed with one real and 3 imaginary axes $i,\ j,\ k$, where $i*j=k$, etc. The algebraic rules for quaternions are well established.

Quaternion numbers are a subset of hypercomplex numbers con-sting of one real and 3 imaginary axes $i,\ j,\ k$. A number $q^{*}$ is defined as

$\qquad \qquad \qquad \qquad \qquad \qquad \qquad ~q^{*}=q_{_{\Re }}+q_{i}i+q_{j}j+q_{k}k$

where $q_{_{\Re }},\ q_{i},\ q_{j},\operatorname {and} q_{k}$ are real coefficients and are imaginary units governed by the following conditions

$\qquad \qquad \qquad \qquad \qquad \qquad \qquad ~~i^{2}=j^{2}=k^{2}=ijk=-1$

$\qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \quad ij=-ji=k$

$\qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad ~~~jk=-kj=i$

$\qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad ~~~ki=-ik=j$

The value for using quaternions is that 3 first order derivates can be computed with a single analysis, whereas 3 CTSE analyses would be required. Of course, a quaternions library is needed since quaternions are typically not intrinsic to most programming languages. As such, the runtime of 3 CTSE runs may in fact be faster than a single quaternion analysis. This aspect will be implementation dependent.

The Quaternion Taylor Series Expansion (QTSE) method uses quaternion algebra to compute up to 3 first order derivatives of a real function evaluated at $x=x_{0},\ y=y_{0},\ z=z_{0}$. In this approach, the variables $x$, $y$ and $z$ are perturbed by a small step size $h$ in imaginary directions $i,\ j$ and $k$ respectively:

$\qquad \qquad \qquad \qquad \qquad \qquad \qquad {\begin{aligned}x^{*}=x_{0}+h{\boldsymbol {i}}+0j+0k,\\y^{*}=y_{0}+0i+h{\boldsymbol {j}}+0k,\\z^{*}=z_{0}+0i+0j+h{\boldsymbol {k}}{\phantom {,}}\end{aligned}}$

where $h<<1$. In a similar manner to CTSE, it is recommended that $h$ be $10^{-10}$ or smaller.

The function $f$ is now evaluated as a quaternion to yield ${\mathcal {F}}(x^{*},y^{*},z^{*})$. The real-variable result and the derivatives can be extracted from the imaginary axes of the quaternion evaluation of the function $f$ as

$\qquad \qquad \qquad \qquad \qquad \qquad \quad ~(x_{0},y_{0},z_{0})\approx \Re \left[{\mathcal {F}}(x^{*},y^{*},z^{*})\right]$

$\qquad \qquad \qquad \qquad \qquad \qquad {f_{,x}(x_{0},y_{0},z_{0})}\approx {\Im _{i}\left[{\mathcal {F}}(x^{*},y^{*},z^{*})\right]}/h$

$\qquad \qquad \qquad \qquad \qquad \qquad {f_{,y}(x_{0},y_{0},z_{0})}\approx {\Im _{j}\left[{\mathcal {F}}(x^{*},y^{*},z^{*})\right]}/h$

$\qquad \qquad \qquad \qquad \qquad \qquad {f_{,z}(x_{0},y_{0},z_{0})}\approx {\Im _{k}\left[{\mathcal {F}}(x^{*},y^{*},z^{*})\right]}/h$

where $\Re \left[\cdot \right]$ extracts the real coefficient of the quaternion number, $\Im _{i}\left[\cdot \right]$ extracts the coefficient associated with the imaginary direction $i$, etc. Since no subtractions are performed during quaternion algebra operations, the method does not suffer from subtraction cancellation error. In practical numerical applications, the minimum possible error is usually bounded by machine error.

Example: $f=\sin(x^{3}*y^{2}*z^{4})$ where $x,\ y,\ z$ are quaternions

Consider the sine function $f=\sin(x^{3}\ y^{2}\ z^{4})$ where $x,\ y,$ and $z$ are quaternions. In this case the quaternion Taylor series method is used to compute the derivatives with respect to$x,\ y,$ and $z$ evaluated at $x=0.5,\ y=1.5,\ z=2.5$. The $x,\ y,$ and $z$ quaternions are defined as shown in equation 37 as

$\qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \quad {\begin{array}{c}x^{*}=0.5+{\boldsymbol {hi}},\\y^{*}=1.5+{\boldsymbol {hj}},\\z^{*}=2.5+{\boldsymbol {hk}}{\phantom {,}}\end{array}}$

where $h=10^{-15}$.

Using the previous equations, the derivatives can be obtained to machine precision using a step size of $h=10^{-15}$ or smaller.

$\qquad \qquad \qquad \qquad \qquad \qquad \qquad {\begin{array}{c}df/dx=Im_{i}(\sin(x^{*^{3}}*y^{*^{2}}*z^{*^{4}}))/h,\\df/dy=Im_{j}(\sin(x^{*^{3}}*y^{*^{2}}*z^{*^{4}}))/h,\\df/dz=Im_{k}(\sin(x^{*^{3}}*y^{*^{2}}*z^{*^{4}}))/h{\phantom {,}}\end{array}}$

where $\Im _{i}$ means to extract the real coefficient of the $i^{th}$ direction, and similarly for $\Im _{j}$ and $\Im _{k}$.

The results are $df/dx=-0.609479570660066$, $df/dy=-0.135439904591125$, and $df/dz=-0.162527885509352$.

These results are accurate to machine precision, extensions to quaternions can also be used to determine multiple first order derivatives. For example, octonions (one real and 7 imaginary axes) can be used to compute seven first order derivatives in a single analysis; sedenions (one real and 15 imaginary axes) can be used to compute seven first order derivatives in a single analysis, etc.

14 Comments on the accuracy of the CTSE method

In many publications, and repeated here, it is stated that CTSE will provide machine precision derivatives if sufficiently small h is used. This is true if the function to which it is being applied provides results that are also machine precision accurate. In general, the accuracy of the derivatives from CTSE will be no more accurate than the algorithm in which it is being employed. In the application to Simpson’s rule, the accuracy is only ∼$10^{-3}$, see Figure 25. This is because the accuracy of Simpson’s rule for the real part is only ∼$10^{-3}$. This behavior is true in general. A general rule of thumb is that the derivatives are slightly less accurate than the real part. However, it is true that CTSE (with small enough h) will provide the most accurate answer possible for the particular application, that is, the accuracy will be the same as if the analytical derivative was taken first, then the numerical algorithm applied.

15 Background reading

There are several excellent references that discuss CTSE and provide examples. As we have indicated, CTSE is not limited to any particular mathematical or science domain. It has been applied successfully in a myriad of disciplines. The article by Squire and Trapp, 1998, provides a concise summary of CTSE. Martins, et al. describe a scripting procedure that can be used to automatically modify Fortran and C/C++ codes to implement CTSE. There have been multiple applications of CTSE within various domains such as fluid mechanics, eigensensitivity analysis, heat transfer, fracture mechanics, 35 nonlinear structural mechanics, among others. More information regarding the use of quaternions for obtaining derivatives can be obtained from the references by Turner and Cano.

Further reading

• W. Squire, G. Trapp, "Using complex variables to estimate derivatives of real functions," SIAM Review 40 (1998) 110–112
• J.J. Martins, "The complex-step derivative approximation," ACM Transactions on Mathematical Software 29 (2003) 245–263
• W.K. Anderson, J.C. Newman, D.L. Whitfield, E.J. Nielsen, "Sensitivity analysis for Navier-Stokes equations on unstructured meshes using complex variables," AIAA Journal 39 (2001) 56–63.
• B.P. Wang, A.P. Apte, "Complex variable method for eigensolution sensitivity analysis," AIAA Journal. 44 (2006) 2958–2961.
• X.W. Gao, M.C. He, "A new inverse analysis for multi-region heat conduction BEM using complex-variable differentiation method," Engineering Analysis with Boundary Elements 29 (2005) 788–795.
• H.R. Millwater, D. Wagner*, A. Baines*, and A. Montoya, “A Virtual Crack Extension Method to Compute Energy Release Rates using a Complex Variable Finite Element Method,” Engineering Fracture Mechanics 162 (2016) 95–111, https://dx.doi.org/10.1016/j.engfracmech.2016.04.002
• A. Montoya, R. Fielder*, A. Gomez-Farias*, H. Millwater, “Finite Element Sensitivity for Plasticity using Complex Variable Methods,” J. Eng. Mech. 141 2 (2014), https://ascelibrary.org/doi/abs/10.1061/%28ASCE%29EM.1943-7889.0000837