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
?