MATH661 Project 1Large system factorization and model reduction

This project presents a typical medium-size problem involving the concepts from approximation using linear combinations: transformation of coordinates, least squares, different types of basis sets.

1Introduction

1.1Modeling 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

Suppose that structure deformation changes the coordinates to

𝒙_{i} + 𝒖_{i}

𝒙_{j} + 𝒖_{j}

. The most important structural response is a force

𝒇_{i}

exerted on point

i

that can be approximated as

in appropriate physical units. An opposite force

𝒇_{j} = - 𝒇_{i}

acts on point

j

by principle of action and reaction. Note the projector along the truss direction

𝑷_{i j} = 𝒍_{i j} 𝒍_{i j}^{T}

, allowing computation of nodal forces as

with

𝒇, 𝒖 \in ℝ^{m}

the forces and displacements,

𝑲 \in ℝ^{m \times m}

the stiffness matrix,

d = 3

the number of dimensions,

N

the number of points, and

m = d N

, the number of degrees of freedom of the structure. This is a simple example of a finite element method. The vectors

𝒇, 𝒖

contain the forces, displacements at all nodes, e.g.,

with each

𝒖_{i} \in ℝ^{d}

. Assembly of the stiffness matrix

𝑲

is carried out by adding contributions from each truss

Such 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 thus equally applicable to physical chemistry, sociology, biology or business management.

1.2Problem data

$•$ Load the data.

The following are typical instructions: open a file, load data, close the file.

∴	ddir = homedir()*"/courses/MATH661/data/EiffelTower";

∴

using MAT

∴

Load the model point coordinates, display the first 3 node coordinates

∴	dfile = matopen(ddir*"/EiffelPoints.mat");

∴	X=read(dfile,"Expression1");

∴	close(dfile);

∴	NX,ndims=size(X);

∴

X[1:3,:]

$[\begin{array}{ccc} 12.563704490661621 & - 12.563704490661621 & 37.89190673828125 \\ - 12.563704490661621 & - 12.563704490661621 & 37.89190673828125 \\ - 14.107048034667969 & - 14.107048034667969 & 37.89190673828125 \end{array}]$ (2)

∴

Load the model edges (links between the nodes) as floating point numbers (Matlab .mat file format)

∴	dfile = matopen(ddir*"/EiffelLines.mat");

∴	flL=read(dfile,"Expression1");

∴	close(dfile);

∴

Display model size

∴	NL,nseg=size(flL); [NX ndims NL nseg]

$[\begin{array}{cccc} 26464 & 3 & 31463 & 2 \end{array}]$ (3)

∴

The model contains $N_{X} = 26464$ nodes (points) each with $d = 3$ coordinates and $N_{L}$ edges (links) between node pairs. Transform the links to integers and display the first 4 edges.

∴	L=zeros(Int64,NL,nseg);

∴	for l=1:NL L[l,:] = floor.(Int64,flL[l,:]) end

∴

L[1:4,:]

$[\begin{array}{cc} 519 & 520 \\ 520 & 521 \\ 521 & 519 \\ 540 & 541 \end{array}]$ (4)

∴

The above states the first edge is from node $i = 519$ to node $j = 520$ .

$•$ Define the system matrix $𝑲$ .

Figure 2. Structure of $𝑲$ at various magnification factors

The data contains $N_{X}$ points in $ℝ^{d}$ , $d = 3$ , and $N_{L}$ line segments defined by 2 endpoints. The total number of degrees of freedom is $m = d N_{X}$ . Of the possible $N_{X}^{2} ≅ 7 \times 10^{8}$ linkages between points in the structure, only $N_{L} = 31463$ are present.

∴	m=3*NX; [NX NX^2 NL m m^2]

$[\begin{array}{ccccc} 26464 & 700343296 & 31463 & 79392 & 6303089664 \end{array}]$ (5)

∴

The matrix $𝑲$ can be formed by loops over the trusses as shown in the algorithm below.

Algorithm $𝑲$ assembly

for $k = 1$ to $N_{L}$

$i = L_{k, 1}$ ; $j = L_{k, 2}$

$𝒍 = 𝑿_{j} - 𝑿_{i}$ ; $𝒍 = 𝒍 / || 𝒍 ||$ ; $𝑷 = 𝒍 𝒍^{T}$

