MATH50003 Numerical Analysis (2021–2022) Practice Computer-based Exam¶
For each problem, replace the # TODO to complete the question.
The unit tests are provided to help you test your answers.
Copyright By PowCoder代写 加微信 powcoder
You have 1 hour to complete the exam, as well as 1 hour for downloading/uploading.
Problems are marked A/B/C to indicate difficulty (“A” being most difficult).
Partial credit will be awarded for reasonable attempts even if the tests
are not passed.
You may use existing code from anywhere
but you are REQUIRED to cite the source if it is not part of the module material,
ideally by including a weblink in a comment. You MUST NOT ask for help online or
communicate with others within or outside the module.
Failure to follow these rules will be considered misconduct.
You should use the following packages:
using LinearAlgebra, SetRounding, Test
WARNING It may be necessary to restart the kernel if issues arise. Remember to reload the packages
when you do so.
1. Numbers¶
Problem 1.1 (C) Complete the following function divideby3(x) that
returns a tuple a,b such that a is the largest Float64 less
than or equal to x/3 and b is the smallest Float64 greater than or equal to x/3.
function divideby3(x)
# TODO: assign a,b so that a ≤ x ≤ b where b is either equal to a or the next float
## SOLUTION
a = setrounding(Float64, RoundDown) do
b = setrounding(Float64, RoundUp) do
divideby3 (generic function with 1 method)
x = 0.1 # arbitary x
a,b = divideby3(x)
@test a ≤ big(x)/3 ≤ b
@test b == a || b == nextfloat(a)
Test Passed
Expression: b == a || b == nextfloat(a)
2. Differentiation¶
Problem 2.1 (C) Use the following off-center finite-difference approximation
f'(x) ≈ {f(x+2h) – f(x-h) \over 3h}
with an appropriately chosen $h$ to approximate
f(x) = \cos(x^2)
at $x = 0.1$ to 5 digits accuracy.
function fd(x)
# TODO: implement a finite-difference rule
# to approximate f'(x)
# for f(x) = cos(x^2)
# with step-size h chosen to get sufficient accuracy
## SOLUTION
h = sqrt(eps())
f = x -> cos(x^2)
(f(x + 2h) – f(x – h)) / (3h)
fd (generic function with 1 method)
@test abs(fd(0.1) + 2*0.1*sin(0.1^2)) ≤ 1E-5
Test Passed
Expression: abs(fd(0.1) + 2 * 0.1 * sin(0.1 ^ 2)) ≤ 1.0e-5
Evaluated: 1.3337952660744112e-10 ≤ 1.0e-5
Problem 2.2 (A) Consider a 2D version of a dual number:
a + b ϵ_x + c ϵ_y
ϵ_x^2 = ϵ_y^2 = ϵ_x ϵ_y = 0.
Complete the following implementation supporting + and * (and
assuming a,b,c are Float64). Hint: you may need to work out on paper
how to multiply (s.a + s.b ϵ_x + s.c ϵ_y)*(t.a + t.b ϵ_x + t.c ϵ_y) using the
relationship above.
import Base: *, +, ^
struct Dual2D
a::Float64
b::Float64
c::Float64
function +(s::Dual2D, t::Dual2D)
## TODO: Implement +, returning a Dual2D
## SOLUTION
Dual2D(s.a + t.a, s.b + t.b, s.c + t.c)
function *(c::Number, s::Dual2D)
## TODO: Implement c * Dual2D(…), returning a Dual2D
## SOLUTION
Dual2D(c * s.a, c * s.b, c * s.c)
function *(s::Dual2D, t::Dual2D)
## TODO: Implement Dual2D(…) * Dual2D(…), returning a Dual2D
## SOLUTION
# we deduce (s.a + s.b ϵ_x + s.c ϵ_y)*(t.a + t.b ϵ_x + t.c ϵ_y) ==
# s.a * t.a + (s.a*t.b + s.b*t.a)*ϵ_x + (s.a*t.c + s.c*t.a)*ϵ_y
Dual2D(s.a * t.a, s.a * t.b + s.b * t.a, s.a * t.c + s.c * t.a)
* (generic function with 366 methods)
f = function (x, y) # (x+2y^2)^3 using only * and +
z = (x + 2y * y)
x,y = 1., 2.
@test f(Dual2D(x,1.,0.), Dual2D(y,0.,1.)) == Dual2D(f(x,y), 3(x+2y^2)^2, 12y*(x+2y^2)^2)
# This has computed the gradient as f(x,y) + f_x*ϵ_x + f_y*ϵ_y
# == (x+2y^2)^3 + 3(x+2y^2)^2*ϵ_x + 12y(x+2y^2)^2*ϵ_y
Test Passed
Expression: f(Dual2D(x, 1.0, 0.0), Dual2D(y, 0.0, 1.0)) == Dual2D(f(x, y), 3 * (x + 2 * y ^ 2) ^ 2, (12y) * (x + 2 * y ^ 2) ^ 2)
Evaluated: Dual2D(729.0, 243.0, 1944.0) == Dual2D(729.0, 243.0, 1944.0)
3. Structured Matrices¶
Problem 3.1 (C) Add an implementation of inv(::PermutationMatrix)
to return the inverse permutation as a PermutationMatrix. Hint: use
import Base: getindex, size, *, inv
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)) getindex(P::PermutationMatrix, k::Int, j::Int) = P.p[k] == j ? 1 : 0 *(P::PermutationMatrix, x::AbstractVector) = x[P.p] function inv(P::PermutationMatrix) # TODO: return a PermutationMatrix representing the inverse permutation ## SOLUTION PermutationMatrix(invperm(P.p)) inv (generic function with 33 methods) P = PermutationMatrix([3,4,2,1]) @test inv(P) isa PermutationMatrix @test P*inv(P) == I Test Passed Expression: P * inv(P) == I Evaluated: [1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1] == UniformScaling{Bool}(true) 4. Decompositions¶ Problem 4.1 (C) For $𝐱 ∈ ℝ^n$ implement the reflection defined by \begin{align*} 𝐲 &:= 𝐱 + \|𝐱\| 𝐞_n \\ 𝐰 &:= 𝐲/\|𝐲\| \\ Q_𝐱 &:= I - 2𝐰𝐰^⊤ \end{align*} in lowerhouseholderreflection(x), which should return a Matrix{Float64}. You may assume that x is a Vector{Float64}. function lowerhouseholderreflection(x) ## TODO: implement the householder reflector defined above ## SOLUTION y = copy(x) y[end] += norm(x) w = y/norm(y) lowerhouseholderreflection (generic function with 1 method) x = [1.0,2,3,4] Q = lowerhouseholderreflection(x) @test Q*x ≈ [zeros(3); -norm(x)] @test Q'Q ≈ I @test Q ≈ Q' Test Passed Expression: Q ≈ Q' Evaluated: [0.9807354816671587 -0.03852903666568264 -0.05779355499852396 -0.18257418583505536; -0.03852903666568264 0.9229419266686347 -0.11558710999704792 -0.3651483716701107; -0.05779355499852396 -0.11558710999704792 0.8266193350044281 -0.5477225575051661; -0.18257418583505536 -0.3651483716701107 -0.5477225575051661 -0.7302967433402217] ≈ [0.9807354816671587 -0.03852903666568264 -0.05779355499852396 -0.18257418583505536; -0.03852903666568264 0.9229419266686347 -0.11558710999704792 -0.3651483716701107; -0.05779355499852396 -0.11558710999704792 0.8266193350044281 -0.5477225575051661; -0.18257418583505536 -0.3651483716701107 -0.5477225575051661 -0.7302967433402217] Problem 4.2 (A) Complete the function ql(A) that returns a QL decomposition, that is, Q is an orthogonal matrix and L is lower triangular. You may assume that A is a square Matrix{Float64}. Hint: use Problem 4.1 to lower triangularise. function ql(A) m,n = size(A) m == n || error("not square") ## TODO Create Q and L such that Q'Q == I and L is lower triangular ## SOLUTION L = copy(A) Q = Matrix(1.0I, n, n) for j = n:-1:2 Qⱼ = lowerhouseholderreflection(L[1:j, j]) L[1:j, 1:j] = Qⱼ * L[1:j, 1:j] Q[:,1:j] = Q[:,1:j]*Qⱼ ql (generic function with 1 method) A = [1.0 2 3; 1 4 9; 1 1 1] Q,L = ql(A) @test Q'Q ≈ I @test Q*L ≈ A @test L ≈ tril(L) # it is acceptable to have small non-zero entries in L Test Passed Expression: L ≈ tril(L) Evaluated: [0.2540002540003808 -4.440892098500626e-16 5.551115123125783e-17; -1.038431997689249 -0.8254203058555567 9.992007221626409e-16; -1.362770287738494 -4.507624797904249 -9.539392014169454] ≈ [0.2540002540003808 0.0 0.0; -1.038431997689249 -0.8254203058555567 0.0; -1.362770287738494 -4.507624797904249 -9.539392014169454] 5. Singular Value Decomposition¶ Problem 5.1 (C) Find the best rank-4 approximation (in the $2$-norm) to f(x,y) = \cos(x - y)/(x+y+1) evaluated at an evenly spaced 100 × 100 grid on the square $[0,1]^2$. function bestrank4() # TODO: return best rank-4 approximation ## SOLUTION x = y = range(0, 1; length=100) F = cos.(x .- y') ./ (x .+ y' .+ 1) U,σ,V = svd(F) U[:,1:4] * Diagonal(σ[1:4]) * V[:,1:4]' Fr = bestrank4() 100×100 Matrix{Float64}: 0.999993 0.989945 0.979995 0.970141 … 0.281434 0.275776 0.270146 0.989945 0.980195 0.970537 0.960968 0.284208 0.278592 0.273005 0.979995 0.970537 0.961164 0.951874 0.286924 0.281352 0.275807 0.970141 0.960968 0.951874 0.942857 0.289584 0.284055 0.278553 0.96038 0.951486 0.942665 0.933915 0.292187 0.286703 0.281244 0.95071 0.942088 0.933534 0.925045 … 0.294735 0.289295 0.283879 0.941128 0.932772 0.924479 0.916246 0.297228 0.291832 0.286459 0.931631 0.923536 0.915499 0.907516 0.299665 0.294314 0.288985 0.922217 0.914378 0.90659 0.898853 0.302049 0.296742 0.291457 0.912884 0.905295 0.897752 0.890255 0.304378 0.299117 0.293876 0.90363 0.896286 0.888983 0.881721 … 0.306654 0.301438 0.296241 0.894453 0.887348 0.88028 0.873248 0.308877 0.303706 0.298554 0.885351 0.87848 0.871642 0.864836 0.311047 0.305921 0.300814 ⋮ ⋱ 0.333675 0.336013 0.338292 0.340512 0.347152 0.345598 0.34402 0.327754 0.330144 0.332476 0.334749 0.346235 0.344725 0.34319 0.321861 0.324303 0.326687 0.329012 … 0.345289 0.343822 0.34233 0.315998 0.318491 0.320925 0.323302 0.344314 0.342891 0.341442 0.310165 0.312707 0.315191 0.317617 0.343311 0.341931 0.340526 0.30436 0.306951 0.309483 0.311958 0.34228 0.340943 0.339581 0.298585 0.301223 0.303803 0.306325 0.341222 0.339927 0.338607 0.292839 0.295523 0.298149 0.300719 … 0.340135 0.338884 0.337606 0.287122 0.289851 0.292523 0.295138 0.339022 0.337813 0.336578 0.281434 0.284208 0.286924 0.289584 0.337882 0.336715 0.335522 0.275776 0.278592 0.281352 0.284055 0.336715 0.335591 0.334439 0.270146 0.273005 0.275807 0.278553 0.335522 0.334439 0.33333 @test rank(Fr) == 4 @test abs(Fr[10,11] - cos(x - y)/(x + y + 1)) ≤ 2E-6 Test Passed Expression: abs(Fr[10, 11] - cos(x - y) / (x + y + 1)) ≤ 2.0e-6 Evaluated: 1.317880577000885e-6 ≤ 2.0e-6 6. Differential Equations¶ Problem 6.1 (A) Complete the function airyai(n) that returns a length-$n$ Vector{Float64} \begin{bmatrix} \end{bmatrix} such that $u_k$ approximates the solution to the equation \begin{align*} u(0) &= 1 \\ u(1) &= 0 \\ u'' - x u &= 0 \end{align*} at the point $x_k = (k-1)/(n-1)$ using the second order finite-difference approximation: u''(x_k) ≈ {u_{k-1} - 2u_k + u{k+1} \over h^2} for $k = 2, …, n-1$, in $O(n)$ operations. function airy(n) # TODO: return a Vector{Float64} approximating the solution to the ODE ## SOLUTION x = range(0, 1; length = n) h = step(x) A = SymTridiagonal(fill(-2 / h^2, n - 2) .- x[2:end-1], fill(1 / h^2, n - 3)) u = A \ [-1 / h^2; zeros(n - 3)] airy (generic function with 1 method) u = airy(1000) @test u[1] == 1 @test u[end] == 0 # this compares agianst the exact formula @test abs(u[500] - 0.4757167332829094) ≤ 2E-8 Test Passed Expression: abs(u[500] - 0.4757167332829094) ≤ 2.0e-8 Evaluated: 1.8864041551047706e-8 ≤ 2.0e-8 程序代写 CS代考 加微信: powcoder QQ: 1823890830 Email: powcoder@163.com