Python: Differentiation matrix by numpy

differential2_matrix_frontpage

I made snippets for numerical differentiation by matrix.

I think one of the pros using matrix over for-loop is simplicity of code and speed.
Because numpy array is not recommended looping through array, differentiation by multiplying matrix and vector would suit for the proper usage.

But one of the cons using matrix is that it makes very sparse matrix. The resulted element number of matrix is length of list ** 2.
It can be converted into condensed format by scipy.sparse functionality, though.
Another aspect of cons is that making matrix for more complex differential is not intuitive.

By the way, I made functions to make differentiation matrix for 1st order and 2nd order.

 

1st order differentiation matrix:

Function to make differentiation matrix from a given vector:

 
def dif1_matrix(x):
    '''
    In case for len(x) = 5:
    dif = np.array(
    [[1, 0, 0, 0, 0],
     [-1, 1, 0, 0, 0],
     [0, -1, 1, 0, 0],
     [0, 0, -1, 1, 0],
     [0, 0, 0, -1, 1]])
    '''
    # Diagonal elements are 1.
    dif_now = np.diag(np.ones(len(x)))
    
    # Left elements of diagonal are -1.
    dif_pre_ones = np.ones(len(x)-1) * - 1 # -1 vector.
    dif_pre = np.diag(dif_pre_ones, k=-1) # Diagonal matrix shiftedto left.
    
    dif = dif_now + dif_pre
    return dif

Applying 1st order differentiation matrix:

import numpy as np
from matplotlib.pylab import plt
%matplotlib inline

x = np.linspace(-10,10,300)
dx = x[1] - x[0]
y = np.sin(x)

# Differentiate.
dif = dif1_matrix(y)
y_dif = dif.dot(y) / dx

plt.plot(x[1:], y_dif[1:], label="sin(x)'")
plt.plot(x[1:], y[1:], label="sin(x)")
plt.legend()

simplecode_math_differential2_matrix_1st

 

2nd order differentiation matrix:

Function to make differentiation matrix from a given vector:

def dif2_matrix(x):
    '''
    In case for len(x) = 5:
    dif = np.array(
    [[-2, 1, 0, 0, 0],
     [1, -2, 1, 0, 0],
     [0, 1, -2, 1, 0],
     [0, 0, 1, -2, 1],
     [0, 0, 0, 1, -2]])
    '''
    # Diagoenal elements are -2.
    dif_now = np.diag(np.ones(len(x))) * -2
    
    # Left elements of diagonal are 1.
    dif_pre_ones = np.ones(len(x)-1) # 1 vector.
    dif_pre = np.diag(dif_pre_ones, k=-1) # Diagonal matrix shiftedto left.

    # Right elements of diagonal are 1.
    dif_post_ones = np.ones(len(x)-1) # 1 vector.
    dif_post = np.diag(dif_post_ones, k=1) # Diagonal matrix shiftedto right.
    
    dif  =  dif_now + dif_pre + dif_post
    return dif

Applying differential matrix:

import numpy as np
from matplotlib.pylab import plt
%matplotlib inline

x = np.linspace(-10,10,300)
dx = x[1] - x[0]
y = np.sin(x)

# Differentiate.
dif2 = dif2_matrix(y)
y_dif2 = dif2.dot(y) / dx**2

plt.plot(x[1:-1], y_dif2[1:-1],  label="sin(x)''")
plt.plot(x[1:-1], y[1:-1],  label="sin(x)")
plt.legend()

simplecode_math_differential2_matrix_2nd

Leave a Reply