$p = d (i - 1)$ ; $q = d (j - 1)$

for $i_{d} = 1$ to $d$

for $j_{d} = 1$ to $d$

$𝑲 [p + i_{d}, p + j_{d}] = 𝑲 [p + i_{d}, p + j_{d}] - 𝑷 [i_{d}, j_{d}]$

$𝑲 [p + i_{d}, q + j_{d}] = 𝑲 [p + i_{d}, q + j_{d}] + 𝑷 [i_{d}, j_{d}]$

$𝑲 [q + i_{d}, p + j_{d}] = 𝑲 [q + i_{d}, p + j_{d}] + 𝑷 [i_{d}, j_{d}]$

$𝑲 [q + i_{d}, q + j_{d}] = 𝑲 [q + i_{d}, q + j_{d}] - 𝑷 [i_{d}, j_{d}]$

end

$[\begin{array}{l} ⋮ \\ 𝒇_{i} \\ ⋮ \\ 𝒇_{j} \\ ⋮ \end{array}] = [\begin{array}{lll} ⋮ & \dots & ⋮ \\ - 𝑷_{i j} & \dots & 𝑷_{i j} \\ ⋮ & ⋱ & ⋮ \\ 𝑷_{i j} & \dots & - 𝑷_{i j} \\ ⋮ & \dots & ⋮ \end{array}] [\begin{array}{l} ⋮ \\ 𝒖_{i} \\ ⋮ \\ 𝒖_{j} \\ ⋮ \end{array}] .$

Define a function to carry out the assembly process in $d$ dimensions that returns a full matrix $𝑲$ .

∴

function assembleK(d,L,X)
  NL = size(L)[1]; n=maximum(L); m=d*n; K=zeros(m,m)
  for k=1:NL
    i=L[k,1]; j=L[k,2]
    l=X[j,1:d]-X[i,1:d]; l=l/norm(l)
    p=d*(i-1); q=d*(j-1); P=l*l'
    for id=1:d
      for jd=1:d
        K[p+id,p+jd] = K[p+id,p+jd] - P[id,jd]
        K[p+id,q+jd] = K[p+id,q+jd] + P[id,jd]
        K[q+id,p+jd] = K[q+id,p+jd] + P[id,jd]
        K[q+id,q+jd] = K[q+id,q+jd] - P[id,jd]
      end
    end
  end
  return K
end;

∴

XB

$[\begin{array}{c} 0 \\ 1 \end{array}]$ (6)

∴

LB

$[\begin{array}{cc} 1 & 2 \end{array}]$ (7)

∴	assembleK(1,LB,XB)

$[\begin{array}{cc} - 1.0 & 1.0 \\ 1.0 & - 1.0 \end{array}]$ (8)

∴

size(LB)

$[\begin{array}{c} 1 \\ 2 \end{array}]$ (9)

∴

The above has been precomputed for the Eiffel Tower model. It is wasteful, and often impossible due to memory constraints, to store the entire matrix of possible couplings between points $𝑲 \in ℝ^{m \times m}$ ( $m^{2} ≅ 6.3 \times 10^{9}$ ) since most of the matrix would consist of zero entries. Rather, only the nonzero elements corresponding to the linked nodes are stored through what is known as a sparse matrix. There are various techniques to store sparse matrices, one of which is the coordinate format defined by a vector of row indices $i$ , a vector of column indices $j$ and a vector of values $v$ such that the $(i, j)$ element of the matrix $𝑲$ is

k_{i j} = K (i (k), j (k)) = v (k) .

The following instructions load the sparse representation of $𝑲$ .

∴	dfile = matopen(ddir*"/KCOO.mat");

∴	flrows=read(dfile,"rows");

∴	flcols=read(dfile,"cols");

∴	vals=read(dfile,"vals");

∴	close(dfile);

∴	nz=size(flrows)[1]

$370621$

∴	rows=zeros(Int64,nz,1); cols=zeros(Int64,nz,1);

∴	for l=1:nz rows[l] = floor(Int64,flrows[l]) cols[l] = floor(Int64,flcols[l]) end

∴	[minimum(rows) maximum(rows) minimum(cols) maximum(cols)]

$[\begin{array}{cccc} 1 & 79392 & 1 & 79392 \end{array}]$ (10)

∴	using SparseArrays

∴	K=sparse(rows[:,1],cols[:,1],vals[:,1]);

