Homework answers / question archive / Module Code: KC5000 Module Title: Further Computational Mathematics Q1

Module Code: KC5000 Module Title: Further Computational Mathematics Q1

Math

Share With

---

Module Code: KC5000

Module Title: Further Computational Mathematics

Q1. An infinite vibrating string is perturbed at t = 0. The characteristic wave speed of the string is c = 1 (in absolute units). Let u(x, t) describe the displacement of the vibrating string at position x and time t. Then the evolution in time and space of the initial perturbation can be described by means of the following wave problem

 ∂ 2u(x, t) ∂ t2 − c 2 ∂ 2u(x, t) ∂ x2 = 0 , t ≥ 0 , (1) u(x, 0) = 3 exp(−x 2 ), ut(x, 0) = −2 x exp(−x 2 ). Implement in MATLAB and use the Lax-Wendroff method to integrate numerically equation (1) with the given initial condition, u(x, 0) and ut(x, 0).

(a) Report on your exam paper a sketch of the displacement of the infinite vibrating string, u(x, t), for x ∈ [−10, 10], and at t = 0, t = 2, t = 4, t = 6, and t = 8, based on the result of the numerical integration of (1) by means of the Lax-Wendroff method, with periodic boundary conditions in [−10, 10], with Nx = 2000 space-integration intervals and Nt = 2000 time-integration intervals. Report in your exam paper all the MATLAB instructions used (you do not need to comment the code of the Lax-Wendroff method). [10 marks]

(b) Show that: u(x, t) = usol(x, t) = exp −(x − t) 2  + 2 exp −(x + t) 2  (2) is the exact solution of (1). [6 marks]

