Computational Practicals

Numerical and Computational Methods II (M´etodos Num´ericos y Computacionales II) Engineering Physics (UPC - ETSETB) Computational Practicals ` A. Meseguer Enginyeria F´ısica Dept. F´ısica Aplicada (UPC) Numerical and Computational Methods II (2015-2016) A. Meseguer & F. Marques Numerical lab 1: Numerical Linear Algebra I (LU factorization). 1. The figure below shows an increasing series of circuits consisting of identical resistors (R= 3 Ω) connected to a battery V= 1 V. The circuit at the bottom is a generalization of the first 3 circuits. We want to compute the intensities Ikof each resistor as well as the equivalent resistance Re between points aand bof the n-ladder system for arbitrary n≥2. 11 2 123 123n a b a b a b a b I1I3I5I2n−3 I2n−1 I2I4I6 I1I3I5I2n−3 I2n−2 (a) Kirchoff: convince yourself that the top and bottom intensities at each step ladder are the same (so we only have 2n−1 unkowns: I1, . . . , I2n−1), and that: 2R I1+R I2=V I1−I2−I3= 0 −R I2+ 2R I3+R I4= 0 I3−I4−I5= 0 . . .. . .. . .. . . −R I2n−4+ 2R I2n−3+R I2n−2= 0 I2n−3−I2n−2−I2n−1= 0 −R I2n−2+ 3R I2n−1= 0. If we define I,V∈R2n−1as the vector of unknown intensities I≡[I1, I2, . . . , I2n−1]Tand the right hand side vector associated with the battery V≡[V, 0,...,0]T, respectively, we can rewrite the above system in matrix form:                 2R R 0 0 0 0 ··· 1−1−1 0 0 0 ··· 0−R2R R 0 0 ··· 0 0 1 −1−1 0 ··· . . .. . . 0 ......... . . .. . .. . .. . .−R2R R 0 . . .. . .. . .. . .. . . 1 −1−1 . . .. . .. . .. . .. . . 0 −R3R                                I1 I2 I3 I4 I5 . . . I2n−3 I2n−2 I2n−1                =                V 0 0 0 0 . . . 0 0 0                . In Matlab, the matrix above can be generated with the commands repmat and spdiags: B = [repmat([-R;1],2*n-1,1), repmat([2*R;-1],2*n-1,1),repmat([R;-1],2*n-1,1)]; A = full(spdiags(B,[-1:1],2*n-2,2*n-1)) ; A = [A ; [zeros(1,2*n-3) -R 3*R]]; (b) Using your own pa=lu code, solve the system above for n= 10,20 and 30. For each case, plot (k, Ik) and observe that the intensities decay exponentially (try semilogy). (c) Compute Re(n) for the cases above and check that l´ım n→∞ Re(n) = R(1 + √3). Enginyeria F´ısica Dept. F´ısica Aplicada (UPC) Numerical and Computational Methods II (2015-2016) A. Meseguer & F. Marques Numerical lab 2: Numerical Linear Algebra II (qr factorization and the Least Squares Problem). The figure on the right shows the function: f(x) = e−20(x+0.25)2+ 0.25 sin(30 x) e−20(x−0.25)2, within the domain [−1,1]. In mnc1, you studied the polynomial interpolation technique to approximate functions by imposing an interpolant Πnf(x)∈ Rn[x] to adopt the values of the function at n+ 1 given nodes x0, x1,...,xn, i.e.: Πnf(xj) = f(xj), j = 0,1,...,n. x y -1 -0.5 0 0.5 1 -0.2 0 0.5 1 In this case you have to approximate f(x) with polynomial expansions of arbitrary degree by taking an also arbitrary (and greater) number of sampling points. (a) Consider the approximation p(x)∈Rm−1[x] based on a linear combination of the first m−1 monomials 1, x, x2,...,xm−1: p(x) = m−1 X j=0 ajxj=a0+a1x1+a2x2+···+am−1xm−1, and the values fj=f(xj) of the function f(x) evaluated at the set of nequispaced nodes in [−1,1]: xj=−1+2j/(n−1), j = 0,1,...,n−1, with n > m. If we define a= [a0, a1,...,am−1]T and f= [f0, f1,...,fn−1]T, the problem has mcoefficients to be determined and nconditions to be satisfied by the polynomial. Since n > m, imposing p(xk) = f(xk), k = 0,1,...,n−1, leads to an overdetermined Vandermondetype n×msystem Va=fwhose explicit form is:        1x0x2 0··· xm−1 0 1x1x2 1··· xm−1 1 1x2x2 2··· xm−1 2 . . .. . .. . .. . . 1xn−1x2 n−1··· xm−1 n−1               a0 a1 a2 . . . am−1        =        f0 f1 f2 . . . fn−1        The Vandermonde matrix above can be easily constructed by using vander(x) Matlab’s builtin command, where xis the vector of nodes. Then apply the fliplr (flip left-right) on the resulting matrix to reorder the columns and drop the m+ 1,...,n columns with V=V(:,1:m). The vector athat minimizes the norm of the residual kf−Vakcan be computed by means of the qr factorization. Use your own qr-code to solve the Vandermonde system for the cases (n, m) = (14,7),(28,14),(28,20),(64,30). For each case, plot the least squares polynomial on a finer grid. Observe the progressive improvement in capturing the small oscillations on the right. You will also observe a familiar phenomenon near the ends of the interval. (b) Optional: repeat the computation with Chebychev points xj= cos jπ n−1, j = 0,1,...,n−1. Enginyeria F´ısica Dept. F´ısica Aplicada (UPC) Numerical and Computational Methods II (2015-2016) A. Meseguer & F. Marques Numerical lab 3: Numerical Linear Algebra III (Gram-Schmidt reorthogonalization). In this experiment, you will compare three different methods of orthogonalization: classical GramSchmidt (cgs), Gram-Schmidt with reorthogonalization (gsr) and qr. You will orthogonalize Hilbert matrices H(n) that have the explicit elements Hij =1 i+j−1, i, j = 1,2, . . . , n. A Hilbert matrix is generated with the Matlab command hilb(n), where nis the dimension of the matrix. You need three external codes to orthogonalize the matrix Hso they provide the Q orthogonalized columns of the matrix. To measure the ’quality’ of the orthogonalization, you will compute the 2-norm kI−QTQk2, where Iis the identity matrix. If the orthogonalization is well done, this norm should be very small. To do that, you will make use of the norm command (type help norm to see which of the norms you can ask Matlab to compute). (a) cgs: generate H(n) (with n= 9 to n= 15, for example) and monitor the inner products hqj,˜qki for j≤k−1 of the resulting normalized basis. At every stage k, check the norm of ˜qkand observe how it is decreasing. Do you observe a deterioration of kI−QTQk2when you increase n?. (b) gsr: same as in (a) but now monitor the non-orthogonal roundoffs sj=hqj,˜qki, j ≤k. Substract these projections as seen in the lectures and recompute the q(2) kvectors. Re-check orthogonality. Do you observe an improvement when compared with cgs?. (c) qr: repeat with Matlab’s qr. Which of the three methods performs better? (d) Optional: compute the condition number of the Hilbert matrices using the cond(H) command (again, type help cond to check in which norm is this condition number computed). A geometrical interpretation of the illconditioning of a matrix is based on the lack of orthogonality of its columns (that generate the rank of the represented linear map). You can visualize this phenomenon by: (1st) normalizing each column of H, (2nd) computing HTH, i.e., (HTH)ij = cos(θij), where θij is the ’angle’ between the iand jcolumns of H: for j = 1:n H(:,j) = H(:,j)/norm(H(:,j)); % Normalization end CS = flipud(H’*H) ;% CS contains the cosine of the angles contour(CS,[0.9:.01:1.0]); colorbar ; % Contours of cos(angle) in [0.9,1.0] Whereas a well-conditioned matrix should have the contours concentrated nearby the diagonal, an ill-conditioned one will show a spread of this contour lines along all its elements. Repeat the plot with the Qmatrix computed with cgs,gsr and qr. Enginyeria F´ısica Dept. F´ısica Aplicada (UPC) Numerical and Computational Methods II (2012-2013) A. Meseguer & F. Marques Numerical lab 4: Numerical Linear Algebra IV (Krylov iterative methods–GMRES). The figure below shows the same circuit as in the first practical: identical resistors (R= 3 Ω) connected to a battery V= 1 V. We want to compute the intensities Ikup to a fixed tolerance and without explicitly constructing the matrix. As seen in the first practical, the intensities were decreasing from left to right in an exponential fashion. Therefore, we may not be interested in the whole O(m)-vector of intensities [I1I2··· I2m−1]T, but just in its first leading O(n)-components (the “tail ” of this vector will contain Ij’s whose magnitude is below our prescribed tolerance, thus negligible). That is, we are going to look for O(n)-low dimensional approximations of an O(m) dimensional problem. 12 3 m a b I1I3I5I2m−3 I2m−1 I2I4I6 I1I3I5I2m−3 I2m−2 (a) Kirchoff: the solution vector of intensities I= [I1I2··· I2m−1]Thas to satisfy the 2m−1 equations 2R I1+R I2=V, I1−I2−I3= 0, 2R I2j−3+R I2j−2−R I2j−4= 0,(j= 3,...,m), I2j−3−I2j−2−I2j−1= 0,(j= 3,...,m), 3R I2m−1−RI2m−2= 0. In the first practical, the left hand side of the equations above were seen as a matrix-vector product. For the Krylov method, we only need the action of the matrix on the vector: I−→ AI. Generate a Matlab function that for a given vector I= [I1I2··· I2m−1]T, provides the left hand side of the equations above. Example for m= 15: function Ax = Afun(x) R = 3 ; V = 1 ; m = 15; Ax = [] ; Ax = [Ax ; 2*R*x(1) + R*x(2)]; Ax = [Ax ; x(1) - x(2) - x(3)] ; for j = 3:m Ax = [Ax ; 2*R*x(2*j-3) + R*x(2*j-2) - R*x(2*j-4)]; Ax = [Ax ; x(2*j-3) - x(2*j-2) - x(2*j-1)] ; end Ax = [Ax ; 3*R*x(2*m-1) - R*x(2*m-2)] ; The output of this function is the vector Ax. Check that your output is the same as the one provided by the code of the first practical. (b) Using your own gmres code, solve the system above for tolerances kAx−bk ≤ 10−3, 10−4and 10−5. Check that the first components of the solution vector coincide with the ones provided by the matrix version. Also check that the number of Arnoldi iterations grows as long as you decrease the tolerance. Note: your gmres code has to have an input-output syntax like this: [x,k] = mygmres(Afun,b,tol,dimkryl), where Afun is the action (invoked with @Afun) function above, bis the right hand side of the system, tol is the tolerance and dimkryl is the dimension of the Krylov space where your approximation is been computed, i.e., n, so in this case dimkryl ≤2m−1. Optional: type help gmres in Matlab. Enginyeria F´ısica Dept. F´ısica Aplicada (UPC) Numerical and Computational Methods II (2015-2016) A. Meseguer & F. Marques Numerical lab 5: Systems of Nonlinear Equations (Newton method and continuation). The figure on the right shows a double pendulum consisting of two small spheres of mass mconnected by massless rigid rods of equal length `. The system is forced to rotate with a constant angular speed ωaround a vertical axis passing trough the pivot p. The goal of this practical is to get familiar with the concept of multiplicity of solutions arising from nonlinear physical systems. In particular, in this exercise we want to find the possible equilibrium solutions for a given angular speed, i.e., configurations where the angles φ1and φ2remain constant in time. Do it before starting the practical !: apply Newton’s laws and impose equilibrium in a reference frame co-rotating with the two masses. By dividing radial and vertical equations for each mass you should get: F1(φ1, φ2, α) = tan(φ1)−α(2 sin φ1+ sin φ2)=0 F2(φ1, φ2, α) = tan(φ2)−2α(sin φ1+ sin φ2)=0, where α=`ω2/2g is a dimensionless parameter. For simplicity, we will assume that φ1∈[0, π/2) and φ2∈(−π/2, π/2). m m ℓ ℓ φ1 φ2 ~g ω p (a) The configuration (φ1, φ2) = (0,0) is always an equilibrium solution ∀α. For α∈[0,2], compute the determinant of the jacobian det JF(φ1, φ2) of the system F= [F1F2]Tat (φ1, φ2) = (0,0). Plot det JF(0,0) as a function of α. Do you find any indication that new branches of equilibrium solutions may emerge for some value(s) of α?. Compute the jacobian matrix via finite differences (not analytically) and use the Matlab command det to compute the determinant. (b) Newton exploration: for fixed values of α∈[0,2], start your Newton method from initial random seeds (φ(0) 1, φ(0) 2)∈[0, π/2) ×(−π/2, π/2) and check if the iteration converges to other solutions. If that is the case, store only those converged solutions which fall within the interval (φ1, φ2)∈ [0, π/2) ×(−π/2, π/2). Plot the converged solutions φ1(α) and φ2(α) in two different figures. (c) Continuation (optional): starting from any of the new solutions found, continuate them as a function of αand plot them in two different figures. At what values of αare these solutions born? Enginyeria F´ısica Dept. F´ısica Aplicada (UPC) Numerical and Computational Methods II (2013-2014) A. Meseguer & F. Marques Numerical lab 6: Gauss-Legendre approximation. In this practical you are asked to code an approximation scheme for the Gauss-Legendre nodes. This scheme has to provide: (1) gl nodes {x0,···, xn}, (2) quadrature weights {w0,...,wn}, (3) barycentric weights {λ0,...,λn}and (4) differentiation matrix Dij for an arbitrary number N=n+ 1 of nodes. Code a Matlab function [x,w,lamb,D] = myleg(N) that provides the above elements for a given N. (a) Nodes and weights: follow the steps of Problem 3, using the companion matrix technique. (b) Barycentric weights: once you have computed {xj}and {wj}for j= 0,...,n, use the the magical formula to obtain the barycentric weights: λk= (−1)kq(1 −x2 k)wk. That is, if you sample an arbitrary function f(x) at the gl nodes, the interpolating polynomial (in barycentric form) is: Πnf(x) = n X j=0 λjfj(x−xj)−1 n X j=0 λj(x−xj)−1 ,for x6== xjand with fj=f(xj), j = 0,...,n. (c) Differentiation matrix: recall from mnc1 that, in general: Dij =       λj λi 1 xi−xj (i6=j), X k6=j 1 xj−xk (i=j). Test your matrix by numerically differentiating Runge’s function f(x) = (1 + x2)−1and comparing with its exact derivative f′(x) on the interval [−1,1]. If we define Dfjas the numerical derivative of f(x) at node xj, plot on a semilogy figure xjversus |Dfj−f′(xj)|. Increase nand check that there is no Runge instability and that the convergence is geometrical. For higher values of n, test your differentiation matrix on the function you used in Practical 2: f(x) = e−20(x+0.25)2+ 0.25 sin(30 x) e−20(x−0.25)2, whose analytical derivative (in Matlab format) is: df = -((40*(x + .25))./exp(20*(x+.25).^2))-... (10*(x - 1/4).*sin(30*x))./exp(20*(x-1/4).^2)+... ((15/2)*cos(30*x))./exp(20*(x-1/4).^2); (d) Optional: develop a similar routine for the Chebychev nodes. This is easier than before because you have already studied the barycentric weights and differentiation matrix for these nodes in mnc1. But remember that the quadrature weights are only useful for integrals of the form Z1 −1 f(x) √1−x2dx. You will need both routines for the incoming practicals. Enginyeria F´ısica Dept. F´ısica Aplicada (UPC) Numerical and Computational Methods II (2016-2017) A. Meseguer & F. Marques Numerical lab 7: Matrix Exponentiation, Eigenvalues and Discrete Fourier Transform. The figure on the right shows a system of masses connected via linear springs. Let xi(t) and ˙xi≡vi(t) be the displacement of the i-th block from its equilibrium position and its k1k2k3k4 m1m2m3 x1x2x3 instantaneous speed at time t, respectively. You can verify that:         ˙x1 ˙x2 ˙x3 ˙v1 ˙v2 ˙v3         =         0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 β11 β12 0 0 0 0 β21 β22 β23 000 0β32 β33 000                 x1 x2 x3 v1 v2 v3         ,with            β11 =−(k1+k2)/m1 β12 =k2/m1, β21 =k2/m2 β22 =−(k2+k3)/m2 β23 =k3/m2, β32 =k3/m3 β33 =−(k3+k4)/m3. Let z(t) = [x1x2x3v1v2v3]T∈R6be the state vector at time t, and Bthe 6×6 matrix of coefficients of the ode system above. For a given arbitrary initial condition z0=z(0), the solution is z(t) = eBtz0. Random initial positions xi(0) and velocities vi(0) will in general lead to coupled oscillations with some characteristic frequencies. The purpose of this practical is to identify these natural frequencies with two different methods. Numerical values: take, for example: mi= 1, k1= 1, k2=√2, k3=√3, k4= 4,z0= [0.281 0.033 −1.33 1.12 0.35 −0.299]T. (a) Method I: according to what you studied in mechanics, the normal modes or frequencies are the imaginary parts of the spectrum of eigenvalues of B. Using Matlab’s command eig, compute the eigenvalues of B. For the values above, you should get: ω1= 1.02236572722, ω2= 1.8826701636, ω3= 2.5889863102. (b) Method II: for a sampling set of times tj= j∆t, j = 0,1,2,...,N −1, you can compute zj=z(tj) = eBtjz0’exactly’ by using Matlab’s command expm(B*tj) at every tj. However, the command expm is computationally very expensive. You just need to compute eB∆tonce: zj= ej∆tBz0= e∆tBe(j−1)∆tBz0= e∆tBzj−1. Integrate the system above with ∆t= 0.25 and N= 400 and store x2(ti) (or any other of the variables). Plot x2(t) to visualize the oscillation. Using your own dft routine, compute the frequency spectrum of x2and represent it on a semilogy plot. Your dft should show clear peaks at the normal frequencies (eigenvalues of method i) as seen in the plots. Therefore, by analyzing the oscillation of a physical system we can infer its internal properties (here the spring constants kiand masses mj). In practice, this type of time frequency (or spatial wavelength) analysis is used in many branches of Physics in order to identify the inner structure of matter. Of course, this practical is a very simplified model. Enginyeria F´ısica Dept. F´ısica Aplicada (UPC) Numerical and Computational Methods II (2015-2016) A. Meseguer & F. Marques Numerical lab 8: Boundary Value Problem (Quantum Mechanical Square Well) The figure on the right shows the square potential well: V(x) = V0,|x|> a/2 0,|x|< a/2. By introducing the change of variable u=x/a, the Schr¨odinger equation: d2ψ dx2=2m ~2[V(x)−E]ψ, V(x) = V0 V= 0 a/2 −a/20x E can be written as: −d2φ du2+βΘ(|u| − 1/2)φ=βεφ, with: φ(u)≡ψ(x(u)), β =2ma2V0 ~2,ε≡E V0 and Θ(z) = 1, z > 0 0, z < 0. Numerical values: use the same set of parameters as in Eisberg & Resnick (see Appendix G), i.e., β= 64. Set L≈5 in your transformation map from (−1,1) −→ (−∞,+∞). (a) Compute the energies εnof the three bounded eigenstates. Represent the associated φneigenfunctions. (b) Compare your results with the ones given by the analytical solution obtained in Eisberg & Resnick (Appendix H). To do that, remember that the normalized energies Endefined by e & rare: En=rmEna2 2~2=1 2pεnβ, and they have to satisfy either one of the two trascendental equations (that you can solve by means of a simple Newton method): Entan En=rβ 4− E2 nor Encot En=−rβ 4− E2 n. You will need a very high number of points (N∼256 or more). This is due to the fact that the potential is unrealistic and mathematically ill-defined (discontinuous), therefore we cannot expect geometrical convergence.