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:
This equation is interesting since the Jacobian matrix has an analytical form which can be written, in banded matrix notation as:
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 !