Matrices and Determinants
Contents
Matrices and Determinants¶
Matrices¶
In the previous module we saw vectors (arrays). A vector is actually a special case of a more general type of object called a matrix. A has “dimensions” of 1 row and N columns (a row vector) or N rows with 1 column (a column vector). We will now work with the general N row by M column Matrix.
As before there are a range of operations that can be used on Matrices. For a detailed discussion on Matrices see MathChapter G-F in Mcquarrie.
There are two ways to define a matrix.
Lay it out as you would if you were write it on paper.
A = [ 1 2 3
4 5 6
7 8 9]
Use semi-colons to separate rows
A = [ 1 2 3; 4 5 6; 7 8 9]
Operations on One Matrix¶
using LinearAlgebra
A = [0 1 -3; -3 -4 4;-2 -2 1 ]
trace_A = tr(A) # trace of matrix A
display("Matrix trace_A: ")
display(trace_A)
determiant_A = det(A) # determinant of matrix A
display("Matrix determiant_A: ")
display(determiant_A)
transpose_A = transpose(A) # transpose of A
display("Matrix transpose_A: ")
display(transpose_A)
inverse_A = inv(A) # invese of A
display("Matrix inverse_A: ")
println(inverse_A)
"Matrix trace_A: "
-3
"Matrix determiant_A: "
0.9999999999999991
"Matrix transpose_A: "
3×3 transpose(::Matrix{Int64}) with eltype Int64:
0 -3 -2
1 -4 -2
-3 4 1
"Matrix inverse_A: "
[4.000000000000005 5.000000000000006 -8.000000000000009; -5.0000000000000036 -6.000000000000004 9.000000000000007; -2.0000000000000013 -2.0000000000000018 3.0000000000000027]
using LinearAlgebra
A = [1 2 3; 4 5 6; 7 8 9]
row = 1
column = 3
value = A[row, column]
display("Matrix 1,3: ")
display(value)
display("Matrix 3,1: ")
display(A[3,1])
"Matrix 1,3: "
3
"Matrix 3,1: "
7
Below is some code that prints out all of the indecies in a matrix.
using LinearAlgebra
a = [1 2 3; 4 5 6; 7 8 9]
display(a)
b = ["1,1" "1,2" "1,3"; "2,1" "2,2" "2,3"; "3,1" "3,2" "3,3" ]
display(b)
3×3 Matrix{Int64}:
1 2 3
4 5 6
7 8 9
3×3 Matrix{String}:
"1,1" "1,2" "1,3"
"2,1" "2,2" "2,3"
"3,1" "3,2" "3,3"
using LinearAlgebra
a = [1 2 3; 4 5 6; 7 8 9]
display(a)
display(a[:,1])
display(a[2,:])
3×3 Matrix{Int64}:
1 2 3
4 5 6
7 8 9
3-element Vector{Int64}:
1
4
7
3-element Vector{Int64}:
4
5
6
Operations with two matrices:¶
using LinearAlgebra
Identity = [1 0 0; 0 1 0; 0 0 1]
A = [0 1 -3; -3 -4 4;-2 -2 1 ]
B = [0 1 2; 3 4 5;6 7 8]
display("Matrix A: ")
display(A)
Same_A = A * Identity # Multiply by an identity
Same_easier = A * I # I is a reserved built in matrix I, it will automatically scale to the correct size
display("Matrix Same_A: ")
display(Same_A)
display("Matrix Same_easier: ")
display(Same_easier)
A_times_B = A * B # Matrix Multiplication
A_Hadamard_B = A .* B # Element-wise (Hadamard product)
display("Matrix A_times_B: ")
display(A_times_B)
display("Matrix A_Hadamard_B: ")
display(A_Hadamard_B)
A_left_divide_B = A \ B # inverse(A) * B "left-division operator"
display("Matrix A_left_divide_B: ")
display(A_left_divide_B)
"Matrix A: "
3×3 Matrix{Int64}:
0 1 -3
-3 -4 4
-2 -2 1
"Matrix Same_A: "
3×3 Matrix{Int64}:
0 1 -3
-3 -4 4
-2 -2 1
"Matrix Same_easier: "
3×3 Matrix{Int64}:
0 1 -3
-3 -4 4
-2 -2 1
"Matrix A_times_B: "
3×3 Matrix{Int64}:
-15 -17 -19
12 9 6
0 -3 -6
"Matrix A_Hadamard_B: "
3×3 Matrix{Int64}:
0 1 -6
-9 -16 20
-12 -14 8
"Matrix A_left_divide_B: "
3×3 Matrix{Float64}:
-33.0 -32.0 -31.0
36.0 34.0 32.0
12.0 11.0 10.0
Eigen Values and Eigen Vectors¶
Julia will calculate the eigenvalues and eigenvectors of a matrix:
(The kth eigenvector can be obtained from the slice eigen(A).vectors[:, k].)
using LinearAlgebra
B = [0 1 2; 3 4 5;6 7 8]
eigen_of_B = eigen(B)
display("eigen values")
display(eigen_of_B.values)
display("eigen vectors")
display(eigen_of_B.vectors)
display("eigen vector 1")
display(eigen_of_B.vectors[:,1])
"eigen values"
3-element Vector{Float64}:
-1.3484692283495348
-2.484772788451793e-16
13.34846922834954
"eigen vectors"
3×3 Matrix{Float64}:
0.7997 0.408248 0.164764
0.104206 -0.816497 0.505774
-0.591288 0.408248 0.846785
"eigen vector 1"
3-element Vector{Float64}:
0.7996996631118652
0.10420578771917768
-0.5912880876735086
Using Matrices to solve systems of equations¶
Example:
You can turn this into a matrix equation:
A*X = b
A = [1 2 3 ; 2 3 1 ; 2 1 13]
b = [2; 1; 13]
x = A\b
3-element Vector{Float64}:
1.727272727272727
-1.0909090909090908
0.8181818181818182
We can then test to see if our result was correct
A = [1 2 3 ; 2 3 1 ; 2 1 13]
b = [2; 1; 13]
x = A\b
println("b[2]: $(b[1])")
computed_value = A[1,1] * x[1] + A[1,2]*x[2]+ A[1,3]*x[3]
println("computed value: $(computed_value)")
b[2]: 2
computed value: 2.0
This type of operation was exactly what the Atanasoff–Berry computer, built on the Iowa States Campus, was built to do in the 1940s. (Yes the 1940s, that very, very early days in digital computing.)
Exercises¶
Problem 1¶
Find A * B
Find Hadamard Product of A and B
Find 2A + 3B
Hint you can multiply a matrix by a scalar
Solve in: module-3-exercise-1.jl
Problem 3¶
Write a function that checks if a two dimensional (square NxN) matrix is symmetric. Your function should not have any duplicate checks.
i.e. test that for all combinations of i and j, M[i,j] = M[j,i]
Hint: the easiest way to do this will involve nested loops, but don’t loop over all i and all j.
Do not use LinearAlgebra.issymmetric() for this problem (but do use it in real life if you need to know if a matrix is symmetric)
Solve in: module-3-exercise-3.jl
Examples
Not Symmetric $\(A = \begin{bmatrix} -1 & 42 & 23 \\ 44 & -5 & 3 \\ 2 & -5 & 3 \end{bmatrix}\)$
Symmetric $\(A = \begin{bmatrix} 1 & 14 & -5 \\ 14 & 2 & 4 \\ -5 & 4 & 3 \end{bmatrix}\)$
Problem 4¶
Solve H-13 from Mcquarrie (Page 433). Solve the matrix equation and return the x vector with the solution. See section above titled: Using Matrices to solve systems of equations
a will be a parameter passed in
Solve in: module-3-exercise-4.jl
# include("module-3-exercise-1-test-runner.jl") #including the file runs the tests
# include("module-3-exercise-2-test-runner.jl") #including the file runs the tests
# include("module-3-exercise-3-test-runner.jl") #including the file runs the tests
# include("module-3-exercise-4-test-runner.jl") #including the file runs the tests
BONUS Topic Numerical Integration¶
Another application of series and loops (review from last time) are numerical techniques to calculate the area under a curve. One such technique using the midpoint rule divides the area under the curve between the integration limits a and b into n rectangles. The sum of the area of the rectangles is the approximation of the integral value. The larger n becomes, the closer to the exact integral. In the infinite limit, the values are the same. (This process is also called a Riemann Sum)
(source: wikipedia.org)
function function_to_be_integrated(x)
return exp(-2*x)
end
function numerical_integration(number_of_rectangles,
integration_start,
integration_end,
integration_function)
area_under_curve = 0
Δx = (integration_end-integration_start)/number_of_rectangles
for i in 0:number_of_rectangles
area_under_curve += integration_function(i*Δx)*Δx
end
return area_under_curve
end
n = 10000
a = 0
b = 1
## you can pass a function as a parameter to be used in another function.
area_under_curve = numerical_integration(n, a, b, function_to_be_integrated)
analytical_area_under_curve = 0.43233 # calculated with wolfram alpha
println("analytical: $(analytical_area_under_curve)")
println("numerical: $(area_under_curve)")
analytical: 0.43233
numerical: 0.4323891265869634
Problem 3 - Numerical Integration¶
Using similar steps as the numerical integration techniques above (Riemann Sum). Get the area under the curve for the normalized particle in a one-dimensional box wave function squared. Limits of integration are any a ≥ 0 to b ≤ L. L is the box length. Integrate with a given number of rectangles ( \(n_{r}\) ).
Complete this exercise in module-3-optional-exercise.jl
# include("module-3-optional-exercise-test-runner.jl") #including the file runs the tests