October 7, 2002 1Fall 2002. 10.34. Numerical Methods Applied to Chemical EngineeringHomework # 4. Nonlinear algebraic equations and matrix eigenvalue problemsAssigned Monday 10/6/02. Due Friday 10/11/02Problem 1. Calculation of vibrational modes with eigenvalue analysisIn this homework, you will use a program that computes the normal (vibrational) modes of a 2D lattice (see figure below) using eigenvalue analysis. A MATLAB program will be posted on the web site that calculates each vibrational mode and natural frequency for this system. Another program will be provided that takes this information and makes a “movie” to animate a selected vibrational mode. You will be asked to use these programs to compare the natures of the high and low energy vibrational modes of this structure. No programming will be required!FIGURE 1. 2-D lattice of point mass-spring systemIn this problem, we wish to investigate the vibrational modes of an 2-D lattice of point masses connected by springs. We assume that the lattice is infinite, but we only want to simulate a small section of N “rows” and N “columns” of points, shown in the figure above as the black, filled circles. To mimic the effect of having an infinitely large lattice of points, we employ periodic boundary conditions in which we assume that the system is bounded by perfect copies of itself on all sides (see figure below). When one of the “atoms” in our system moves to the left, its image also moves to the left. In the figure below, curved arrows link the atoms in the system (filled circles) to their images outside ofOctober 7, 2002 2the system (open circles). For further details on how the periodic boundary condition is employed, consult the attached MATLAB program. FIGURE 2. Linking of atoms and their images in period boundary conditions to mimic an infinite latticeWe now number the lattice sites from 1 to as shown in the figure below. FIGURE 3. Numbering system for lattice sites, including identification of neighbor site bondingSN2=October 7, 2002 3While we will use the master indices for labelling purposes, it is easier to write the poten-tial energy function for the system if we label the site in “column” i and “row” j as . The position of this site is(EQ 1)We write the total potential energy of the system as a sum of harmonic springs bonding neighboring lattice sites,(EQ 2)A factor of is added before the summation to correct for the overcounting of each spring twice. The distance between site and the image of site with which it interacts is(EQ 3)where we implement the periodic boundary conditions (figure 2) as(EQ 4)(EQ 5)We subtract from each bond distance the preferred (zero-energy) value of so that at mechanical equilibrium, all lattice sites are separated by this distance. The position in the minimum energy state of the lattice site - that in “column” and “row” with master index - isij,()rij,()xij,()yij,()=Ur11,()… rNN,(),,()12---12---Kri 1– j,()ij,()lb–()2ri 1+ j,()ij,()lb–()2+[j 1=N∑i 1=N∑=rij 1–,()ij,()lb–()2rij 1+,()ij,()lb–()2+ ]+12⁄ij,()mn,()rmn,()ij,()∆xmn,()ij,()()2∆ymn,()ij,()()2+[]12⁄=∆xi 1– j,()ij,()xij,()xi 1– j,()– i 2 N,[]∈,x1 j,()xNj,()– Nlb+ i, 1==∆xi 1+ j,()ij,()xij,()xi 1+ j,()– i 1 N 1–,[]∈,xNj,()x1 j,()– Nlb– i, N==∆yij 1–,()ij,()yij,()yij 1–,()– j 2 N,[]∈,yi 1,()yiN,()– Nlb+ j, 1==∆yij 1+,()ij,()yij,()yij 1+,()– j 1 N 1–,[]∈,yiN,()yi 1,()– Nlb– j, N==lbij,()ijki1–()Nj+=October 7, 2002 4(EQ 6)We wish to compute the vibrational modes of the lattice and their natural frequencies. To do so, we pack all of the degrees of freedom of the system into a single vector of length , . We also define a vector of the forces on each degree of freedom that is the negative of the gradient of the potential energy. (EQ 7)For the potential energy function above, the derivatives of the potential energy with respect to the lattice site positions are computed in a straight-forward, if tedious, manner.(EQ 8)rˆij,()xˆij,()yˆij,()i 1–()lbj 1–()lb==F 2S=qfqq1q2q3q4:qF 1–qFx11,()y11,()x21,()y22,():xNN,()yNN,()== f ∇U–q1∂∂Uq2∂∂Uq3∂∂Uq4∂∂U:qF 1–∂∂UqF∂∂U–x11,()∂∂Uy11,()∂∂Ux21,()∂∂Uy21,()∂∂U:xNN,()∂∂UyNN,()∂∂U–== =xij,()∂∂UK2----xij,()∂∂ri 1+ j,()ij,()lb–()2ri 1– j,()ij,()lb–()2rij 1+,()ij,()lb–()2rij 1–,()ij,()lb–()2+++[]=xij,()∂∂UK2----2 ri 1+ j,()ij,()lb–()xij,()∂∂ri 1+ j,()ij,()2 ri 1– j,()ij,()lb–()xij,()∂∂ri 1– j,()ij,()+=2 rij 1+,()ij,()lb–()xij,()∂∂rij 1+,()ij,()2 rij 1–,()ij,()lb–()xij,()∂∂rij 1–,()ij,()++xij,()∂∂UKri 1+ j,()ij,()lb–()∆xi 1+ j,()ij,()ri 1+ j,()ij,()----------------------ri 1– j,()ij,()lb–()∆xi 1– j,()ij,()ri 1– j,()ij,()----------------------+=rij 1+,()ij,()lb–()∆xij 1+,()ij,()rij 1+,()ij,()----------------------rij 1–,()ij,()lb–()∆xij 1–,()ij,()rij 1–,()ij,()----------------------++October 7, 2002 5(EQ 9)We define also a matrix of size of the second derivatives of the energy with respect to each degree of freedom, the so-called Hessian matrix , with the elements(EQ 10)Here we see that because the order of differentiation does not matter, the Hessian matrix is necessarily real and symmetric, so that all eigenvalues of are real and all eigenvectors of are mutually orthogonal. (EQ 11)Rather than derive the Hessian analytically, we will compute the value of the Hessian - evaluated at using finite difference approximations,(EQ 12)For each of the degrees of freedom in the system, we change the value of by some small value and evaluate again the vector of first derivatives. Application of the formula above provides a numerical estimate of the values in the column of the Hessian matrix.In general, we are given the potential energy function, and must find the state of minimum energy numerically. Below, we show in an optional discussion section how one may mod-ify Newton’s method to find the minimum. For this system, we can identify a priori the values of the generalized coordinates at the minimum energy state, . yij,()∂∂UKri 1+ j,()ij,()lb–()∆yi 1+ j,()ij,()ri 1+ j,()ij,()----------------------ri 1– j,()ij,()lb–()∆yi 1– j,()ij,()ri 1– j,()ij,()----------------------+=rij 1+,()ij,()lb–()∆yij 1+,()ij,()rij 1+,()ij,()----------------------rij 1–,()ij,()lb–()∆yij