Numerical Stability of a Generalized Spatial Discretization Scheme

Question

After reading the matrix stability chapter (10) of Hirsch [1], I decided to dive in the reference list of the chapter. One of the papers [2], which is cited as reference shows an very interesting approach to stability analysis.

The paper shows an stability evaluation of a real world flow solution of a NACA0012 airfoil with artificial dissipation and all the bells and witless which are expected to solve this type of flowfield.

Basically, the matrix stability method, as explained in Hirsch and, apparently used in Ref. [2], is defined as follows:

If you have an discretized PDE, you end up with.

$\frac{d \bf U }{d t} = \bf S \bf U + \bf Q$.

If you compute the eigenvalues of the Jacobian $\bf S$, the maximum real part of the eigenvalues say to you if you spatial discretization scheme is stable or not. A positive real eigenvalue amplifies in time and ends up messing up with your stability.

I am trying to study the stability of some very simple cases, however, I do not want to use the defined banded matrix definitions of the matrices $\bf S$. I want to be able to derive it numerically so, when the cases become less trivial, I will still be able to evaluate the stability of my case. The problem I am facing right now is how to derive the matrix $\bf S$, since it is defined as:

$ \bf S = \frac{\partial\, \bf RHS}{\partial \, \bf U} $

Where $\bf RHS$ represent my spatial discretization scheme, and $\bf U$ is my solution vector.

For a 1-D case, $\bf RHS$ is a vector with the same size as $\bf U$. I tried to use simple finite differences, however, they do not seen to be robust enough. Or, at least, trying to compute the derivatives with:

for i in range(1,n_points):
    drhs_du = (rhs[i+1] - rhs[i-1]) / 2.0*(u[i+1] - u[i-1])

gives me a whole lot of by-zero-divisions.

I saw some references using Frechet derivatives to compute the $\bf S$ Jacobian, but the authors who used it, which I know in person, also reported that this type of numerical derivative is quite troublesome.

My question, finally, is related to robust suggestions to compute this Jacobian numerically. Which algorithms would you use or, even, which derivative definition would you guys use ? Are there reliable high performance libraries with the algorithms implemented ?

Thank you !

[1] Hirsch.C "Numerical Computation of Internal and External Flows Vol.01"

[2] Eriksson.L.E. Rizzi "A Computer Aided Analysis of the Convergence to Steady State of Discrete Approximations to the Euler Equations"

Edit # 01

Can the proper implementation of the method described in the second paragraph of the comment made below be translated as:

def frhs(unn,unm1,dx):
    return (unn - unm1)/dx

eps = 0.0001

drhs_du = np.zeros((nx,nx))

for i in range(1,nx-1):
    for j in range(1,nx-1):
        drhs_du[i,j] = (frhs(un[i] + eps,un[i-1] + eps,dx) - frhs(un[j],un[j-1],dx))/eps

?

Answer 1

So it looks to me like there is a misunderstanding. your residual, or right hand side, presumably depends on multiple cells, so when building your jacobian, you should ask yourself how each entry in your residual depends on each cell. Now this should give you an nxn sparse matrix. If you in fact linearize this by hand, you should be able to store it in a vector that scales linearly with the number of cells. Given that you're 1d, I'm assuming it's size 3n, but that depends on your discretization. This is rather simple to do but rather tedious.

Another method is to calculate numerical derivatives, which you mentioned, and seem to be the source of your misunderstanding. It is not evaluated as you suggest, which is looking at different cells and taking a difference of the residuals and taking a difference of the cell values. Rather it is done through a cell by cell approach, where you perturb each cell's value individually and see how this affects the residual, and you take the difference between the residual from the perturbed and unperturbed states and divide by the perturbation. When doing this, you will discover that the stencil is pretty limited, and you should get a great deal of zeros as you fill out your matrix, and you will be able to store this sparse matrix in a vector as mentioned above.

The last possibility, is that you in fact do not need the jacobian matrix, you need the jacobian vector product. You can instead perturb the entire field by an infinitesimal complex perturbation in the direction of the vector you are multiplying the jacobian by, evaluate the residual of this perturbed state, and you will have obtained the frechet derivative without round-off error.