Copyright (C) 2020 Andreas Kloeckner
import numpy as np
Here's an upper-triangular matrix $A$ and two vectors $x$ and $b$ so that $Ax=b$.
See if you can find $x$ by computation.
n = 5
A = np.random.randn(n, n) * np.tri(n).T
print(A)
x = np.random.randn(n)
print(x)
b = A @ x
[[ 0.53359808 -1.66129881 0.31267643 -0.07384466 1.20957795] [-0. -1.1204435 -1.5348203 1.38270361 -0.34971611] [ 0. 0. -0.47187693 0.9763103 0.55054242] [ 0. -0. -0. -0.16929913 0.21209806] [-0. 0. 0. -0. -0.52165269]] [-0.9170418 1.40215838 1.41534372 -0.53305575 -1.02625922]
xcomp = np.zeros(n)
for i in range(n-1, -1, -1):
tmp = b[i]
for j in range(n-1, i, -1):
tmp -= xcomp[j]*A[i,j]
xcomp[i] = tmp/A[i,i]
Now compare the computed $x$ against the reference solution.
print(x)
print(xcomp)
[-0.9170418 1.40215838 1.41534372 -0.53305575 -1.02625922] [-0.9170418 1.40215838 1.41534372 -0.53305575 -1.02625922]
Questions/comments: