Demo notebook for linear operators as matrices
Demo notebook for linear operators as matrices#
In this notebook we explore implementing linear operators in Python as matrices (using numpy).
import numpy as np
First make a vector with specified spacing:
t = np.arange[0,10,1]
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
<ipython-input-2-f09977b052e8> in <module>
----> 1 t = np.arange[0,10,1]
TypeError: 'builtin_function_or_method' object is not subscriptable
Oops, I thought I was using Mathematica. Fix the brackets and check that we got what we wanted:
t = np.arange(0,10,1)
print(t)
[0 1 2 3 4 5 6 7 8 9]
Useful matrix functions include np.identity and np.ones:
print( np.identity(5) )
[[1. 0. 0. 0. 0.]
[0. 1. 0. 0. 0.]
[0. 0. 1. 0. 0.]
[0. 0. 0. 1. 0.]
[0. 0. 0. 0. 1.]]
print( np.ones(5) )
[1. 1. 1. 1. 1.]
print( np.ones((5,5)) )
[[1. 1. 1. 1. 1.]
[1. 1. 1. 1. 1.]
[1. 1. 1. 1. 1.]
[1. 1. 1. 1. 1.]
[1. 1. 1. 1. 1.]]
Ok, now try an operator that multiplies by \(\omega_0^2\):
omega0 = 2
D_op = omega0**2 * np.identity(10)
print(D_op)
[[4. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
[0. 4. 0. 0. 0. 0. 0. 0. 0. 0.]
[0. 0. 4. 0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 4. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 4. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 4. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 0. 4. 0. 0. 0.]
[0. 0. 0. 0. 0. 0. 0. 4. 0. 0.]
[0. 0. 0. 0. 0. 0. 0. 0. 4. 0.]
[0. 0. 0. 0. 0. 0. 0. 0. 0. 4.]]
Try it out! It is tempting to use * but for linear algebra we use @. Compare:
print(D_op * t)
[[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
[ 0. 4. 0. 0. 0. 0. 0. 0. 0. 0.]
[ 0. 0. 8. 0. 0. 0. 0. 0. 0. 0.]
[ 0. 0. 0. 12. 0. 0. 0. 0. 0. 0.]
[ 0. 0. 0. 0. 16. 0. 0. 0. 0. 0.]
[ 0. 0. 0. 0. 0. 20. 0. 0. 0. 0.]
[ 0. 0. 0. 0. 0. 0. 24. 0. 0. 0.]
[ 0. 0. 0. 0. 0. 0. 0. 28. 0. 0.]
[ 0. 0. 0. 0. 0. 0. 0. 0. 32. 0.]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 36.]]
print(D_op @ t)
[ 0. 4. 8. 12. 16. 20. 24. 28. 32. 36.]
Ok, let’s make it more general (note the difference between these two):
t_min = 0
t_max = 10
delta_t = 1
t = np.arange(t_min, t_max, delta_t)
print(t)
[0 1 2 3 4 5 6 7 8 9]
t_min = 0
t_max = 10
delta_t = 1
t = np.arange(t_min, t_max, delta_t)
print(t)
num_t = len(t)
print(num_t)
[0 1 2 3 4 5 6 7 8 9]
10
omega0 = 2
D_op = omega0**2 * np.identity(num_t)
print(D_op)
[[4. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
[0. 4. 0. 0. 0. 0. 0. 0. 0. 0.]
[0. 0. 4. 0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 4. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 4. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 4. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 0. 4. 0. 0. 0.]
[0. 0. 0. 0. 0. 0. 0. 4. 0. 0.]
[0. 0. 0. 0. 0. 0. 0. 0. 4. 0.]
[0. 0. 0. 0. 0. 0. 0. 0. 0. 4.]]
Now try the simplest derivative operator, building it from shifted diagonal matrices.
print( np.diag(np.ones(5), 0) )
[[1. 0. 0. 0. 0.]
[0. 1. 0. 0. 0.]
[0. 0. 1. 0. 0.]
[0. 0. 0. 1. 0.]
[0. 0. 0. 0. 1.]]
print( np.diag(np.ones(5), 1) )
[[0. 1. 0. 0. 0. 0.]
[0. 0. 1. 0. 0. 0.]
[0. 0. 0. 1. 0. 0.]
[0. 0. 0. 0. 1. 0.]
[0. 0. 0. 0. 0. 1.]
[0. 0. 0. 0. 0. 0.]]
print( np.diag(np.ones(5), -1) )
[[0. 0. 0. 0. 0. 0.]
[1. 0. 0. 0. 0. 0.]
[0. 1. 0. 0. 0. 0.]
[0. 0. 1. 0. 0. 0.]
[0. 0. 0. 1. 0. 0.]
[0. 0. 0. 0. 1. 0.]]
Diff_op = (1 * np.diag(np.ones(num_t-1), 1) + (-1) * np.diag(np.ones(num_t), 0)) / delta_t
print(Diff_op)
[[-1. 1. 0. 0. 0. 0. 0. 0. 0. 0.]
[ 0. -1. 1. 0. 0. 0. 0. 0. 0. 0.]
[ 0. 0. -1. 1. 0. 0. 0. 0. 0. 0.]
[ 0. 0. 0. -1. 1. 0. 0. 0. 0. 0.]
[ 0. 0. 0. 0. -1. 1. 0. 0. 0. 0.]
[ 0. 0. 0. 0. 0. -1. 1. 0. 0. 0.]
[ 0. 0. 0. 0. 0. 0. -1. 1. 0. 0.]
[ 0. 0. 0. 0. 0. 0. 0. -1. 1. 0.]
[ 0. 0. 0. 0. 0. 0. 0. 0. -1. 1.]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. -1.]]
Try it!
print(Diff_op @ t)
[ 1. 1. 1. 1. 1. 1. 1. 1. 1. -9.]
print(Diff_op @ t**2, '\n', 2*t)
[ 1. 3. 5. 7. 9. 11. 13. 15. 17. -81.]
[ 0 2 4 6 8 10 12 14 16 18]
Build a better derivative operator by making it symmetric:
Diff_sym_op = (1 * np.diag(np.ones(num_t-1), 1) + (-1) * np.diag(np.ones(num_t-1), -1)) / (2*delta_t)
print(Diff_sym_op)
[[ 0. 0.5 0. 0. 0. 0. 0. 0. 0. 0. ]
[-0.5 0. 0.5 0. 0. 0. 0. 0. 0. 0. ]
[ 0. -0.5 0. 0.5 0. 0. 0. 0. 0. 0. ]
[ 0. 0. -0.5 0. 0.5 0. 0. 0. 0. 0. ]
[ 0. 0. 0. -0.5 0. 0.5 0. 0. 0. 0. ]
[ 0. 0. 0. 0. -0.5 0. 0.5 0. 0. 0. ]
[ 0. 0. 0. 0. 0. -0.5 0. 0.5 0. 0. ]
[ 0. 0. 0. 0. 0. 0. -0.5 0. 0.5 0. ]
[ 0. 0. 0. 0. 0. 0. 0. -0.5 0. 0.5]
[ 0. 0. 0. 0. 0. 0. 0. 0. -0.5 0. ]]
print(Diff_sym_op @ t**2, '\n', 2*t)
[ 0.5 2. 4. 6. 8. 10. 12. 14. 16. -32. ]
[ 0 2 4 6 8 10 12 14 16 18]
print(Diff_sym_op @ t**3, '\n', 3*t**2)
[ 0.5 4. 13. 28. 49. 76. 109. 148. 193. -256. ]
[ 0 3 12 27 48 75 108 147 192 243]
Try with smaller spacing delta_t to get more accuracy.