TPS Sizing with Complex Step Method

TPS Sizing Optimization Using Complex Variable Differentiation Sensitivity

I stumbled upon an interesting old presentation that shows a neat application of the complex step method of calculating numerical derivatives for use in optimizing thermal protection system (TPS) thickness. The great thing about the method is that it is minimally intrusive.

Uncertainty Quantification with Stochastic Collocation

I’ve been posting lots of links to uncertainty quantification (UQ) references lately in comments (eg here, here, here), but I don’t really have a post dedicated only to UQ. So I figured I’d remedy that situaton. This post will be a simple example problem based on methods presented in some useful recent papers on UQ (along with a little twist of my own):

• A good overview of several of the related UQ methods [1]
• Discussion of some of the sampling approaches for high dimensional random spaces, along with specifics about the method for CFD [2]
• Comparison of speed-up over naive Monte Carlo sampling approaches [3]
• A worked example of stochastic collocation for a simple ocean model [4]

The example will be based on a nonlinear function of a single variable

(1)

The first step is to define a probability distribution for our input (Gaussian in this case, shown in Figure 1).

(2)

---

Figure 1: Uncertain Input Distribution

---

Then, since our example is a simple function, we can analytically calculate the resulting probability density for the output. To do that we need the derivative of the model with respect to the input (in the multi-variate case we’d need the Jacobian).

(3)

Then we divide the input density by the slope to get the new output density.

(4)

We’ll use equation 4 to measure the convergence of our collocation method.

The idea of stochastic collocation methods is that the points which the model (equation 1 in our example) is evaluated at are chosen so that they are orthogonal with the probability distributions on the inputs as a weighting function. For normally distributed inputs this results in choosing points at the roots of the Hermite polynomials. Luckily, Scipy has a nice set of special functions, which includes he_roots to give us the roots of the Hermite polynomials. What we really want is not the function evaluation at the collocation points, but the slope at the collocation points. One way to get this information in a nearly non-intrusive way is to use a complex step method (this is that little twist I mentioned). So, once we have values for this slope at the collocation points we fit a polynomial to those points and use this polynomial as a surrogate for equation 3. Figure 2 shows the results for several orders of surrogate slopes.

---

Figure 2: Stochastic Collocation Estimates of Output Uncertainty

---

The convergence of the method is shown in Figure 3. For low dimensional random spaces (only a few random inputs) this method will be much more efficient than random sampling. Eventually the curse of dimensionality takes its toll though, and random sampling becomes the faster approach.

---

Figure 3: Convergence of the Collocation Method

---

The cool thing about this method is that the complex step method which I used to estimate the Jacobian is also useful for sensitivity quantification which can be used for design optimization as well. On a related note, Dan Hughes has a post describing an interesting approach to bounding uncertainties in the output of a calculation given bounds on the inputs by using interval arithmetic.

References

[1]   Loeven, G.J.A., Witteveen, J.A.S., Bijl, H., Probabilistic Collocation: An Efficient Non-Intrusive Approach For Arbitrarily Distributed Parametric Uncertainties, 45th Aerospace Sciences Meeting, Reno, NV, 2007.

[2]   Hosder, S., Walters, R.W., Non-Intrusive Polynomial Chaos Methods for Uncertainty Quantification in Fluid Dynamics, 48th Aerospace Sciences Meeting, Orlando, FL, 2010.

[3]   Xiu, D., Fast Numerical Methods for Stochastic Computations, Comm. Comp. Phys., 2009.

[4]   Webster, M., Tatang, M.A., McRae, G.J., Application of the Probabilistic Collocation Method for an Uncertainty Analysis of a Simple Ocean Model, MIT JPSPGC Report 4, 2000.

Complex Step Jacobian-Free Newton-Krylov

Newton's method are really useful for solving non-linear systems.

Each iteration of the Newton's method requires solving a linear equation:

where the Jacobian is the matrix of partial derivatives

and then applying the update to the solution vector:

For small systems it is tractable to analytically calculate the Jacobian and invert it directly. This rapidly becomes intractable for realistic engineering problems because the number of unknowns can be on the order of 1e6 (if you have a good workstation) to 1e9 (if you have a big cluster). Rather than explicitly creating a 1e6-squared matrix and inverting it at every iteration, an iterative method that only requires the action of the matrix on a vector can be used (such as conjugate gradient or GMRES). Then we can approximate this matrix vector product by

which is a good improvement in storage over the matrix approach, but we can improve it still further by using a complex step

