1.MATH347 Homework 4

1.1.Background: modeling binary interaction in science and engineering

Many systems consist of point interactions. The iconic Eiffel Tower can be modeled through the trusses linking two points with coordinates ${𝒙}_i,{𝒙}_j\in {ℝ}^3$. The unit vector along direction from point $i$ to point $j$ is

Suppose that structure deformation changes the coordinates to ${𝒙}_i+{𝒖}_i$, ${𝒙}_j+{𝒖}_j$. The most important structural response is that on point $i$ a force ${𝒇}_{ij}$ is exerted that can be approximated as

in appropriate physical units. An opposite force ${𝒇}_{ji}=-{𝒇}_{ij}$ acts on point $j$. Note the projector along the truss direction ${𝒍}_{ij}{𝒍}_{ij}^T$. Adding forces on a point from all trusses leads to the linear relation

with $𝒇,𝒖\in {ℝ}^m$, $𝑲\in {ℝ}^{m\times m}$, $d=3$ the number of dimensions, $N$ the number of points, and $m=dN$, the number of degrees of freedom of the structure. Point interactions arise in molecular or social interactions, cellular motility, and economic exchange much in the same form (1). The techniques in this homework are equally applicable to physical chemistry, sociology, biology or business management.

Load the Eiffel Tower data.

octave] 

cd /home/student/courses/MATH347DS/data/EiffelTower;
load EiffelPoints; X=Expression1; [NX,ndims]=size(X);
load EiffelLines; L=Expression1; [NL,nseg]=size(L);

octave] 

m=3*NX; disp([NX NX^2 NL m])

octave] 

Of the possible $N^2\cong 7\times {10}^8$ linkages between points in the structure, only $N_l=31462$ are present. It is wasteful, and often impossible, to store the entire matrix of couplings between points $𝑲\in {ℝ}^{m\times m}$ ($m^2\cong 6\times {10}^9$), but storing only the nonzero elements corresponding to the interacting can be carried out through a sparse matrix. After allocating space with spalloc, the matrix is initialized with zeros along the truss directions.

octave] 

K=spalloc(m,m,3*NL);
for k=1:NL
  i = L(k,1); j = L(k,2);
  p = 3*(i-1)+1; q = 3*(j-1)+1;
  for id=1:3
    for jd=1:3
      K(i+id,j+jd) = 0;
    end;
  end;
end

octave] 

The matrix $𝑲$ can be formed by loops over the trusses. The following computation takes a few minutes. The matrix $K$ has been pre-computed, and can be loaded from a disk file when answering the questions below, rather than executing the following loops.

octave] 

for k=1:NL
  i = L(k,1); j = L(k,2);
  lij = (X(j,:)-X(i,:))'; Pij = lij*lij';
  p = 3*(i-1)+1; q = 3*(j-1)+1;
  for id=1:3
    for jd=1:3
      K(i+id,j+jd) = K(i+id,j+jd) + Pij(id,jd);
    end;
  end;
end

octave] 

save "K.mat" K

octave] 

1.2.Numerical experimentation of theoretical questions

Use small matrices in Octave to find evidence or counter-examples for these true or false questions.

1. If $𝑹$ is an upper-triangular matrix, its singular values are the same as the diagonal values $r_{ii}$.

2. If $𝑨={𝑨}^T$, the singular values of $𝑨$ are the same as its eigenvalues.

3. If $𝑻$ is of full rank, $𝑨$ and $𝑩={𝑻}^{-1}𝑨𝑻$ have the same singular values.

4. For $𝑨,𝑩\in {ℝ}^{m\times m}$, the matrices $𝑨𝑩$ and $𝑩𝑨$ have the same eigenvalues.

1.3.Reduced order modeling

Load data, and define functions to draw the deformed structure.

octave] 

size(K)

octave] 

octave] 

cd /home/student/courses/MATH347DS/data/EiffelTower;
load EiffelPoints; X=Expression1; [NX,ndims]=size(X);
load EiffelLines; L=Expression1; [NL,nseg]=size(L);
m=3*NX; load K;

octave] 

