Examination projects

1. Yet another cubic sub-spline

   Consider a data set {x_i, y_i}_i=1,..,n which represents a tabulated function.

   Implement a sub-spline of this data using the following algorithm:

   • for each inner point x_i, i=2,..,n-1, build a quadratic interpolating polynomial through the points x_i-1, x_i, x_i+1 and estimate, using this polynomial, the first derivative, p_i, of the function at the point x_i;
   • estimate p₁ and p_n from the the same polynomial that you used to estimate correspondingly p₂ and p_n-1.
   • now you have a data set {x_i, y_i, p_i}_i=1,..,n where the function is tabulated together with its derivative: build a cubic sub-spline for it, S_i(x)= y_i +b_i(x-x_i) +c_i(x-x_i)² +d_i(x-x_i)³,
     where for each interval the three coefficients b_i, c_i, d_i are determined by the three conditions, S_i(x_i+1)=y_i+1,
     S'_i(x_i)=p_i,
     S'_i(x_i+1)=p_i+1.

   See the subsection "Akima sub-spline interpolation" for the inspiration

   Check whether it makes smaller wiggles than quadratic and cubic splines.

2. Cubic sub-spline for data with derivatives

   For a data set, {x_i, y_i, y'_i}_i=1,..,n , where the function is tabulated together with its derivative — for example, as a result of a numerical integration of a second order ordinary differential equation — build a cubic sub-spline, S_i(x)= y_i +b_i(x-x_i) +c_i(x-x_i)² +d_i(x-x_i)³, where for each interval the three coefficients b_i, c_i, d_i are determined by the three conditions, S_i(x_i+1)=y_i+1,
   S'_i(x_i)=y'_i,
   S'_i(x_i+1)=y'_i+1.

   See the subsection "Akima sub-spline interpolation" for the inspiration.

   A possible extension: make the second derivative of the spline continuous by increasing the order of the spline to four, for example, in the form

   S_i(x)= y_i +b_i(x-x_i) +c_i(x-x_i)² +d_i(x-x_i)³ +e_i(x-x_i)²(x-x_i+1)² ,

   and choosing the coefficients e_i such that the spline has continuous second derivative.

3. Akima sub-spline

   Implement the Akima sub-spline. See the section "Akima sub-spline interpolation" in the book for the inspiration.

4. Bi-linear interpolation on a rectilinear grid in two dimensions

   A rectilinear grid (note that rectilinear is not necessarily cartesian nor regular) in two dimensions is a set of n_x×n_y points where each point can be adressed by a double index (i,j) where 1≤i≤n_x, 1≤j≤n_y and the coordinates of the point (i,j) are given as (x_i,y_i), where x and y are vectors with sizes n_x and n_y correspondingly. The values of the tabulated function F at the grid points can then be arranged as a matrix {F_i,j=F(x_i,y_i)}.

   Build an interpolating routine which takes as the input the vectors x_i and y_i, and the matrix F_i,j and returns the bi-linear interpolated value of the function at a given 2D-point p=(p_x,p_y).

   See the chapter "Bi-linear interpolation" in the book.

   The signature of the interpolating subroutine can be

   double bilinear(int nx, int ny, double* x, double* y, double** F, double px, double py)

   or

   double bilinear(gsl_vector* x, gsl_vector* y, gsl_matrix* F, double px, double py)

5. Generalized eigenvalue problem

   A generalized eigenvalue problem is the problem of finding all generalized eigenvectors x and the corresponding generalized eigenvalues λ which obey Ax=λNx, where A is a real symmetric matrix and N is a real symmetric positive-defined matrix (all eigenvalues are positive).

   Implement a routine for solving the generalized eigenvalue problem using the following algorithm:

   1. Using your Jacobi-eigenvalue routine find the eigenvalue decomposition of the matrix N, N=VDV^T, where D is the diagonal matrix with (positive) eigenvalues of the matrix N on the diagonal and V is the matrix with the corresponding eigenvectors as its columns. Now the original generalized eigenvalue problem can be represented as an ordinary eigenvalue problem (prove it analytically), By=λy, where the real symmetric matrix B=√D^-1V^TAV√D^-1, and y=√DV^Tx.
   2. Using your Jacobi-eigenvalue routine find the the eigenvalues and eigenvectors of the matrix B.
   3. Restore the original eigenvectors x.

6. Jacobi eigenvalue algorithm: several smallest (largest) eigenvalues.

   Implement the variant of the Jacobi eigenvalue algorithm which calculates the given number n of the smallest (largest) eigenvalues and corresponding eigenvectors by consecutively zeroing the off-diagonal elements in the first (last) n rows of a real symmetrix matrix.

   Hints: see exercise "matrix diagonalization" section B.

7. Tridiagonalization of a symmetric matrix: Jacobi (Givens) rotations

   Consider a Jacobi rotation with the matrix J(p,q),

   A→J(p,q)^TAJ(p,q)

   where the rotation angle is chosen not to zero the element A_p,q but rather to zero the element A_p-1,q. Argue that the following sequence of such Jacobi rotations,

   J(2,3), J(2,4), ..., J(2,n); J(3,4), ..., J(3,n); ...; J(n-1,n)

   reduces a symmetric matrix to a tridiagonal form. Implement this algorithm. Remember to accumulate the total transformation matrix.

8. Jacobi diagonalization with classic sweeps and indexing

   Implement the classic Jacobi diagonalization algorithm:

   • do:
     • find the largest off-diagonal element;
     • zero it using the Jacobi rotation;
     until converged;

   A possible improvement -- indexing:

   • Before starting diagonalization find out the largest off-diagonal elements in each row and store their indices in a separate index array. It's O(n²) operations—expensive—but you only do it only once in the beginning.
   • Find the largest off-diagonal element from the pre-found largest elements in each row (only O(n) operations) and zero it by a Jacobi rotation.
   • Update your index. Since the Jacobi rotation only affects elements in rows and columns with indices p and q the update is also only O(n) operations.
   • Repeat the last two steps until converged

9. Inverse iteration algorithm for eigenvalues (and eigenvectors)

   Well, just implement it.

   See the chapter "Power iteration methids and Krylov subspaces" in the book.

10. Lanczos tridiagonalization algorithm

    Implement the Lanczos algorithm (Lanczos interation) for real symmetric matrices.

    Possible extention: Arnoldi iteration and GMRES.

11. Symmetric rank-1 update of a size-n symmetric eigenvalue problem

    The matrix A to diagonalize is given in the form

    A = D +uu^T,

    where D is a diagonal matrix and u is a column-vector.

    Given the diagonal elements of the matrix D and the elements of the vector u find the eigenvalues of the matrix A using only O(n²) operations (see section "Eigenvalues of updated matrix" in the book).

12. Symmetric row/column update of a size-n symmetric eigenvalue problem

    The matrix to diagonalize is given in the form

    A = D + e(p) u^T + u e(p)^T

    where D is a diagonal matrix with diagonal elements {d_k, k=1,...,n}, u is a given column-vector, and the vector e(p) with components e(p)_i=δ_ip is a unit vector in the direction p where 1≤p≤n.

    Given the diagonal elements of the matrix D, the vector u, and the integer p, calculate the eigenvalues of the matrix A using O(n²) operations (see section "Eigenvalues of updated matrix" in the book).

13. ODE: a two-step method

    Implement a two-step method for solving ODE (as in the book).

14. ODE: a more advanced step-size control

    In the homework we used the following simple condition to accept the trial step:

    ‖δy‖ < τ≡(δ+ε‖y‖)√(h/(b-a)) .

    Implement an ODE driver that accepts the step if the set of following more refined conditions is satisfied,

    |δy_k| < τ_k≡(δ+ε|y_k|)√(h/(b-a)) , k=1,...,n . where n is the number of components of the vector-function y.

    When estimating the next step-size, instead of the factor (τ/‖δy‖) use the factor (min_k |τ_k/δy_k|).

    See the chapter "ODE → Adaptive step-size control" in the book.

    A possible extension: check whether this works better than our "simple condition" on a system where one of the components is much smaller that another component. For example, the predator-prey model where the parameters are chosen such that the number of predators is much smaller than the number of prey.

15. ODE: scaling of the errors

    Investigate whether the errors on idividual steps are indeed uncorrelated and scale as square root:

    • Implement an ode-driver which caclulates the step-tolerance as

      τ=(ε‖y‖+δ)(h/(b-a))

      instead of more optimistic tolerance

      τ=(ε‖y‖+δ)√(h/(b-a))

      which we used for the homeworks.

    • Investigate which driver provides a better estimate of the integration error by integrating an ODE with known solution and comparing the error estimated by the driver with the actual error.

16. ODE: Bogacki-Shampine stepper (with FSAL property)

    Implement the (embedded) Bogacki-Shampine stepper (make sure that you utilise its FSAL property).

    See also the ODE chapter in the book.

17. Quasi-Newton's minimization with DFP update.

18. Quasi-Newton's minimization with BFGS update.

19. ODE with complex-valued functions of complex variable

    Generalize the ODE solver of your choice to solve ODE with complex-valued functions of complex variable.

20. Adaptive integration of complex-valued functions

    Implement an adaptive integrator which calculates the integral of a complex-valued function f(z) of a complex variable z along a straight line between two points in the complex plane.

21. Adaptive integrator with subdivision into three subintervals.

    Implement a (one-dimensional) adaptive integrator which at each iteration subdivides the interval not into two, but into three sub-intervals.

22. Adaptive numerical integration: error scaling

    Investigate whether the errors on sub-intervals in the adaptive numerical integration are indeed uncorrelated.

    • Implement an adaptive integrator which under recursion rescales the absolute accuracy goal with factor 2 (rather than √2).

    • Investigate which algorithm (which factor 2 or with factor √2) provides a better estimate of the error by computing some integrals with known values and comparing the errors reported by the integrator with the actual errors.

23. Clenshaw-Curtis variable transformation quadrature

    Implement an adaptive integrator which employs Clenshaw-Curtis variable substitution algorithm. Check whether it is more effective than your ordinary adaptive integrator.

24. Two-dimensional integrator

    Implement a two-dimensional integrator for integrals in the form

    ∫abdx ∫d(x)u(x)dy f(x,y)

    which consecutively applies your favourite adaptive one-dimensional integrator along each of the two dimensions. The signature might be something like

    double integ2D(double a, double b, double d(double x), double u(double x), double f(double x, double y), double acc, double eps)
    

25. Multidimensional pseudo-random adaptive recursive integrator in a rectangular volume.

    The integrator should be multi-dimensional but, for simplicity, non-stratified: simply subdivide the integration volume in all dimensions.

    Do not waste points.

    Check how far in the number of dimensions you can go before you overwhelm the computer (most probably running out of stack, but in any case up to 100M allocated memory and up to 5min "nice" run time).

26. Multidimensional pseudo-random (plain Monte Carlo) vs quasi-random (Halton and/or lattice sequence) integrators

    Investigate the convergence rates (of some interesting integrals in different dimensions) as function of the number of sample points.

27. One-sided Jacobi Singular Value Decomposition (SVD) algorithm

    Implement the one-sided Jacobi SVD algorithm as described in the "Eigenvalues" chapter of the book.

    Here is a short summary:

    • The SVD of a (real square, for simplicity) matrix A is a representation of the matrix in the form

      A = U S V^T , where matrix S is diagonal and matrices U and V are orghogonal.

    • In the cyclic Jacobi eigenvalue algorithm for symmetric matrices we applied the elementary Jacobi transformation

      A → J^T A J

      where J ≐ J(p,q) is the Jacobi matrix where the angle is chosen to eliminate the off-diagonal elements A_pq=A_qp to all upper off-diagonal matrix elements in cyclic sweeps until the matrix becomes diagonal.

    • The one-sided Jacobi SVD algorithm for general real square matrices is the same as the Jacobi eigenvalue algorithm, except that the elementary transformation is slightly different,

      A → J^T G^T A J ,

      where

      • G ≐ G(p,q) is the Givens (Jacobi) matrix where the angle is chosen such that the matrix A' ≐ G^T(p,q) A has identical off-diagonal elements, A'_pq=A'_qp: θ = atan2(A_pq - A_qp , A_qq + A_pp) (check that this is correct).

      • J ≐ J(p,q) is the Jacobi matrix where the angle is chosen to eliminate the off-diagonal elements A'_pq=A'_qp.

      The matrices U and V are updated after each transformation as

      U → U G J , V → V J .

    Of course you should not multiply Jacobi matrices, which costs O(n³) operations, but rather perform updates which only cost O(n) operations.

28. Rootfinding: 1D complex vs. 2D real

    Implement a (quasi) Newton method for rootfinding for a complex function f of complex variable z,

    f(z) = 0 .

    Hints: do it the ordinary way, just with complex numebers:

    • The usual Newton's step Δz, f(z+Δz) = 0 ,   ⇒     f(z)+f'(z)Δz = 0 ,   ⇒     Δz = -f(z)/f'(z) .
    • Instead of full Newton's step do backtracking linesearch. In the simplest incarnation, |f(z+λΔz)| < |f(z)| .

    Compare the effectiveness of your complex implementation with your homework multi-dimensional implementation of real rootfinding applied to the equivalent 2D system of two real equation with two real variables x and y,

    Re f(x + iy) = 0 ,
    Im f(x + iy) = 0 .

The functions you could use to test your 1D integrators

1. ∫01 exp(x) dx
2. ∫01 ln(x) dx
3. ∫01 x1/2 dx
4. ∫01 x3/2 dx
5. ∫01 x-1/2 dx
6. ∫01 (x>0.3?1:0) dx
7. ∫03 int(exp(x)) dx