xuebaunion@vip.163.com

3551 Trousdale Rkwy, University Park, Los Angeles, CA

留学生论文指导和课程辅导

无忧GPA：https://www.essaygpa.com

工作时间：全年无休-早上8点到凌晨3点

微信客服：xiaoxionga100

微信客服：ITCS521

python代写-MATH0058

时间：2021-05-01

Final Assessment - MATH0058

All questions must be attempted.

You must submit your solution as a single pdf file. The following options are allowed:

• Write all formulas and code together into one Jupyter Notebook and export to pdf. This may

not work on every platform where Jupyter is installed and you should test it.

• Separately write the theory part into a Word or similar document, and add the codes and

output as screenshots to it. Then convert to pdf. This is slightly more complex but should

work on most setups.

Question 1

In this question we want to explore the effects of floating point accuracy and the sensitivity of

problems with respect to small perturbations.

Q1a State the definitions of backward error, forward error, and condition number. How are these

quantities connected? Which of these are due to the sensitivity of a problem and which are due to

a specific algorithm?

Q1b We consider the root finding problem P(a, z) = 0 with P(a, z) = z2a2 + za1 + a0. Here,

a =

[

a0 a1 a2

]T is the vector of coefficients of the polynomial. Let z be a solution of P(a, z) = 0

for given coefficients a. Show that if we perturb a single coefficient aj by an amount of ∆aj to

first order the zero z is perturbed by ∆z = −z

j∆aj

P′(z) . As a slightly more advanced question try to

derive this perturbation results for polynomials of arbitrary order and not just degree 2. Based

on the perturbation result, what is the resulting condition number of a zero of P with respect to

a perturbation in the coefficient aj? (Hint: You need not give a formal proof. Consider that the

condition number is the ratio of relative forward and relative backward error.) What happens

with the condition number if the polynomial has two zeros that are close together? Comment on

the ability to accurately compute the roots of polynomials with nearby roots.

Q1c We now want to visualize the effect of small perturbations in the coefficients. Implement the

following function.

[1]: def perturbation_roots(coeffs, eps=1E-10, nrepeats=10):

"""

This function takes a polynomial z**n * coffs[0] + ... z**1 * coeffs[n-1] +

↪→coeffs[n]

and computes roots of random perturbations of this polynomial using the

np.roots function.

The perturbation is a vector of normally distributed random numbers with

2-norm eps. The function does nrepeats repetitions and at the end

returns a single one-dimensional Numpy array of type complex128

that containes all the different nrepeats * n roots of

all perturbed polynomials.

"""

1

# Delete the following line to implement your own code.

pass

Once you have implemented your function you should test it with the famous Wilkin-

son polynomial. It is defined as P(z) = ∏20j=1(z − j). Hence, its zeros are the integers

from 1 to 20. Use the perturbation_roots function to generate a single plot that shows

all the roots of this polynomial and its perturbations for the value eps=1E-10. For this

you will need the coefficients of the Wilkinson polynomial. You find them for example at

https://blogs.mathworks.com/cleve/2013/03/04/wilkinsons-polynomials/, a blog post on the

Wilkinson polynomial by Cleve Moler (the inventor of Matlab). We also give them below. The

zeros should be plotted as blue crosses. Also plot the exact zeros of the Wilkinson polynomial as

red circles.

Numerically evaluate the condition number of the root z = 15 with respect to perturbations in the

coefficients aj, separately and print the result. Interpret what you see in the graph with respect

to the computed condition numbers. To evaluate the derivative you can use the Numpy function

numpy.polyder.

[150]: # The following code gives the coefficients of the leading order monomial first,

↪→i.e. P(z) = z^{20}...

wilkinson_coeffs = [

1,

-210,

20615,

-1256850,

53327946,

-1672280820,

40171771630,

-756111184500,

11310276995381,

-135585182899530,

1307535010540395,

-10142299865511450,

63030812099294896,

-311333643161390640,

1206647803780373360,

-3599979517947607200,

8037811822645051776,

-12870931245150988800,

13803759753640704000,

-8752948036761600000,

2432902008176640000,

]