(c) Compare your numerical solution with the exact solution by plotting the difference between usol(x, t) and the numerical solution at t = 8 for x ∈ [−10, 10]. Include all MATLAB instructions in your submission. [4 marks] [Total: 20 mark
Q2. (a) Write and report in full details on your exam paper a MATLAB function file called ralston that implements the four-stage Runge-Kutta method associated to the following Butcher array 0 2 5 2 5 3 5 − 3 20 3 4 1 19 44 − 15 44 40 44 11 72 25 72 25 72 11 72 for the generic one-dimensional Cauchy problem  dy(t) dt = f (t, y(t)) y(t0) = y0 for t ∈ [t0, tend] ⊂ R , with y(t) : R → R and f(t, y) : R × R → R. The MATLAB function ralston has to take as input the generic right-hand side f(t, y) of the ordinary differential equation, the initial condition y0 = y(t0), the range of the independent variable (t0, tend), and the number of intervals N in the range, and it has to return as output the column vectors with the values of the nodes tn = t0 + n ?t, n = 0, ..., N, for the independent variable t, and the approximation un of the function y(t) at the nodes tn. [10 marks]

(b) Use the function ralston defined above to integrate numerically the following Cauchy problem ? ? ? dy dt = cos(t 2 ) y + sin(cos(y)) cos(y), for t ∈ [−5, 5] y(0) = 0 with N = 4 · 103 discretisation intervals in [0, 5]. Report on your exam paper a sketch of the function y(t) based on the numerical approximation. Use MATLAB to find and report, in format long, the position and the value of the first two local maxima featured by the numerical solution in the interval [0, 5]. Include in your submission all the MATLAB instructions used. [10 marks]
Q3. Consider the matrix function M(x) : R → R 3×3 defined as M(x) = ? ? 1 0 0 0 0 0 0 0 0 ? ? + u v t + v u t , with u = ? ? 1 x 1 ? ? and v = ? ? 1 2 −x ? ? . (3) (a) For an arbitrary diagonalisable matrix A (assumed to have all distinct eigenvalues, for the sake of simplicity), the shifted power method returns the eigenvalue of A that is the closest in modulus to a value µ (called “the shift”) specified by the user, along with its associated eigenvector. Use the shifted power method, with 2·103 iterations, tolerance 10−9 and u(0) = (1, 1, 1)t as initial guess for the eigenvector, to compute the smallest modulus eigenvalue of M(x) for x ∈ [−3, 2], taking 5000 values of x in the interval [−3, 2]. Report on your exam paper a sketch of the plot, as well as all the MATLAB instructions utilised (you do not need to include comments in the code of the shifted power method). [Hint: it is better adding the ’.’ option in the plot command to avoid joining points with segments in the plot.] [7 marks]

(b) Observe that the smallest modulus eigenvalue has a minimum for x ≈ −0.3305. Set x¯ = −0.3305. For the matrix M(¯x), plot the eigenvalue returned by the shifted power method as a function of the shift µ ∈ [−4, 4], with 2 · 103 iterations, tolerance 10−9 and u(0) = (1, 1, 1)t as initial guess for the eigenvector. Report on your exam paper a sketch of the plot obtained, as well as all the MATLAB instructions utilised (you do not need to include again the code of the shifted power method). [Hint: again, it is better adding the ’.’ option in the plot command to avoid joining points with segments in the plot.] [7 marks]

(c) Looking at the plot generated in the previous question Q3(b), choose 3 different integer values of the shift µ corresponding to the 3 different eigenvalues featured by M(¯x). Report these values in your exam paper, along with the corresponding approximation in format short of the eigenvalues and eigenvectors returned by the shifted power method, with 2 · 103 iterations, tolerance 10−9 and u(0) = (1, 1, 1)t as initial guess for the eigenvector. Include in your submission all the MATLAB instructions written to generate the output (you do not need to include again the code of the shifted power method). [Hint: after choosing the three integer values of the shift, it would be just more convenient to re-run the shifted power method for these shifts.] [If you were unable to generate the plot in the previous question Q3(b), use the shifts −3, 0, and 3.]
Q4. Given a vector function F : R n → R n , we denote by JF (x) its Jacobian evaluated at x. The problem of finding α ∈ R n such that F(α) = 0 can be solved by means of the multidimensional Newton method: given x (0) ∈ R n , for k = 0, 1, ..., until convergence solve for δx(k) : JF x (k)  δx(k) = −F x (k)  , set: x (k+1) = x (k) + δx(k) , provided that det JF x (k)  6= 0 for all k ≥ 0.

(a) Given the 2-dimensional nonlinear system  x 2 1 + x 3 2 = 1 − c 2 x 3 1 + (x2 − 1)2 = 2 − c , (4) where c is an arbitrary real parameter, find F : R 2 → R 2 such that system (4) can be written in the form F(x) = 0, and compute the Jacobian JF (x). [6 marks]

(b) In MATLAB, use the multi-dimensional Newton method for solving system (4) for c ∈ [−0.9, 0.9] with 103 iterations and tolerance ε = 10−9 . Start with c = −0.9 and set the initial condition for the Newton method to be x (0) = (0.01, −0.01)t . Take 103 values of c in the interval [−0.9, 0.9]. For each of these values of c in the interval [−0.9, 0.9], use the multi-dimensional Newton method to solve system (4) with initial condition given by the numerical solution found for the previous value of c. For each value of c in the interval [−0.9, 0.9], store the numerical solution in a 2 × 103 array. As the vector x has two components, at the end you will have a 2 × 103 array of values, representing a curve on the plane. Plot this curve and report on your submission a sketch of the plot of the curve obtained. Include in your submission all the MATLAB instructions utilised (you do not need to write line comments for the code of the Newton method).

Q5. Given the generic one-dimensional Cauchy problem y ′ (t) = f (t, y(t)) , with y(0) = y0 , the backward Euler method can be written as un+1 = un + ?t fn+1 , u0 = y0 , where un is the approximation of y(t) at the nodes tn = n ?t, ?t is the discretization step, and fn = f (tn, un). Let the functions y1(t) and y2(t) be defined by means of the following Cauchy problems: y ′ 1 (t) = cos t 3  y1(t), with y1(0) = 1 , (5a) y ′ 2 (t) = sin t 3  y2(t), with y2(0) = 2 . (5b) In the following, on MATLAB, implement and use the backward Euler method to integrate numerically (5a) and (5b).

(a) Choose N = 5000 discretization steps, and integrate numerically the Cauchy problems (5a) and (5b) in the same interval t ∈ [−5, 5]. Report on your exam paper a sketch of the numerical solution in the (y1, y2)-plane. In other words, if {u (1) n } and {u (2) n } are the numerical solutions of the Cauchy problems (5a) and (5b) respectively, then plot all the pairs (u (1) n , u (2) n ), for n = 0, 1, ..., N. Report also on your exam paper all the MATLAB instructions utilised (you do not need to comment the code of the backward Euler method). [10 marks]

(b) By plotting the numerical solutions of (5) with respect to t, we observe that y1 features several local maxima and minima for t ∈ [0, 5]. Using the numerical solution obtained above, plot on the same figure the solution y1(t) and its envelope for t ∈ [0, 5]. Include in your submission all the MATLAB instructions used. [6 marks]

(c) By overlapping y1 and y2 on the same plot, we observe that there exists a value t ∗ of t in the interval t ∈ (1/2, 1) such that y1(t ∗ ) = y2(t ∗ ). Using the numerical solution obtained above, find and report in format long on your exam paper an approximation of t ∗ , and the correspondin approximations of y1 and y2 at that point. Include in your submission all the MATLAB instructions used
Q6. Let the following Cauchy problem define the function y(t) =  y1(t), y2(t), y3(t) t , ? ??????? ??????? y ′ 1 (t) = −167 y1(t) + 134 y2(t) + 155 y3(t) + exp  cos(t) y ′ 2 (t) = −52 y1(t) + 26 y2(t) + 62 y3(t) + exp  sin(t) y ′ 3 (t) = −40 y1(t) + 93 y2(t) − 18 y3(t) , with y(0) = y 0 = ? ? 2 1 0 ? ? . (6)

(a) Write (6) in the form y ′ (t) = A y(t) + b(t), with A ∈ R 3×3 , and b(t) : R → R 3 , and show, by computing the eigenvalues of A, that the problem is stiff. [6 marks]

(b) Find the value ?tc such that, if ?t < ?tc, then the forward (explicit) Euler method, u (n+1) = u (n) + ?t f(n) = u (n) + ?t A u(n) + b , u(0) = y 0 , for solving (6) numerically is absolutely stable. [4 marks]

(c) Let Nc be the number of discretization intervals N in the interval t ∈ [0, tmax] such that tmax Nc ≥ ?tc and tmax Nc + 1 < ?tc , where ?tc is as computed above. For the Cauchy problem (6), find Nc for t in the interval t ∈ [0, 5]. Implement in MATLAB and use the forward (explicit) Euler method for solving numerically problem (6) in the interval t ∈ [0, 5] with N = Nc − 10, and N = Nc + 10 discretisation intervals. In each case, report on your exam paper the approximation found, in format long, for y(5). Include in your submission all the MATLAB instructions used (you do need to comment the code of the forward Euler method).
Q7. Let A x = b, with A = ? ? 10 −3/2 5 10 2 10 2 3/2 7 ? ? , and b = ? ? 2 5 12 ? ? , be a linear system in R 3 . In the framework of the so-called “relaxation” iterative method for solving systems of linear equations, the preconditioner P is P ≡ P(ω) = 1 ω − 1 D + L , ω 6= 0 , where ω 6= 0 is a real parameter fixed by the user, D ≡ D(A) is the diagonal matrix containing the diagonal elements of A, and L ≡ L(A) is the lower triangular matrix containing the elements in the lower triangular part of

A. Observe that, when ω = 1, the preconditioner coincides with that of the Gauss-Seidel method, P = L. For any choice of ω 6= 0, let B(ω) be the iteration matrix corresponding to the preconditioner P(ω), B(ω) = P(ω) −1  P(ω) − A . (a) Let Λ(B(ω)) indicate the spectrum of B(ω). Sketch the plot of the spectral radius ρ(B(ω)) = max λ∈Λ(B(ω)) |λ| as a function of ω ∈ [0.2, 1.12], taking 500 values of ω in that interval. Include all the MATLAB instructions written to generate the output in your submission. [Hint: Construct the two MATLAB functions PR = @(w) ((1/w)-1)*diag(diag(A))+tril(A) and BR = @(w) inv(PR(w))*(PR(w)-A); divide the interval in 500 parts; the spectral radius can be computed by means of the instruction max(abs(eig(BR(w)))), which has to be used inside a for loop for generating the data for the plot]. [7 marks] (b) Use the MATLAB built-in command fminsearch (with the default set of options for accuracy, termination tolerances, and maximum number of iterations) for evaluating the minimum of ρ(B(ω)) as a function of ω ∈ [0.2, 1.2], using the point ω = 0.95 as initial condition. Let ω¯ be the minimizer. Report in format long the value of the minimizer ω¯ and the corresponding minimum ρ(B(¯ω)). Include in your submission all the MATLAB instructions written to generate the output. [3 marks] (c) Implement in MATLAB and apply the “relaxation” iterative method for ω¯, in the form x (k+1) = B(¯ω) x (k) + g , with g = P(¯ω) −1 b , for solving the linear system A x = b, with 1000 iterations, initial condition x (0) = (1, 1, 1)t , and stopping criterion given by the relative residual   
A x(n) − b     2 / ||b||2 with tolerance ε = 10−12. Report on the exam paper your MATLAB code in full details, and the numerical solution in format long. Report also the number of iterations required by the algorithm to converge. [If you were unable to obtain an estimate for ω¯ in Exercise 7(b), solve the present Exercise 7(c) setting ω¯ = 0.95]