MATH50003 Numerical Analysis: Problem 3¶
This problem sheet explores implementation of triangular solves,
supporting a matrix with two super-diagonals, as well as
Copyright By PowCoder代写 加微信 powcoder
permutation and Householder reflections that can be applied to a vector in
$O(n)$ complexity.
using LinearAlgebra, Test
# We will override these functions below
import Base: getindex, setindex!, size, *, \
1. Dense Matrices¶
Problem 1.1 Show that A*x is not
implemented as mul(A, x) from the lecture notes
by finding a Float64 example where the bits do not match.
2. Triangular Matrices¶
Problem 2.1 Complete the following functions for solving linear systems with
triangular systems by implementing back and forward-substitution:
function ldiv(U::UpperTriangular, b)
n = size(U,1)
if length(b) != n
error(“The system is not compatible”)
x = zeros(n) # the solution vector
## TODO: populate x using back-substitution
function ldiv(U::LowerTriangular, b)
n = size(U,1)
if length(b) != n
error(“The system is not compatible”)
x = zeros(n) # the solution vector
## TODO: populate x using forward-substitution
Problem 2.2⋆ Given $𝐱 ∈ ℝ^n$, find a lower triangular matrix of the form
L = I – 2 𝐯 𝐞_1^⊤
such that:
L 𝐱 = x_1 𝐞_1.
What does $L𝐲$ equal if $𝐲 ∈ ℝ^n$ satisfies $y_1 = 𝐞_1^⊤ 𝐲 = 0$?
3. Banded matrices¶
Problem 3.1 Complete the implementation of UpperTridiagonal which represents a banded matrix with
bandwidths $(l,u) = (0,2)$:
struct UpperTridiagonal{T} <: AbstractMatrix{T}
d::Vector{T} # diagonal entries
du::Vector{T} # super-diagonal enries
du2::Vector{T} # second-super-diagonal entries
size(U::UpperTridiagonal) = (length(U.d),length(U.d))
function getindex(U::UpperTridiagonal, k::Int, j::Int)
d,du,du2 = U.d,U.du,U.du2
# TODO: return U[k,j]
function setindex!(U::UpperTridiagonal, v, k::Int, j::Int)
d,du,du2 = U.d,U.du,U.du2
if j > k+2
error(“Cannot modify off-band”)
# TODO: modify d,du,du2 so that U[k,j] == v
U # by convention we return the matrix
Problem 3.2 Complete the following implementations of * and \ for UpperTridiagonal so that
they take only $O(n)$ operations.
function *(U::UpperTridiagonal, x::AbstractVector)
T = promote_type(eltype(U), eltype(x)) # make a type that contains both the element type of U and x
b = zeros(T, size(U,1)) # returned vector
# TODO: populate b so that U*x == b (up to rounding)
function \(U::UpperTridiagonal, b::AbstractVector)
T = promote_type(eltype(U), eltype(b)) # make a type that contains both the element type of U and b
x = zeros(T, size(U,2)) # returned vector
# TODO: populate x so that U*x == b (up to rounding)
\ (generic function with 136 methods)
4. Permutations¶
Problem 4.1⋆ What are the permutation matrices corresponding to the following permutations?
\begin{pmatrix}
1 & 2 & 3 \\
\end{pmatrix}, \begin{pmatrix}
1 & 2 & 3 & 4 & 5 & 6\\
2 & 1 & 4 & 3 & 6 & 5
\end{pmatrix}.
Problem 4.2 Complete the implementation of a type representing
permutation matrices that supports P[k,j] and such that * takes $O(n)$ operations.
struct PermutationMatrix <: AbstractMatrix{Int} p::Vector{Int} # represents the permutation whose action is v[p] function PermutationMatrix(p::Vector) sort(p) == 1:length(p) || error("input is not a valid permutation") size(P::PermutationMatrix) = (length(P.p),length(P.p)) function getindex(P::PermutationMatrix, k::Int, j::Int) # TODO: Return P[k,j] function *(P::PermutationMatrix, x::AbstractVector) # TODO: permute the entries of x # If your code is correct, this "unit test" will succeed p = [1, 4, 2, 5, 3] P = PermutationMatrix(p) @test P == I(5)[p,:] 5. Orthogonal matrices¶ Problem 5.1⋆ Show that orthogonal matrices preserve the 2-norm of vectors: \|Q 𝐯\| = \|𝐯\|. Problem 5.2⋆ Show that the eigenvalues $λ$ of an orthogonal matrix $Q$ are on the unit circle: $|λ| = 1$. Problem 5.3⋆ Explain why an orthogonal matrix $Q$ must be equal to $I$ if all its eigenvalues are 1. Problem 5.4 Complete the implementation of a type representing reflections that supports Q[k,j] and such that * takes $O(n)$ operations. # Represents I - 2v*v' struct Reflection{T} <: AbstractMatrix{T} v::Vector{T} Reflection(x::Vector{T}) where T = Reflection{T}(x/norm(x)) size(Q::Reflection) = (length(Q.v),length(Q.v)) function getindex(Q::Reflection, k::Int, j::Int) # TODO: Return Q[k,j] function *(Q::Reflection, x::AbstractVector) # TODO: permute the entries of x # If your code is correct, these "unit tests" will succeed x = randn(5) Q = Reflection(x) v = x/norm(x) @test Q == I-2v*v' @test Q*v ≈ -v @test Q'Q ≈ I Problem 5.5 Complete the following implementation of a Housholder reflection so that the unit tests pass. Here s == true means the Householder reflection is sent to the positive axis and s == false is the negative axis. function householderreflection(s::Bool, x::AbstractVector) # TODO: implement Householder reflection, returning a Reflection x = randn(5) Q = householderreflection(true, x) @test Q isa Reflection @test Q*x ≈ norm(x)*[1; zeros(length(x)-1)] Q̃ = householderreflection(false, x) @test Q̃ isa Reflection @test Q̃*x ≈ -norm(x)*[1; zeros(length(x)-1)] 程序代写 CS代考 加微信: powcoder QQ: 1823890830 Email: powcoder@163.com