I am trying to study the stability of numerical discretization schemes using the Jacobian matrix of the residues with respect to the vector of conserved variables.
For a simple diffusion equation discretized with central difference scheme, the RHS can be written as:
$RHS=\nu\frac{u_{i+1}^{n}-2u_{i}^{n}+u_{i-1}^{n}}{\Delta x^2}$
This equation is interesting since the Jacobian matrix has an analytical form which can be written, in banded matrix notation as:
$J=\nu/dx^2 (+1,-2,+1)$
Basically, what I am trying to do is to compare the eigenvalues of the numerical approximated Jacobian with the eigenvalues of the analytical solution.
However, I am having some troubles because the results are not making sense. The numerically approximated Jacobian is presenting completly different eigenvalues when compared with the exact Jacobian. I am suspecting that the problem is in the way I am trying to compute the Frechet derivatives to compute the Jacobian. The code fragment below shows how I am doing so:
# Residue function of a simple centered three point scheme.
def frhs(um1,uo1,up1,dx,nu):
return nu*( (up1 - 2.0 * uo1 + um1)/dx**2.0 )
# Computes the derivative of the residues with respect to the solution vector.
eps = 0.00001
drhs_du = np.zeros((nx-1,nx-1)) # In order to take the eigenvalues, this shall be a matrix.
e = np.zeros(nx-1)
# This loop computes the jacobian matrix according to http://www.netlib.org/math/docpdf/ch08-04.pdf
for i in range(1,nx-1):
for j in range(1,nx-1):
e[j] = 1.0
drhs_du[i,j] = ( frhs(un[i-1]+eps*e[j], un[i]+eps*e[j], un[i+1]+eps*e[j],dx,nu) - frhs(un[i-1], un[i],un[i+1],dx,nu) )/eps
e[j] = 0.0
The complete code can be found in: https://github.com/lmarmotta/disc_stability/blob/master/1dDiffusionEq.py
Am I doing anything wrong ?
How would you guys approximate the Jacobian numerically ? And Are there any other algorithms to approximate Jacobian matrices for these type of cases ?
Thank You !