SciComp Project 1 Solved

SciComp Project 1

1        Background
 Many of you will know the phenomenon that a prism refracts light, i.e. splits it up in different colors, because the refractive index of the prism varies with the frequency ω of light. The underlying property of the molecules forming the material is the frequency dependent polarizability, α(ω). The polarizability, like many other properties of molecules and materials, can be calculated from the basic laws of physics, here the the timedependent Schr¨odinger equation.

In the end, the polarizability for a given frequency ω of the incoming light is obtained as the following scalar product of two column vectors z and x:

α(ω) = zTx                                                    (1)
Figure 1: A dispersive prism slows light at different rates depending on the wave-length, causing refraction.
where z is a vector that can be calculated from the Schr¨odinger equation, and x is the solution to the following system of linear equations:

                                                              (E − ωS)x = z                                                                                       (2)

Here, E and S are two square matrices, and ω is the frequency of the incoming light. Like the column vector z, the matrices E and S are calculated from the Schr¨odinger equation, which we will not discuss further in this course.

2        Data
In this project we consider the water molecule, H2O, and its frequency dependent polarizability, α(ω). It turns out that the matrices E and S, and the column vector z, have a structure that lets us decompose the matrices into submatrices as follows:

A B

                                   E =                       ,                             S,                                                                                            z(3)

B A

The Python-file watermatrices.py and the Matlab-file watermatrices.mat contain the submatrices A and B, and the subvector y for a water molecule, obtained by an approximate solution to the Schr¨odinger equation.

3        Questions for Week 1
a.    (1) Write a small function that computes the condition number of a matrix under the max-norm:

cond 

Use a library matrix inversion routine[1] for M−1, but do program the max-norm yourself using sum, abs, and max. (2) For three frequencies, ω = {0.800,1.146,1.400}, calculate the condition number for the matrix E − ωS. The right-hand-side z is given with 8 significant digits. How many digits in the solution x could we trust if everything else were exact? Why?

b.   (1) For each of the three ω, determine a bound on the relative forward error in the max-norm:

 cond∞ 

for the perturbation that the frequency ω is changed by is given with

3 digits after the comma, how many digits in x could we be guarantee are valid if everything else were exact? Why?

c.    Implement three separate functions

1.   L,U = lu factorize(A), which takes a square matrix M as input and returns two square matrices: A triangular matrix L and upper triangular matrix U such that

M = LU.

2.   y = forward substitute(L,z), which takes a square lower triangular matrix L and a vector b as input, and returns the solution vector y to Ly = b.

3.   x = back substitute(U,y), which takes a square upper triangular matrix U and a vector y as input, and returns the solution vector x to Ux = y.

and test them with the linear equation

 

You can test your solution against a library routine.[2]

If you know how, try to use vector operations instead of for-loops where possible (orders of magnitude faster in Python and Matlab).

d.   Implement a function α = solve alpha(omega) for calculating the frequency-dependent polarizability α(ω) = zTx for water in the given approximation. This routine should solve Equation (2) by LU-factorization using your own three routines from (c). (1) Using your routine, make a table of the polarizabilities for the frequencies given in (a) and their perturbations, i.e. for ω = {0.800 ± δω,1.146 ± δω,1.400 ± δω}, with  as before. (2)

Which error-bound is the correct one to understand the variation of the calculated polarizabilities due to the perturbation: (a) or (b) or both? Do your calculated values fall within the bounds you calculated above?

e.    (1) Compute a table of α(ω) for 1000 evenly spaced values in the interval [0.7,1.5][3] using your routine from (d), and plot the values. (2) What happens to the linear system of Equation

(2) around the frequency ω = 1.146307999, and how is this reflected in α(ω)?

4        Questions for Week 2
f.     Implement two routines and test them with the matrix Atest and right-hand-side b:

Atest   and btest  

1.   Q,R = householder QR(A), which takes as input a rectangular matrix A: m × n and uses the Householder method to compute its QR decomposition. Check that Q: m×m is orthogonal, i.e., QTQ = QQT = I, and that the upper triangular matrix R: m × n satisfies A = QR.

2.   A more efficient version of (f.1) that doesn’t compute the m×m matrix Q, but just stores the reflection vectors v. You can either give it the interface V,R = householder fast(A) or (as is done in practice) let the result be a combined (m+1)×n matrix VR, in which the upper n × n triangle contains R, and the kth column below the diagonal is vk (of length m − k).

3.   ˜x = least squares(A,b), which combines this routine with your back-substitution from (c.3) to compute a linear least squares fitting. It should take as input a rectangular m×n matrix A and an m×1 right-hand-side vector b, returning an n×1 approximate solution vector x˜ to Ax˜ ' b as output. If you did not get (f.2) to work, you can just use your solution to (f.1).

g.    We now want to approximate α(ω) in the interval [0.7,ωp] using a polynomial

                                                                                                                                            (4)

1.   Suggest a suitable value of ωp < 1.5. What makes this choice reasonable? (2) Find the coefficients of the polynomial using your linear least squares routine from (f)[4] , the table computed above, and n = 4. (3) Repeat the computation for n = 6 and compare the accuracies of the two polynomials: Plot the magnitude of the relative error (in a log10-scale) of the polynomial approximation as the difference between the P(ω) and α(ω) values. (4) How many correct digits does each approximation yield?

h.   We would now like to approximate α(ω) in the entire interval [0.7,1.5], and choose the rational approximating function, which is able to represent singularities:

                                                                                                                                 (5)

1.   Why will this fail with the polynomial approximation of Problem (g)? And why can Eq. (5) do a better job? (2) Find the coefficients aj and bj using your linear least squares routine, the table of α-values computed above, and n = 2. You need to reformulate the expression as an linear approximation so that you can use a linear least squares fitting. Plot the error of the the rational-function approximation Q(ω) compared to α(ω) calculated by Equations (1) and (2).[5] (3) Repeat the computation for n = 4 and compare the accuracies of the two approximations quantitatively.


 
[1] E.g. inv from numpy.linalg for Python, inv for Matlab. But only here.
[2] E.g. numpy.linalg.solve in Python, or linsolve in Matlab.
[3] Use linspace in both Python and Matlab.
[4] If you could not get your own least squares solver to work, you can use a library version for the remaining problms - but make sure to write clearly that you have done so.
[5] Once you have computed the coefficients a and b, be sure to finish the construction of Q(ω) using Eq. (5), and to use this for the error. It is tempting (but wrong) to just use the linear approximation you used to find a and b, but that is linear and hence cannot represent singularities.