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 !


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.