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 and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.