2

Question 2

Let A ∈ Rn×n be symmetric and positive definite, that is A = AT and all eigenvalues of A are

positive. In this question we want to investigate a special variant of the LU Decomposition, called

the Cholesky Decomposition, which exists for symmetric positive definite matrices.

The Cholesky Decomposition is defined as

A = CTC

for C a nonsingular upper triangular matrix with positive entries on its diagonal.

Q2a Show that the Cholesky Decomposition exists if and only if A is positive definite. In order

to show that the Cholesky Decomposition exists if A is positive definite proceed by induction.

Assume that the Cholesky Decomposition exists for matrices of dimenion n − 1 and show that

you can construct a Cholesky Decomposition for a matrix of dimension n of the form

[

Aˆ b

bT α

]

=

[

CˆT 0

cT δ

] [

Cˆ c

0 δ

]

with δ > 0, where Aˆ = CˆTCˆ is a Cholesky Decomposition of dimension n− 1. Be careful to also

mention the induction start. In order to show that δ > 0 you may find the following formula for

the determinant of a block matrix helpful:

det

([

A B

C D

])

= det A det(D− CA−1B).

Q2b Use the Singular Value Decomposition to show that ‖A‖2 = ‖C‖22.

Q2c From the lectures you know that the LU Decomposition of a matrix can have a large backward

error if the growth factor ρ is large. Based on this, give a heuristical argument why you expect the

Cholesky Decomposition to be backward stable.

Q2d The decomposition in a) suggests an algorithm for computing the Cholesky Decomposition

by starting at the upper left element a11 of A and successively computing the Cholesky Decompo-

sition for larger and larger submatrices until you have computed the Cholesky Decomposition of

the whole matrix.

Implement a function def cholesky(A) that for a given positive definite matrix A implements

this algorithm. The skeleton code below will help you with the implementation and testing. You

want to achieve that the printed residual is a small multiple of machine precision.

[ ]: import numpy

def cholesky(A):

"""

Compute the Cholesky decomposition of A.

Return the upper triangular matrix C of

the Cholesky decomposition A = C^TC

3

of a symmetric positive definite matrix A.

"""

# Delete the following line for your implementation.

pass

n = 50

B = numpy.random.randn(n, n)

A = B.T @ B # Computes a positive definite random matrix.

C = cholesky(A)

relative_residual = numpy.linalg.norm(C.T @ C - A) / numpy.linalg.norm(A)

print(f"The relative residual is: {relative_residual}")

Question 3

Consider the following algorithm:

1.) Start with random x0 ∈ Rn, ‖x0‖2 = 1

2.) Until convergence do the following steps

• Compute y = Axj

• Compute z = ATy

• Compute γj = xTj z

• Compute xj+1 = z/‖z‖2.

Q3a Use the Singular Value Decomposition of A to investigate the convergence of this algorithm.

In particular, show under which conditions the algorithm converges and what xj and γj converge

to.

Q3b Implement the algorithm. Use the following skeleton for your implementation.

[ ]:

[ ]: def fun(A, nsteps=100):

"""

Perform nsteps steps of the algorithm.

Return a vector containing the values

\gamma_1, \gamma_2, \dots

"""

# Delete the following line for your implementation.

pass

Now plot the convergence for the algorithm in the following two cases

4

• A random 50× 50 matrix wwhose entries are uniformly distributed random numbers.

• A random 50× 50 matrix whose entries are normally distributed random numbers.

For the convergence plot the values |γj − γ|/|γ| using a semilogarithmic plot. Here γ is the theo-

retical limit of the iteration as j → ∞. If you don’t know what this is you can use instead the last

value γend that was returned in the output vector from fun.