∴

size(K)

$[\begin{array}{c} 79392 \\ 79392 \end{array}]$ (11)

1.3Utility routines

$\circ$ Draw the deformed structure in figure $n_{f}$ , given node coordinates $𝑿 \in ℝ^{N_{X} \times d}$ , displacements $𝑼 \in ℝ^{N_{X} \times d}$ , $m = d N_{X}$ , skipping $k_{s}$ nodes.

∴

function drawXU(X,U,L,ks,nf)
  fig=figure(nf); ax = fig[:add_subplot](1,1,1,projection="3d")
  NX,d = size(X); NL,ns = size(L)
  for k=1:ks:NL
    i=L[k,1]; j=L[k,2]
    x=[X[i,1]+U[i,1] X[j,1]+U[j,1]]'
    y=[X[i,2]+U[i,2] X[j,2]+U[j,2]]'
    z=[X[i,3]+U[i,3] X[j,3]+U[j,3]]'
    ax[:plot](x,y,z,"k")
  end
  axis("equal")
end;

∴

Here are instructions to draw the structure with no node displacement.

∴	U=zeros(NX,3);

∴	drawXU(X,U,L,20,1);

∴

These are instructions to draw the structure with displacements $u$ along the $x_{}$ direction proportional to the $z$ coordinate $u = c z$

∴	U=zeros(NX,3); c=0.25; U[:,1]=c*X[:,3];

∴	drawXU(X,U,L,20,1);

∴

1.4Simplified “A” model

The full Eiffel Tower model is quite complex. It is convenient to build familiarity with modeling binary interactions using the much simpler two-dimensional (

d = 2

) 8-node model shown in Fig. 3

$•$ Define the data for this model and compute $𝑲$ .

∴	x5 = (2.0/3.0)2; x7 = (1.0/3.0)2; [x5 x7]

$[\begin{array}{cc} 1.3333333333333333 & 0.6666666666666666 \end{array}]$ (12)

∴	X=[-2 0; 2 0; -x5 2; 0 2; x5 2; -x7 4; x7 4; 0 6]

$[\begin{array}{cc} - 2.0 & 0.0 \\ 2.0 & 0.0 \\ - 1.3333333333333333 & 2.0 \\ 0.0 & 2.0 \\ 1.3333333333333333 & 2.0 \\ - 0.6666666666666666 & 4.0 \\ 0.6666666666666666 & 4.0 \\ 0.0 & 6.0 \end{array}]$ (13)

∴	L=[ 1 3; 3 6; 6 8; 8 7; 7 5; 5 2; 6 7; 3 4; 4 5; 6 4; 4 7]

$[\begin{array}{cc} 1 & 3 \\ 3 & 6 \\ 6 & 8 \\ 8 & 7 \\ 7 & 5 \\ 5 & 2 \\ 6 & 7 \\ 3 & 4 \\ 4 & 5 \\ 6 & 4 \\ 4 & 7 \end{array}]$ (14)

∴	KA=assembleK(2,L,X);

∴

1.5Enforcing boundary conditions

$\circ$ When structures are attached to their environment, boundary conditions have to be enforced in the linear system $𝒇 = 𝑲 𝒖$ . A common technique is set the diagonal element of $𝑲$ for the node $i$ subjected to boundary conditions to a large value $B$ such that it effectively decouples from the system, and set $f_{i} = B g_{i}$ where $g_{i}$ is the desired boundary value.

∴	maximum(K)

$1592.0704340194352$

∴	f=rand(m,1);

∴	i=500; B=1.0e9; K[i,i]=B; f[i]=B*0.2;

∴	i=501; B=1.0e9; K[i,i]=B; f[i]=B*0.2;

∴	i=502; B=1.0e9; K[i,i]=B; f[i]=B*0.2;

∴	i=1000; B=1.0e9; K[i,i]=B; f[i]=B*0.2;

∴	i=1001; B=1.0e9; K[i,i]=B; f[i]=B*0.2;

∴	i=1002; B=1.0e9; K[i,i]=B; f[i]=B*0.2;

2Track 1 and 2 common problems

3Track 2 task

Carry out a bibliographic search for review papers on reduced order modeling. Choose a paper you find appealing and draft a two-paragraph summary. Compare the approach you find in the literature to problems 3 to 5 above.