This means we only need one system evaluation to approximate the Jacobian, and much more accurate estimates can be made. As shown in a previous post it is second order accurate, and not as limited by subtractive cancellation as the standard difference approach is.

Here's a simple example of the idea in Maxima

/* trying out complex step finite differences for jacobian free newton
krylov methods */

/* define a simple system of equations */
F : zeromatrix(2,1);
F[1,1] : x[1]*x[2] + 2*x[2];
F[2,1] : 2*x[1] + x[1]*x[2];

/* analytic Jacobian */
J : zeromatrix(2,2);
J[1,1] : diff(F[1,1],x[1]);
J[1,2] : diff(F[1,1],x[2]);
J[2,1] : diff(F[2,1],x[1]);
J[2,2] : diff(F[2,1],x[2]);

u : zeromatrix(2,1);
v : zeromatrix(2,1);
u[1,1] : 1.5;
u[2,1] : 1.0;
v[1,1] : 0.2;
v[2,1] : 0.1;

c_step : 0.0 + 1.0e-200*%i; /* complex step */

/* approximate Jacobian using complex step */
Jv_approx : (
ev(F, x[1] = u[1,1] + c_step*v[1,1], x[2] = u[2,1] + c_step*v[2,1])
) / imagpart(c_step);
imagpart(Jv_approx);

/* evaluate the analytic Jacobian to compare */
Jv_analyt : ev(J.v,x[1]=u[1,1],x[2]=u[2,1]);

A good survey article introducing the JFNK method and describing some applications:
Jacobian-free Newton-Krylov methods: a survey of approaches and applications

Nice short paper about using complex steps in finite differences:
Using Complex Variables to Estimate Derivatives of Real Functions

Asymmetric Complex Step

To solve Laplace's equation with a complex step finite difference we needed to develop asymmetric difference relations to solve on the +/- imaginary planes.

On the plane in the negative imaginary direction we difference centrally in the real direction and forward two steps in the imaginary direction:

And on the plane in the positive imaginary direction we difference centrally in the real direction and backward two steps in the imaginary direction:

We still want the real part of the answer as we did before, but now we have imaginary coefficients for our function values in the real direction. While we have a compact stencil on the real plane (only nearest neighbors are used), we lose that property on either imaginary plane of the solution. This means our system matrix bandwidth is larger than it would be if we had a purely compact stencil for all parts of the solution.

It's fairly easy to calculate arbitrary derivative approximation formulas using Maxima.

/* calculate the coefficients for a finite difference approximation of
a derivative, can include complex steps */

steps : matrix([1,2,%i,2*%i]); /* multiple of the step size h */

/* sum df d2f d3f d4f */
derivs : matrix([0], [1], [0], [0], [0]); /* the rhs for our system */

n : matrix_size( steps );

A : zeromatrix(n[2]+1,n[2]+1);

for j : 1 thru n[2]+1 do
( A[1,j] : 1,
for i : 2 thru n[2]+1 do
(
if j = 1 then A[i,j] : 0 else A[i,j] : steps[1,j-1]^(i-1)/(i-1)!
)
);

invA : invert(A);

coeffs : expand( invA.derivs );
expand( 20*coeffs );

Where the vector steps tells which points will be included in the stencil, and the vector derivs tells which derivative we want to approximate. This example is calculating a real and complex forward difference approximation for the first derivative.

You can download this Maxima batch file here. This example can be run from the command line:

$ maxima -b fd.mac

or using the batch command from an interactive Maxima session:

(%i1) batch("fd.mac")

Complex Step: Subtractive Cancellation

We showed in an earlier post that approximating first derivatives with complex steps leads to a second order accurate formula that avoids the subtractive cancellation associated with a normal finite difference.

Following the simple convergence study for second derivative approximations we did for cos(). We'll look at the convergence of some first derivative approximation formulas.

The simple complex step approximation is given by

The standard second order finite difference is given by

The 4th order combined complex and real step approximation is

Notice that this formula is not able to avoid subtractive cancellation the way the second order approximation does, but we get a fourth order accurate approximation on a compact stencil.

The convergence of these three formulas for approximating the derivative of cosine is shown in the figure below.

The 4th order scheme runs into subtractive cancellation earlier than the standard 2nd order scheme, but it remains more accurate over the examined range of resolution. The complex step approximation continues to converge as the step size decreases even after the other two formulas diverge.