Q3c Explain the observed convergence behaviour as accurately as possible.

Question 4

In this question we want to investigate quadratic eigenvalue problems of the form

T(λ)x = 0

with T(λ) being a complex function T : C→ Cn×n defined by

T(λ) = A2λ2 + A1λ+ A0,

where A0, A1, A2 are complex n× n matrices.

We say that λj is an eigenvalue of T if there exists an associated vector xj 6= 0 such that T(λj)xj = 0.

Q4a Show that the quadratic eigenvalue problem T(λ)x = 0 is identical to the following gener-

alised linear eigenvalue problem of dimension 2n.

[

0 I

−A0 −A1

] [

x

u

]

− λ

[

I 0

0 A2

] [

x

u

]

= 0.

(Note that a generalized linear eigenvalue problem is an eigenvalue problem of the form Ax =

λBx. Hence, a standard eigenvalue problem is obtained for B being the identity matrix. The Scipy

eig and eigvals functions can solve generalized linear eigenvalue problems.)

Q4b In all the following parts we investigate a specific quadratic eigenvalue problem defined by

the matrices in the following Python code.

[131]: A0 = numpy.array(

[[121, 18.9, 15.9],

[0, 2.7, 0.145],

[11.9, 3.64, 15.5],

]

)

A1 = numpy.array(

[[7.66, 2.45, 2.1],

[0.23, 1.04, 0.223],

[0.6, 0.756, 0.658],

]

)

5

A2 = numpy.array(

[[17.6, 1.28, 2.89],

[1.28, 0.824, 0.413],

[2.89, 0.413, 0.725]

]

)

Numerically compute these eigenvalues and plot them in the complex plane.

Q4c One way to interpet the numerical stability of computed eigenvalues for quadratic eigenvalue

problems is to consider the resolvent function f : z → T(z)−1 in the complex plane. In order to

visualize this function evaluate the function ‖ f (z)‖2 for complex z = x+ iy with x ∈ [−2, 2] and

y ∈ [−10, 10]. Use a plot grid of 100 points in the x-direction and 200 points in the y-direction.

Then plot log10 ‖ f (z)‖2 (the log allows you to distinguish magnitudes in the values). You can cre-

ate the plot grid with numpy.meshgrid, and plot the result with matplotlib.pyplot.imshow. Also

display a colorbar (matplotlib.pyplot.colorbar). Explain heuristically why your numerical re-

sults suggest that the eigenvalues are stable with respect to perturbations.

Q4d Finally, we want to investigate a simple iteration to solve the quadratic eigenvalue problem

that can also be generalised to other types of matrix functions T(λ) beyond the quadratic case.

The following method goes back to Ruhe.

Instead of the problem T(λ)x = 0 we perform a Taylor expansion around some approximate

eigenvalue λ(j) and solve the linear problem, resulting in

T(λ)x ≈

(

T(λ(j))− ηT′(λ(j))

)

x = 0

with λ = λ(j) − η. In order to compute η we find the smallest magnitude eigenvalue of the

generalized eigenvalue problem

T(λ(j))x = ηT′(λ(j))x

and update λ(j+1) = λ(j) − η.

Implement this iteration in Python using the following skeleton.

[155]: def ruhe_method(A0, A1, A2, lam0, nsteps=100):

"""

Implement the Ruhe method to find a single eigenvalue of a quadratic

↪→eigenvalue problem.

A0, A1, A2 are the coefficient matrices for the quadratic eigenvalue problem.

↪→ The value

lam0 is the initial eigenvalue approximation and steps is the number of

↪→iteration steps.

6

The method returns a vector containing the eigenvalue approximation

↪→\lambda^{(j)} in each

iteration step.

"""

# Delete the following line for your implementation.

pass

Now pick an exact eigenvalue λ of the given quadratic eigenvalue problem and as starting value

for the iteration choose λ(0) = λ + .1. Run Ruhe’s method and make a semilog plot of the con-

