Convergence for Falkner-Skan Solutions

About 6 months ago Dan Hughes sent me a link to an interesting paper on "chaotic" behavior in the trajectory of iterates in a numerical Falkner-Skan solution. It struck me that the novel results reported in that paper were an artifact of the numerical method, and had little to do with any "chaotic" physics that might be going on in boundary layers or other systems that might be well described by this equation. This is similar to the point I made in the Fun with Filip post: the choice of numerical method matters. Do not rush to judgment about problems until you have brought the most appropriate methods to bear.

There are some things about the paper that are not novel, and others that seem to be nonsense. It is well-known that there can be multiple solutions at given parameter values (non-uniqueness) for this equation, see White. There is the odd claim that "the flow starts to create shock waves in the medium [above the critical wedge angle], which is a representation of chaotic behavior in the flow field." Weak solutions (solutions with discontinuities/shocks) and chaotic dynamics are two different things. They use the fact that the method they choose does not converge when two solutions are possible as evidence of chaotic dynamics. Perhaps the iterates really do exhibit chaos, but this is purely an artifact of the method (i.e. there is no physical time in this problem, only the pseudo-time of the iterative scheme). By using a different approach you will get different "dynamics", and with proper choice of method, can get convergence (spectral even!) to any of the multiple solutions depending on what initial condition you give your iterative scheme. They introduce a parameter, \(\eta_{\infty}\), for the finite value of the independent variable at "infinity" (i.e. the domain is truncated). There is nothing wrong with this (actually it's a commonly used approach for this problem), but it is not a good idea to solve for this parameter as well as the shear at the wall in your Newton iteration. A more careful approach of mapping the boundary point "to infinity" as the grid resolution is increased (following one of Boyd's suggested mappings) removes the need to solve for this parameter, and gives spectral convergence for this problem even in the presence of non-uniqueness and the not uncommon vexation of a boundary condition defined at infinity (all of external aerodynamics has this helpful feature).

Chaotic Method of Manufactured Solutions: Lorenz 63

I quickly brushed over the method of manufactured solutions (MMS) aspects in the previous post on Lorenz 63, you might wonder why. Well, it turns out that applying the MMS to chaotic systems includes an additional consideration that is not normally of concern. Here is the ’standard’ sort of checklist for choosing manufactured solutions (taken from Verification of Computer Codes in Computational Science and Engineering, emphasis mine):

