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()
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()