vergence of the relative error |λ(j) − λ|/|λ|. By inspecting the values of the relative errors make a

conjecture about the rate of convergence of Ruhe’s method.

7

学霸联盟

All questions must be attempted.

You must submit your solution as a single pdf file. The following options are allowed:

• Write all formulas and code together into one Jupyter Notebook and export to pdf. This may

not work on every platform where Jupyter is installed and you should test it.

• Separately write the theory part into a Word or similar document, and add the codes and

output as screenshots to it. Then convert to pdf. This is slightly more complex but should

work on most setups.

Question 1

In this question we want to explore the effects of floating point accuracy and the sensitivity of

problems with respect to small perturbations.

Q1a State the definitions of backward error, forward error, and condition number. How are these

quantities connected? Which of these are due to the sensitivity of a problem and which are due to

a specific algorithm?

Q1b We consider the root finding problem P(a, z) = 0 with P(a, z) = z2a2 + za1 + a0. Here,

a =

[

a0 a1 a2

]T is the vector of coefficients of the polynomial. Let z be a solution of P(a, z) = 0

for given coefficients a. Show that if we perturb a single coefficient aj by an amount of ∆aj to

first order the zero z is perturbed by ∆z = −z

j∆aj

P′(z) . As a slightly more advanced question try to

derive this perturbation results for polynomials of arbitrary order and not just degree 2. Based

on the perturbation result, what is the resulting condition number of a zero of P with respect to

a perturbation in the coefficient aj? (Hint: You need not give a formal proof. Consider that the

condition number is the ratio of relative forward and relative backward error.) What happens

with the condition number if the polynomial has two zeros that are close together? Comment on

the ability to accurately compute the roots of polynomials with nearby roots.

Q1c We now want to visualize the effect of small perturbations in the coefficients. Implement the

following function.

[1]: def perturbation_roots(coeffs, eps=1E-10, nrepeats=10):

"""

This function takes a polynomial z**n * coffs[0] + ... z**1 * coeffs[n-1] +

↪→coeffs[n]

and computes roots of random perturbations of this polynomial using the

np.roots function.

The perturbation is a vector of normally distributed random numbers with

2-norm eps. The function does nrepeats repetitions and at the end

returns a single one-dimensional Numpy array of type complex128

that containes all the different nrepeats * n roots of

all perturbed polynomials.

"""

1

# Delete the following line to implement your own code.

pass

Once you have implemented your function you should test it with the famous Wilkin-

son polynomial. It is defined as P(z) = ∏20j=1(z − j). Hence, its zeros are the integers

from 1 to 20. Use the perturbation_roots function to generate a single plot that shows

all the roots of this polynomial and its perturbations for the value eps=1E-10. For this

you will need the coefficients of the Wilkinson polynomial. You find them for example at

https://blogs.mathworks.com/cleve/2013/03/04/wilkinsons-polynomials/, a blog post on the

Wilkinson polynomial by Cleve Moler (the inventor of Matlab). We also give them below. The

zeros should be plotted as blue crosses. Also plot the exact zeros of the Wilkinson polynomial as

red circles.

Numerically evaluate the condition number of the root z = 15 with respect to perturbations in the

coefficients aj, separately and print the result. Interpret what you see in the graph with respect

to the computed condition numbers. To evaluate the derivative you can use the Numpy function

numpy.polyder.

[150]: # The following code gives the coefficients of the leading order monomial first,

↪→i.e. P(z) = z^{20}...

wilkinson_coeffs = [

1,

-210,

20615,

-1256850,

53327946,

-1672280820,

40171771630,

-756111184500,

11310276995381,

-135585182899530,

1307535010540395,

-10142299865511450,

63030812099294896,

-311333643161390640,

1206647803780373360,

-3599979517947607200,

8037811822645051776,

-12870931245150988800,

13803759753640704000,

-8752948036761600000,

2432902008176640000,

]