function drawXu(X,u,L,is)
  clf; view(-30,45); [NX,nd]=size(X);
  [nL,nseg]=size(L); x=[]; y=[]; z=[];
  U=reshape(u,3,NX)';
  for k=1:is:nL
    x=[x X(L(k,1),1)+U(L(k,1),1) X(L(k,2),1)+U(L(k,2),1)];
    y=[y X(L(k,1),2)+U(L(k,1),2) X(L(k,2),2)+U(L(k,2),2)];
    z=[z X(L(k,1),3)+U(L(k,1),3) X(L(k,2),3)+U(L(k,2),3)];
  end;
  plot3(x,y,z);
end

octave] 

drawXu(X,zeros(m,1),L,20); axis 'equal'

octave] 

function drawXU(X,U,L,is)
  clf; view(-30,45); [NX,nd]=size(X);
  [nL,nseg]=size(L); x=[]; y=[]; z=[];
  for k=1:is:nL
    x=[x X(L(k,1),1)+U(L(k,1),1) X(L(k,2),1)+U(L(k,2),1)];
    y=[y X(L(k,1),2)+U(L(k,1),2) X(L(k,2),2)+U(L(k,2),2)];
    z=[z X(L(k,1),3)+U(L(k,1),3) X(L(k,2),3)+U(L(k,2),3)];
  end;
  plot3(x,y,z);
end

octave] 

drawXU(X,zeros(NX,3),L,20); axis 'equal'

octave] 

octave] 

The following questions only require a few Octave lines to solve correctly. Be careful to not request output of the large vectors and matrices that arise. Use visualization with the draw function to represent results, and concentrate on the mathematical formulation. For example, imposing displacement $𝒖=(z,0,0)$ at each node, i.e., a horizontal displacement along the $x$ direction proportional to the $z$ coordinate is accomplished by:

octave] 

U=zeros(NX,3); U(1:NX,1)=X(1:NX,3);

octave] 

drawXU(X,U,L,20)

octave] 

1.3.1.Least squares solution on singular modes

Find the linear combination $𝒔$ of the first 10 singular modes that best approximates a force $𝒇=(z\cos \theta ,z\sin \theta ,0)$ on the structure. Randomly choose $\theta$. Show the first 3 singular modes and the linear combination you find.

1. The displacements $𝒖$ resulting from the application of a force $𝒇$ are found by solving the linear system

   ${\displaystyle 𝒇=𝑲𝒖,}$

   which however is costly since $𝑲$ is of large size.

2. Find an approximation based upon the dominant part of $𝑲=𝑼𝚺{𝑽}^T$, namely choosing $n=10$ dominant singular modes

   octave]

   n=10; [U S V] = svds(K,n);

3. Find displacements for the $n=10$ dominant modes by least squares

   ${\displaystyle {min}_{𝒔\in {ℝ}^n}||𝒇-𝑼𝒔||}$

4. Solve by projection

octave] 

1.3.2.Least squares solution on eigenmodes

Find the linear combination $𝒕$ of the first 10 eigenmodes that best approximates a force $𝒇=(z\cos \theta ,z\sin \theta ,0)$ on the structure. Show the first 3 eigenmodes and the linear combination you find.

1.3.3.Least squares solution on Krylov modes

Replace singular modes or eigenmodes by orthogonal $𝑸$, with $𝑸𝑹=\left[\begin{array}{ccccc}𝒇& 𝑲𝒇& {𝑲}^2𝒇& \dots & {𝑲}^9𝒇\end{array}\right]$ (known as Krylov modes). Find the linear combination $𝒘$ of the first 10 Krylov modes that best approximates a force $𝒇=(z\cos \theta ,z\sin \theta ,0)$ on the structure. Show the first 3 Krylov modes and the linear combination you find.

1.3.4.Reduced-rank singular mode approximation

From the SVD $𝑲=𝒀𝚺{𝒁}^T$ construct an approximate matrix

and compare the force $𝒇$ with $\stackrel{\sim }{𝒇}=\stackrel{\sim }{𝑲}𝒔$, $𝒔$ from question 1.3.1.

1.3.5.Reduced-rank eigenmode approximation

From the eigendecomposition $𝑲𝑹=𝑹𝚲$ construct an approximate matrix

and compare the force $𝒇$ with $\stackrel{\sim }{𝒇}=\stackrel{\sim }{𝑲}𝒕$, $𝒕$ from question 1.3.2.

1.3.6.Reduced model convergence

Choose one of the reduced models from 1.3.1 to 1.3.4. Repeat the calculations for 15,20,25 modes, and comment on whether the solutions converge.