1. Manufactured solutions should be sufficiently smooth on the problem domain so that the theoretical order-of-accuracy can be matched by the observed order-of-accuracy obtained from the test. A manufactured solution that is not sufficiently smooth may decrease the observed order-of-accuracy (however, see # 7 in this list).
2. The solution should be general enough that it exercises every term in the governing equation. For example, one should not choose temperature T in the unsteady heat equation to be independent of time. If the governing equation contains spatial cross-derivatives, make sure the manufactured solution has a non-zero cross derivative.
3. The solution should have a sufficient number of nontrivial derivatives. For example, if the code that solves the heat conduction equation is second order in space, picking T in the heat equation to be a linear function of time and space will not provide a sufficient test because the discretization error would be zero (to within round-off) even on coarse grids.
4. Solution derivatives should be bounded by a small constant. This ensures that the solution is not a strongly varying function of space or time or both. If this guideline is not met, then one may not be able to demonstrate the required asymptotic order-of-accuracy using practical grid sizes. Usually, the free constants that are part of the manufactured solution can be selected to meet this guideline.
5. The manufactured solution should not prevent the code from running successfully to completion during testing. Robustness issues are not a part of code order verification. For example, if the code (explicitly or implicitly) assumes the solution is positive, make sure the manufactured solution is positive; or if the heat conduction code expects time units of seconds, do not give it a solution whose units are nanoseconds.
6. Manufactured solutions should be composed of simple analytical functions like polynomials, trigonometric, or exponential functions so that the exact solution can be conveniently and accurately computed. An exact solution composed of infinite series or integrals of functions with singularities is not convenient.
7. The solution should be constructed in such a manner that the differential operators in the governing equations make sense. For example, in the heat conduction equation, the flux is required to be differentiable. Therefore, if one desires to test the code for the case of discontinuous thermal conductivity, then the manufactured solution for temperature must be constructed in such a way that the flux is differentiable. The resultant manufactured solution for temperature is nondifferentiable.
8. Avoid manufactured solutions that grow exponentially in time to avoid confusion with numerical instability.

As emphasized in # 4 we often have to use the free constants in our manufactured solution to bound the derivatives in the domain. For chaotic systems we’ll need to use not only the free constants, but the initial conditions as well, to achieve another end: ensure that the chosen solution is in the convergent region at all times. This is a sufficient but not necessary condition for the source term to be in the basin of entrainment. When the source term is in the basin of entrainment then we can expect to be able to verify the ordered convergence of our implementation as we would in the non-chaotic case. Otherwises we can expect the truncation error to eventually lead to an arbitrarily large departure from the manufactured solution. This terminology comes from the nonlinear dynamics and control literature [1] [3]. A worked example for the Lorenz ’63 system showing the boundary of the convergent region is in Chapter 6 of [2].

The convergent region is defined by the region of phase space for which the eigenvalues of the system Jacobian has negative real parts. As in any other application of the MMS, we start by choosing some convenient solution consistent with the desiderata enumerated above. Trigonometric functions are convenient in this case:

(1)

(2)

(3)

Now we need the Jacobian of the system (and then its eigenvalues), since our governing equations are definied in Maxima, that is easy enough:

/* Need the Jacobian of the system to find the ’convergent region’ so 
   we can select a useful IC/amplitude for our manufactured solutions */ 
sys_eqns : [eqn_1, eqn_2, eqn_3] $ 
J_cr : zeromatrix(3,3) $ 
for j : 1 thru 3 do 
( for i : 1 thru 3 do 
  ( 
    J_cr[i,j] : diff(rhs(sys_eqns[i]), unknowns[j]) 
    ) 
  ) $ 
[J_cr_eigvals, J_cr_eigmult] : eigenvalues(J_cr) $ 
J_cr_eigvals : fullratsimp(J_cr_eigvals) $ $

As before these expressions are output in Fortran90 format for compilation. Here are the eigenvalues in the complex plane for b_x = 1.0, b_y = -1.0, b_z = 10.0.

---

Figure 1: Jacobian eigenvalues for manufactured solution, bx = 1.0, by = -1.0, bz = 10.0

---

We can immediately see the problem from Figure 1: our manufactured solution is in the region of phase space where one of the Jacobian’s eigenvalues has a positive real part. We can adjust the initial conditions of our solution to get all of our eigenvalues on the left half-plane.

---

Figure 2: Jacobian eigenvalues for manufactured solution, bx = 10.0, by = -10.0, bz = 40.0

---

Being limited to only certain portions of the phase space for conducting verification is fine, because a single well-designed manufactured solution verifies the correctness of the entire implementation (as long as all of the terms are activated). This is the same reason that it is often emphasized in the V&V literature that there is no need for the manufactured solution to be physically meaningful.

References

[1]   Jackson, E.A., “Controls of dynamic flows with attractors,” Physics Review A, Vol. 44, Is. 8, 1991.

[2]   Wu, W., “Analytical and Numerical Methods Applied to Nonlinear Vessel Dynamics and Code Verification for Chaotic Systems,” Disertation, Virginia Polytechnic Institute and State University, Dec, 2009.

[3]   W. Wu, L.S. McCue, and C.J. Roy. “The method of manufactured solutions applied to chaotic systems,” Nonlinear Dynamics, 2009, under review.

Lorenz 63 Ensemble

Dan Hughes has started an ambitious series of posts about formally verifying (using the method of manufactured solutions (MMS)) a range of chaotic ordinary differential equations (ODE) systems.

One of the classic ones he’s looking at is Lorenz’s three-dimensional toy system from 1963:

(1)

(2)

(3)

where the system exhibits chaotic behaviour depending on the choice of parameter values (Ra = 28, Pr = 10.0, b = 8∕3 will be used throughout).

The first thing to do is define our governing equations in a computer algebra system (I like Maxima) because that makes writing bug free numerical software and complicated MMS forcing functions much easier.

  unknowns : [x,y,z] $ 
  depends(unknowns, t) $ 
 
  /* define the governing equations: */ 
  expr_1 : -Pr * x + Pr * y - diff(x,t) $ 
  eqn_1 : solve(expr_1, diff(x,t))[1] $ 
  expr_2 : -y + Ra * x - x * z - diff(y,t) $ 
  eqn_2 : solve(expr_2, diff(y,t))[1] $ 
  expr_3 : -b * z + x * y - diff(z,t) $ 
  eqn_3 : solve(expr_3, diff(z,t))[1] $

Then we need to define our discrete approximation to the governing equations (I’ll skip the MMS stuff for now)

/* elements of the Butcher tableau: */ 
A : matrix([1/4, -1/4], 
           [1/4, 5/12]) $ 
bT : matrix([1/4, 3/4]) $ 
c : matrix([0], [2/3]) $ 
 
/* identity matrix for doing Kronecker products, dimension depends on 
   the size of the system you are trying to solve */ 
sys_I : diagmatrix(3, 1) $ 
 
K : transpose([k[1], k[2], k[3], k[4], k[5], k[6]]) $ 
Xn : transpose([x[n],y[n],z[n],x[n],y[n],z[n]]) $ 
tn : transpose([t[n], t[n], t[n], t[n], t[n], t[n]]) $ 
F : transpose([f[1], f[2], f[3], f[4], f[5], f[6]]) $ 
 
A_rk : kronecker_product(A, sys_I) $ 
bT_rk : kronecker_product(bT, sys_I) $ 
hc_rk : kronecker_product(t[n] + h*c, sys_I) $ 
 
/* argument to our system right-hand-side operator: */ 
nonlin_arg : Xn + h * A_rk . K $ 
/* set up the two-stage system: */ 
rk_sys : zeromatrix(6,1) $ 
q_sys : zeromatrix(6,1) $ 
rk_sys[1,1] : k[1] - subst(x=nonlin_arg[1,1], subst(y=nonlin_arg[2,1], 
    subst(z=nonlin_arg[3,1], rhs(eqn_1)))) $ 
rk_sys[2,1] : k[2] - subst(x=nonlin_arg[1,1], subst(y=nonlin_arg[2,1], 
    subst(z=nonlin_arg[3,1], rhs(eqn_2)))) $ 
rk_sys[3,1] : k[3] - subst(x=nonlin_arg[1,1], subst(y=nonlin_arg[2,1], 
    subst(z=nonlin_arg[3,1], rhs(eqn_3)))) $ 
rk_sys[4,1] : k[4] - subst(x=nonlin_arg[4,1], subst(y=nonlin_arg[5,1], 
    subst(z=nonlin_arg[6,1], rhs(eqn_1)))) $ 
rk_sys[5,1] : k[5] - subst(x=nonlin_arg[4,1], subst(y=nonlin_arg[5,1], 
    subst(z=nonlin_arg[6,1], rhs(eqn_2)))) $ 
rk_sys[6,1] : k[6] - subst(x=nonlin_arg[4,1], subst(y=nonlin_arg[5,1], 
    subst(z=nonlin_arg[6,1], rhs(eqn_3)))) $

This is basically the same approached that I applied to the simple harmonic oscillator, but this time our system is non-linear, so we can’t invert the time-integration operator directly. We need to use an iterative Newton’s method to solve for each time-step, for which we need to define a Jacobian:

/* calculate the system Jacobian for use in a Newton’s method: */ 
J : zeromatrix(6,6) $ 
for j : 1 thru 6 do 
( for i : 1 thru 6 do 
  ( 
    J[i,j] : diff(rk_sys[i,1], k[j]) 
    ) 
  ) $ 
 
time_update : factor(h* bT_rk . K) $ $

The really nice thing about putting in a bit of work up-front with the symbol manipulation software is you can just dump the resulting expressions to compile-able code:

load(f90) $ 
file_output_append : false $ 
fname : ”rhs.f90” $ 
with_stdout(fname, print(”! generated by lorentz63_mms.mac”)) $ 
file_output_append : true $ 
with_stdout(fname, f90(f[1] = rk_sys[1,1])) $ 
/* ...much more boring file IO follows... */

With all of our useful expressions safely (and correctly) in Fortran, we can compile with F2py and play with them in Python. Integrating the system with a third-order implicit Runge-Kutta method gives the following time-step convergence behaviour:

---

(a) X component (b) Y component (c) Z component Figure 1: Time-step convergence for Lorenz ’63 system with 2-stage Radau-IA RK scheme

---

Figure 1 shows that as the time-step is decreased the trajectory stays with the ’pack’ for longer and longer times, but when the answers do diverge significantly the change is quite abrupt.

The most practically interesting thing with chaotic systems is the high sensitivity to initial conditions. This sensitivity competes with our uncertain knowledge of the true state of things in determining the limit on our ability to predict the trajectory of (some) dynamic systems (even if our knowledge were perfect, the finite precision of our machines still causes us significant trouble). Part of the way we deal with uncertainty is to do averaging over ensembles. To illustrate this idea I’ll start 7 realizations of the previous calculation, each adjusted by a small amount (on the order of the floating point round-off error for my machine).

nt = 400000 
T = 40.0 # integration period 
Pr = 10.0 
Ra = 28.0 
b = 8.0 / 3.0 
# knobs for the Newton’s method: 
tol = 1e-12 
maxits = 20 
nic = 7 # number of different IC’s in our ensemble 
h = T / float(nt) 
x = sp.zeros((nt,nic), dtype=float, order=’F’) 
y = sp.zeros((nt,nic), dtype=float, order=’F’) 
z = sp.zeros((nt,nic), dtype=float, order=’F’) 
t = sp.linspace(0.0, T, nt) 
# some perturbations: 
eps = 1e-16 * sp.array([0, 1, -1, 2, -2, 3, -3]) 
for i in xrange(nic): 
    x[0][i] = 1.0 + eps[i] 
    y[0][i] = -1.0 + eps[i] 
    z[0][i] = 10.0 + eps[i]

I’ll then calculate the pair-wise difference between each trajectory (thats 7 choose 2 or 21 differences). The combination function from the itertools module comes in handy here

from itertools import combinations 
from scipy.misc import comb 
 
# calculate our ensemble of trajectories: 
for i in xrange(nic): 
    L63.time_loop(x[:,i],y[:,i],z[:,i],t,Ra,Pr,b,h,tol,maxits) 
 
# calculate all the pair-wise differences between trajectories (nic 
# choose 2 of em): 
npairs = comb(nic, 2, exact=1) 
xdiffs = sp.zeros((nt,npairs), dtype=float, order=’F’) 
ydiffs = sp.zeros((nt,npairs), dtype=float, order=’F’) 
zdiffs = sp.zeros((nt,npairs), dtype=float, order=’F’) 
for i, pair in enumerate(combinations(xrange(nic), 2)): 
    xdiffs[:,i] = x[:,pair[0]] - x[:,pair[1]] 
    ydiffs[:,i] = y[:,pair[0]] - y[:,pair[1]] 
    zdiffs[:,i] = z[:,pair[0]] - z[:,pair[1]]

Here’s the resulting differences:

---

(a) X component (b) Y component (c) Z component Figure 2: Differences between members of the ensemble

---

You can see that things go along well in the early time, but eventually the differences between ensemble members grows suddenly to be on the same order of magnitude as the things we care about predicting.

Radau-type Runge-Kutta

Implicit Runge-Kutta methods are interesting because you can get high-order accuracy with relatively few stages. A method based on Radau quadrature that is third order accurate and L-stable has only two stages. A comparable multi-step method based on a backward difference that requires saving two steps is only second order accurate.

The Butcher tableau for the method is




The method is applied by calculating several slopes at various points and then adding a weighted sum of them to the initial condition.

The stability of the method can be understood by taking a look at the amplification factor, which for this Runge-Kutta method is given by

The amplification factor in the complex plane for two third-order A-stable Radau-based methods are shown below

The one on the left (which is the one we're implementing) exhibits stiff decay, otherwise known as L-stability. The amplification factor goes to zero as the max eigenvalue scaled by the time step goes to infinity. We can confirm this with Maxima:

limit((2*z+6)/(z**2-4*z+6), z, minf, plus);

L-stable methods are particularly important for stiff problems. These often occur in chemical kinetics or turbulent flows. Also, Chebyshev pseudospectral methods for PDEs introduce really large eigenvalues which are not very accurately approximated, so damping them is advantageous.

It's easy to apply the method to a linear ODE, such as the second order equation describing the motion of a simple harmonic oscillator

which has the analytical solution

and can be written as a first order system

We can use Maxima's kronecker_product() function to calculate the system operator for both stages

and similarly the operators for the Runge-Kutta integration


Since the system is linear we can actually invert the operator analytically and solve for the update.

Which in our case gives the solution at the new time step as

Of course the reason this type of time-integration is not very popular is because often you can't analytically invert the system. For nonlinear problems the size of the system that must be iteratively solved is twice as large (for a two-stage method) as the previously mentioned second order backward difference method.

We can dump the update expression to Fortran using Maxima's f90() function, and then compile it to a Python module using F2Py. This makes doing a convergence study in Python pretty easy:

f = 1.0 # frequency
w0 = 2.0*sp.pi*f # angular frequency
np = 1.0 # number of periods
nt = [16, 32, 64, 128, 256] # number of time-steps
h = []
x = []
dxdt = []
t = []
xa = []
err = []
for i in xrange(5):
h.append(np * (1.0/f) / nt[i])
x.append(sp.ones(nt[i], dtype=float))
dxdt.append(sp.zeros(nt[i], dtype=float))
t.append(sp.linspace(0, np-h[i], nt[i]))
sho.radau3(x[i], dxdt[i], h[i], w0)
xa.append(sho.analytical(x[i][0], dxdt[i][0], w0, t[i]))
err.append(sp.sqrt(sum((xa[i] - x[i])**2)/nt[i]))

The numerical solution is compared to the analytical solution for several time-step sizes. The observed order of convergence (2.97) matches the design order (3) pretty closely.

Check out Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations and Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems for some pretty detailed coverage of these and other advanced integration methods. Also see this page with Fortran implementations of several of these types of methods.

Derivative of Unequally Spaced Points

I posted a while back on solving Burgers' equation on a moving grid. The same coordinate transformation techniques applied there can be used to calculate the derivative of unequally spaced data points. Just as in the moving grid application, the grid transformation takes the points in the physical space to the computational space. Once in computational space, the sky's the limit!

One way to do this would be just a simple, low-order finite difference for the first derivative, maybe using an average of two one-sided differences (that would reduce to a central difference if the points were equally spaced). We can easily do better than that though since there are so many great, optimized transforms available in the FFTW library. Spectral methods give results that are accurate down to the precision of the machine for more than about 20 points (for functions without discontinuities).

The cost of using a grid-transformation is that we have to calculate two derivatives and then multiply them together, and the cost of a spectral method is O(n*log(n)) instead of O(n) as in a finite difference approximation. In one-dimension things are pretty easy, the chain rule shows us the two derivatives we need to calculate:

It is not straight-forward to directly calculate the additional derivative, so we calculate it by inverting the derivative of the physical coordinate with respect to the computational one:

Using the discrete cosine transform from FFTW (Chebyshev psuedospectral method) is straightforward to do by using F2Py and Python:

n = 17 # number of data points
a = 0.0 # left side of the interval
b = 3 * np.pi # right side of the interval
alpha = 0.25 # parameter to scale the perturbation
dx = (b - a) / float(n) # this is the equally spaced delta x
half_interval = (b - a) / 2.0 # so the dct diff is properly normalized
n_perturb = 3 # number of perturbed sets to try

dx_perturb = alpha * dx * (sp.random.random_sample((n,n_perturb)) - 0.5)

x = np.linspace(a, b, n)
dxdxi_perturb = np.zeros((n,n_perturb),dtype=float)
dydxi_perturb = np.zeros((n,n_perturb),dtype=float)
y_perturb = np.zeros((n,n_perturb),dtype=float)
x_perturb = np.zeros((n,n_perturb),dtype=float)

for i in xrange(n_perturb):
x_perturb[:,i] = x + dx_perturb[:,i] # add in the random fluctuations

y = np.cos(x)
y_perturb = np.cos(x_perturb)
dxdxi = cs.dct_diff(x, half_interval)
dy_analytical = - np.sin(x)
dydxi = cs.dct_diff(y, half_interval)
# derivatives with respect to the computational coordinate:
for i in xrange(n_perturb):
dxdxi_perturb[:,i] = cs.dct_diff(x_perturb[:,i], half_interval)
dydxi_perturb[:,i] = cs.dct_diff(y_perturb[:,i], half_interval)

# now actually calculate the derivatives we're after:
dydx = dydxi / dxdxi
dydx_perturb = dydxi_perturb / dxdxi_perturb

The derivative of the physical coordinate with respect to the computational coordinate for the equally spaced case and a couple of cases with some random perturbations added:

The actual derivative that we're after:

This example shows the flexibility of the approach, but it also indicates the importance of smooth changes in grid spacing for accurate derivative approximations. Abrupt changes in grid spacing translates into noise in the derivative.

Central Schemes and Discontinuities

I've been posting lately on solving Burgers' equation on moving grids, but haven't said much about difficulties involved in properly approximating the non-linear term involving the spatial derivative. A naive approach to solving Burgers' equation with central differences leads to problems, termed "blow-up" by early practitioners doing climate modeling. Boyd discusses this problem as it relates to spectral methods.

A backward-time centered-space (BTCS) discretization of Burgers equation leads to a non-linear equation for each grid point that depends on it's nearest neighbors in the spatial coordinate and it's predecessor in the time coordinate. Backward time:

Centered space:

Resulting discrete Burgers' equation:

We can use a Newton method to solve the resulting system, and to do that we'll need the Jacobian:

To solve for the update in each iteration of the Newton's method, we'll generally need an expression for the action of the Jacobian on the update vector:

All of these symbolic manipulations of the governing equation can be done easily in Maxima (see here and here and here for examples).

Here's a graph of several time integration steps of the BTCS Burgers equation starting with part of a cosine:

This shows the characteristic behaviour of non-linear wave equations known as wave steepening, which causes the formation of shock waves (discontinuities). Unfortunately, if we take just a few more time-steps we get the beginnings of "blow-up":


Another interesting case is a small advecting step shown below.

In this case the size of the discontinuity is much smaller compared to the grid resolution, so the solution does not blow-up, but some weird things happen. Burgers' equation has one eigenvalue, or wave speed. All of the waves should travel in the plus x-direction, because u is positive everywhere in this case. That's not what we see though. We see the main step travel in the proper direction (after some smearing due to the first order time integration), but we see wiggly errors going in the opposite direction. Because the central differencing scheme does not respect the eigenvalues of the problem (the way an upwind scheme would), we introduced a new dispersive error wave that travels the wrong way.

Implicit-Time Burgers' Equation on a Moving Grid

In the last post on solving Burgers' equation on a moving grid we ended up with the semi-discrete equation for the rates:

where the spatial derivatives (for the solution and the grid) are approximated by simple central differences. We still haven't committed to a time-integration scheme. This is one of the nice things about the method of lines, we can play with different time integration schemes independent of our spatial discretization.

A good way to start out with implicit integration schemes is with the basic Backward Euler formula:

where the big X represents the unknown vector and the F(X) represents the right-hand side of our semi-discrete equation. The superscripts represent the time level. The reason this is an implicit time integration method is that the spatial terms are evaluated at the next time level rather than the present time level. This substitution is easily made in Maxima (just like the spatial derivatives):

/* substitute in a backward time difference: */
Discrete_Rate : ratsubst((u[i] - u[i,n-1])/dt, 'diff(u,t), Rate_Vector);
Discrete_Rate : ratsubst((x[i] - x[i,n-1])/dt, 'diff(x,t), Discrete_Rate);

Once the choice of time integration scheme is made, we have a fully-discrete problem that we can solve.

The Newton expression, G, in Maxima:

/* expression for Newton's method (zero finding): */
Newton_expression : fullratsimp(Discrete_Rate - Rate_expression);

Since G is non-linear, we'll need to use a Newton's method to solve for the zero (which will give us the solution and grid point locations at the next time-step). This requires the Jacobian:

/* Jacobian of the Newton expression: */
J : diff(Newton_expression, u[i-1]);
J : addcol(J, diff(Newton_expression, u[i]));
J : addcol(J, diff(Newton_expression, u[i+1]));
J : addcol(J, diff(Newton_expression, x[i-1]));
J : addcol(J, diff(Newton_expression, x[i]));
J : addcol(J, diff(Newton_expression, x[i+1]));

The update for each iteration of the Newton's method is given by:


Of course we generally don't invert the Jacobian directly to find the update vector, usually an iterative method such as Successive Over Relaxation or a Krylov subspace method is preferred. Most of these methods require a function which returns the application of the Jacobian on the update vector, so we go ahead and find that in Maxima as well:

/* Jacobian applied on the update vector: */
update_vector : transpose([du[i-1], du[i], du[i+1], dx[i-1], dx[i], dx[i+1]]);
Jdu : factor(J.update_vector);

The great thing about defining the numerical scheme in Maxima this way is that the f90() function will output these (increasingly complex) expressions to Fortran with the proper array indexing. If you've done this upfront part properly they are ready to go into loops as is.

Chebyshev Spectral Method Shock Filtering

Now that I've got a nice handy spectral derivative routine, time to test it out on some practical problems (ie non-smooth).

A useful scalar test function is a square wave with a little Gaussian blip added to it:

This is good test function because it contains discontinuities as well as small smooth features that we'd like to be able to resolve.

There are tons of options to filter out high frequency "ringing" around discontinuities (Gibbs phenomenon). One nice option is the Butterworth filter, the gain for which is:

Where the cut-off frequency and the order of the filter are free parameters you can "fiddle" with. For this set of derivatives I'm going to choose a cut-off frequency of half the highest frequency in my grid and see what the effect of the order parameter is. Here's a plot of the gain for a couple different order filters:

Notice they all have the same gain at the cut-off frequency, they just pass through it differently.

The un-smoothed spectral derivative compared to the derivative after different order Butterworth filters are applied to the Chebyshev coefficients before the inverse transform are shown below.

The presence of the discontinuities causes "ringing" throughout the entire grid, but the filters do a pretty good job of containing it to the area immediately around the shocks. Here's a zoomed in view of the smooth feature area:

Without filtering the derivatives of the smooth feature are totally swamped by ringing, but you actually get a nice result out after some filtering.