2

Question 2

Let A ∈ Rn×n be symmetric and positive definite, that is A = AT and all eigenvalues of A are

positive. In this question we want to investigate a special variant of the LU Decomposition, called

the Cholesky Decomposition, which exists for symmetric positive definite matrices.

The Cholesky Decomposition is defined as

A = CTC

for C a nonsingular upper triangular matrix with positive entries on its diagonal.

Q2a Show that the Cholesky Decomposition exists if and only if A is positive definite. In order

to show that the Cholesky Decomposition exists if A is positive definite proceed by induction.

Assume that the Cholesky Decomposition exists for matrices of dimenion n − 1 and show that

you can construct a Cholesky Decomposition for a matrix of dimension n of the form

[

Aˆ b

bT α

]

=

[

CˆT 0

cT δ

] [

Cˆ c

0 δ

]

with δ > 0, where Aˆ = CˆTCˆ is a Cholesky Decomposition of dimension n− 1. Be careful to also

mention the induction start. In order to show that δ > 0 you may find the following formula for

the determinant of a block matrix helpful:

det

([

A B

C D

])

= det A det(D− CA−1B).

Q2b Use the Singular Value Decomposition to show that ‖A‖2 = ‖C‖22.

Q2c From the lectures you know that the LU Decomposition of a matrix can have a large backward

error if the growth factor ρ is large. Based on this, give a heuristical argument why you expect the

Cholesky Decomposition to be backward stable.

Q2d The decomposition in a) suggests an algorithm for computing the Cholesky Decomposition

by starting at the upper left element a11 of A and successively computing the Cholesky Decompo-

sition for larger and larger submatrices until you have computed the Cholesky Decomposition of

the whole matrix.

Implement a function def cholesky(A) that for a given positive definite matrix A implements

this algorithm. The skeleton code below will help you with the implementation and testing. You

want to achieve that the printed residual is a small multiple of machine precision.

[ ]: import numpy

def cholesky(A):

"""

Compute the Cholesky decomposition of A.

Return the upper triangular matrix C of

the Cholesky decomposition A = C^TC

3

of a symmetric positive definite matrix A.

"""

# Delete the following line for your implementation.

pass

n = 50

B = numpy.random.randn(n, n)

A = B.T @ B # Computes a positive definite random matrix.

C = cholesky(A)

relative_residual = numpy.linalg.norm(C.T @ C - A) / numpy.linalg.norm(A)

print(f"The relative residual is: {relative_residual}")

Question 3

Consider the following algorithm:

1.) Start with random x0 ∈ Rn, ‖x0‖2 = 1

2.) Until convergence do the following steps

• Compute y = Axj

• Compute z = ATy

• Compute γj = xTj z

• Compute xj+1 = z/‖z‖2.

Q3a Use the Singular Value Decomposition of A to investigate the convergence of this algorithm.

In particular, show under which conditions the algorithm converges and what xj and γj converge

to.

Q3b Implement the algorithm. Use the following skeleton for your implementation.

[ ]:

[ ]: def fun(A, nsteps=100):

"""

Perform nsteps steps of the algorithm.

Return a vector containing the values

\gamma_1, \gamma_2, \dots

"""

# Delete the following line for your implementation.

pass

Now plot the convergence for the algorithm in the following two cases

4

• A random 50× 50 matrix wwhose entries are uniformly distributed random numbers.

• A random 50× 50 matrix whose entries are normally distributed random numbers.

For the convergence plot the values |γj − γ|/|γ| using a semilogarithmic plot. Here γ is the theo-

retical limit of the iteration as j → ∞. If you don’t know what this is you can use instead the last

value γend that was returned in the output vector from fun.

Q3c Explain the observed convergence behaviour as accurately as possible.

Question 4

In this question we want to investigate quadratic eigenvalue problems of the form

T(λ)x = 0

with T(λ) being a complex function T : C→ Cn×n defined by