So for doing sensitivity analysis the second order complex step formula should be used, because only one step needs to be taken, and arbitrarily high precision can be achieved. If you need to approximate a derivative over a region, then the higher order method will pay-off because you can get better accuracy with fewer points.

Is it really O(h**4)?

Over the course of a couple posts (1 and 2) we've shown that you can extend the idea of finite difference approximations to include differences in the imaginary direction. That is you can take complex steps. I presented a formula for the second derivative on a compact stencil (using only nearest neighbours) which I claimed was fourth order accurate. But is it really?

Here's a simple convergence study using cos(), whose second derivative is -cos(). We'll use python to loop over an array with the number of points per wavelength we want to use in our approximation and then calculate a residual, or error, between the analytical solution and our approximation.

for ni in n:
x = 2*pi*arange(0,ni)/(ni-1)
h = x[1]-x[0]
f = cos(x)
analytical_d2f = -cos(x)
# take a positive imaginary step
x_step = x + complex(0,h)
f_pc = cos(x_step)
# take a negative imaginary step
x_step = x - complex(0,h)
f_mc = cos(x_step)
# take a positive real step
x_step = x + h
f_ph = cos(x_step)
# take a negative real step
x_step = x - h
f_mh = cos(x_step)

approx_d2f_c = ( f_ph + f_mh - f_pc.real - f_mc.real)/(2*h*h)
approx_d2f_h = ( f_ph + f_mh - 2*f )/(h*h)

resid_c[i] = sqrt( sum( (analytical_d2f - approx_d2f_c)**2 ) )
resid_h[i] = sqrt( sum( (analytical_d2f - approx_d2f_h)**2 ) )
i += 1

The approximated and analytical second derivatives are shown below for 10 points per wavelength using both a standard central difference with only a real step as well as the complex step formula.

The convergence behaviour for the two approximations is shown below.

Obviously we're getting much better convergence with our complex step formula.

Complex Step Laplace

We can use the results from the previous post about complex step finite differences to solve the 2D Laplace's equation. Instead of solving the equation on a single rectangular grid as we would with normal finite differences, we now solve on three parallel grids (shown below).

One of the grids lies on the real x-y plane, one of the grids is offset slightly in the positive imaginary direction (blue arrow in the figure above), and one of the grids is offset slightly in the negative imaginary direction. We care about the real part of those two imaginary solutions for our second derivative approximation:


The central difference takes care of our solution on the real grid, but what about the two imaginary grids? For those we need asymmetric differences, and I'll put those up in another post.

You might recognize the look of the example mesh if you've ever played around with the 3D modeling/animation program Blender.

Complex Step**2

I wrote a while back about using a complex step to approximate the first derivative of a function. The neat thing being that there was no subtractive cancellation in the approximation formula as there would be in a normal finite difference. This meant that we could choose very small step sizes without suffering the effects of magnified round-off error.

What about higher order derivatives? The simple result for the first derivative was a second order approximation that only required one function evaluation. Can we achieve the same benefits for a second derivative? The answer is not quite, but it's still pretty good.

The usual way to approximate a second derivative with finite differences is to take a central difference. This uses Taylor series expansions about the location of interest in two directions:


Adding the idea of using complex steps gives two additional points to include in the approximation formula:


Now we just have to solve a linear system to figure out what the coefficients are in our approximation formula:

Based on the Taylor series expansions above, we want the real part of the answer (rather than the imaginary part as was the case with the first derivative). The coefficients are 1/2 for the functions with the real step and -1/2 for the functions with the imaginary step (and zero for the function at x). So we weren't able to avoid subtractive cancellation, but we got an explicit 4th order accurate approximation for the second derivative on a compact stencil!

Complex Step

Finite differences are cool, but you are limited by subtractive cancellation (section 3.4: subtractive cancellation exercises). What if you want to make a numerical estimate of the derivative of a function, like a 2D Gaussian:


A finite difference approximation would be (in Octave):

dfdx = ( f(x+h) - f(x) )/ h;

Complex step is even cooler, because you don't have any subtraction (there's no difference), so you can choose a very small step size without loosing accuracy due to subtractive cancelation:

complex_step = complex( 0, 1e-320 );

The derivative is approximated by just the imaginary part:

x = x + complex_step;
dfdx = imag( f(x) )/imag( complex_step );

So the derivative with respect to x looks like this:

And the estimated derivative with respect to y:


This approach is really useful for design sensitivity analysis, and since modern Fortran supports complex types we can even use this method for serious number crunching!