MATH661 Project 2Large system factorization and model reduction:Time-dependent case

This project extends the study of static binary interactions from Project 1 to time-dependent interactions.

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

$\circ$ 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$ .

$\circ$ 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;

∴

∴

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}]$ (6)

∴

∴	using SparseArrays

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

∴

size(K)

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

∴

K[5,6]

$0.0$

∴	[rows[56] cols[56] vals[56]]

$[\begin{array}{ccc} 22.0 & 14.0 & 1.9999999999999996 \end{array}]$ (8)

∴

K[22,14]

$1.9999999999999996$

∴

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

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;

2Time-dependent case

Newton's second law of motion states that the rate of change of an object's momentum is equal to the externally applied forces. For a single point mass moving along the

x

-axis, this is stated as

where

x (t)

is the point mass position,

\dot{x} (t)

its velocity,

m

its mass (inertia, resistance to change of motion), and

F

the external force. If the mass remains constant (not the case for systems that eject/intake matter such as stars or rockets, neither for motion at speeds close to the speed of light when the resistance to change of motion goes to infinity), the familiar

For binary interaction models (Eiffel tower model, “A” model) a common technique to account for inertia is the so-called “lumped mass” model in which the mass of a truss is allocated as point masses concentrated at the nodes, i.e., the truss end-points. Assembly of contributions from all trusses leads to the system

The advantage of a lumped-mass model is that

𝑴

is a diagonal matrix. Project this system onto the orthonormal basis

𝑩 \in ℝ^{m \times n}

with

𝒗 (t)

the projection of the system trajectory

𝒖 (t) = 𝑩 𝒗 (t)

onto

C (𝑩)

. Obtain

with

𝑵 = 𝑩^{T} 𝑴 𝑩 \in ℝ^{n \times n}

𝑳 = 𝑩^{T} 𝑲 𝑩 \in ℝ^{n \times n}

𝒈 = 𝑩^{T} 𝒇 \in ℝ^{n}

. While

𝑲, 𝑴

are large and require sparse representations,

𝑵, 𝑳

are small and can be stored as full matrices. The above second-order ODE system with

n

equations can be rewritten as a first-order system with

2 n

equations

The notation

𝑵^{- 1} (𝒈 - 𝑳 𝒗)

is meant to signify the solution of the system

𝑵 𝒔 = 𝒈 - 𝑳 𝒗

to find

𝒔

by some numerical algorithm, e.g.,

L U

Q R

decomposition. Since the lumped mass matrix

𝑵 = 𝑩^{T} 𝑴 𝑩

does not vary in time, the decomposition is computed just once, prior to time iteration,

𝑸 𝑹 = 𝑵

. At some given time

t

𝒔 (t) = 𝑵^{- 1} (𝒈 (t) - 𝑳 𝒗 (t))

is found by solving the triangular system

d

dimensions. This corresponds to a force increasing with structure height

x_{d}

, and is a reasonable model of structure interaction with a wind.

1Introduction

1.1Modeling binary interaction in science and engineering

1.2Problem data

1.3Utility routines

1.4Simplified “A” model

1.5Enforcing boundary conditions

2Time-dependent case

3Tasks