T(λ) = A2λ2 + A1λ+ A0,

where A0, A1, A2 are complex n× n matrices.

We say that λj is an eigenvalue of T if there exists an associated vector xj 6= 0 such that T(λj)xj = 0.

Q4a Show that the quadratic eigenvalue problem T(λ)x = 0 is identical to the following gener-

alised linear eigenvalue problem of dimension 2n.

[

0 I

−A0 −A1

] [

x

u

]

− λ

[

I 0

0 A2

] [

x

u

]

= 0.

(Note that a generalized linear eigenvalue problem is an eigenvalue problem of the form Ax =

λBx. Hence, a standard eigenvalue problem is obtained for B being the identity matrix. The Scipy

eig and eigvals functions can solve generalized linear eigenvalue problems.)

Q4b In all the following parts we investigate a specific quadratic eigenvalue problem defined by

the matrices in the following Python code.

[131]: A0 = numpy.array(

[[121, 18.9, 15.9],

[0, 2.7, 0.145],

[11.9, 3.64, 15.5],

]

)

A1 = numpy.array(

[[7.66, 2.45, 2.1],

[0.23, 1.04, 0.223],

[0.6, 0.756, 0.658],

]

)

5

A2 = numpy.array(

[[17.6, 1.28, 2.89],

[1.28, 0.824, 0.413],

[2.89, 0.413, 0.725]

]

)

Numerically compute these eigenvalues and plot them in the complex plane.

Q4c One way to interpet the numerical stability of computed eigenvalues for quadratic eigenvalue

problems is to consider the resolvent function f : z → T(z)−1 in the complex plane. In order to

visualize this function evaluate the function ‖ f (z)‖2 for complex z = x+ iy with x ∈ [−2, 2] and

y ∈ [−10, 10]. Use a plot grid of 100 points in the x-direction and 200 points in the y-direction.

Then plot log10 ‖ f (z)‖2 (the log allows you to distinguish magnitudes in the values). You can cre-

ate the plot grid with numpy.meshgrid, and plot the result with matplotlib.pyplot.imshow. Also

display a colorbar (matplotlib.pyplot.colorbar). Explain heuristically why your numerical re-

sults suggest that the eigenvalues are stable with respect to perturbations.

Q4d Finally, we want to investigate a simple iteration to solve the quadratic eigenvalue problem

that can also be generalised to other types of matrix functions T(λ) beyond the quadratic case.

The following method goes back to Ruhe.

Instead of the problem T(λ)x = 0 we perform a Taylor expansion around some approximate

eigenvalue λ(j) and solve the linear problem, resulting in

T(λ)x ≈

(

T(λ(j))− ηT′(λ(j))

)

x = 0

with λ = λ(j) − η. In order to compute η we find the smallest magnitude eigenvalue of the

generalized eigenvalue problem

T(λ(j))x = ηT′(λ(j))x

and update λ(j+1) = λ(j) − η.

Implement this iteration in Python using the following skeleton.

[155]: def ruhe_method(A0, A1, A2, lam0, nsteps=100):

"""

Implement the Ruhe method to find a single eigenvalue of a quadratic

↪→eigenvalue problem.

A0, A1, A2 are the coefficient matrices for the quadratic eigenvalue problem.

↪→ The value

lam0 is the initial eigenvalue approximation and steps is the number of

↪→iteration steps.

6

The method returns a vector containing the eigenvalue approximation

↪→\lambda^{(j)} in each

iteration step.

"""

# Delete the following line for your implementation.

pass

Now pick an exact eigenvalue λ of the given quadratic eigenvalue problem and as starting value

for the iteration choose λ(0) = λ + .1. Run Ruhe’s method and make a semilog plot of the con-

vergence of the relative error |λ(j) − λ|/|λ|. By inspecting the values of the relative errors make a

conjecture about the rate of convergence of Ruhe’s method.

7

学霸联盟