A fundamental mode of a complex spiderweb. The motion of vibrating strings (such as musical instrument strings or, in this case, spiderwebs) can be described by the one dimensional wave equation on an interval x∈[0,ℓ], with u(t,0)=u(t,ℓ)=0, where u is the displacement of the string and ℓ is the strings length. The eigenvalues derived from this model progress in a well-known linear fashion, similar to the Western scale, leading to a pleasant sound when the string is plucked. A network of connected strings can be expressed in a similar manner; however, the progression of eigenvalues is much less regular and depends largely on the topology of the network. We examine these eigenvalues and their associated eigenvectors using a finite element discretization of such networks, then compare these results to closed form eigensolutions based on Joachim Von Below's examination of networks of strings in “A Characteristic Equation Associated to an Eigenvalue Problem on c2-Networks", Linear Algebra and its Applications, Volume 71 (1985), p309-325.

Introduction The purpose of the Physics of Strings seminar has traditionally been to study the motion of a vibrating string by analyzing its eigenfunctions and eigenvalues, equivalent to the string's fundamental modes and fundamental frequencies, respectively. The progression of these eigenvalues and eigenvectors tells us a great deal about the string; for example, given eigenvalues of a string, we can determine how quickly its vibrations decay, and whether the frequency of a vibration affects how quickly it's damped. The properties of the string, likewise, can tell us something about the eigenvalues. Physical constants, such as the length of the string, are proportionally related to the eigenvalues. Given data on the vibration of a string, there are also methods for reverse-engineering the eigenvalues of that string. There are several models of a vibrating string, and the most detailed ones can reproduce eigenvalues that accurately match the reverse-engineered string eigenvalues. However, while much research has been done on several models of a single string, the behavior of networks of strings is less well understood. We seek to mathematically model and investigate the motion of networks of strings, specifically by understanding eigenvalues and the corresponding modes of vibration. We study these behaviors within the context of the tritar (a guitar-like instrument based upon a Y-shaped network of 3 strings) and in the vibrations of more complex networks such as spiderwebs.

The wave equation The vibration of a string in one dimension can be understood through the standard wave equation, given by ∂ 2 u ∂ t 2 = c 2 ∂ 2 u ∂ x 2 , u ( 0 , t ) = u ( ℓ , t ) = 0 where u describes string displacement, c is a constant describing wave speed and ℓ is the length of the string. The string is fixed at displacement 0 at the endpoints and assume without loss of generality c=1. This equation is derived in much more detail in . This second order partial differential equation can likewise be rewritten as a system of two ordinary differential equations in time ∂ u ∂ t = v ∂ v ∂ t = c 2 ∂ 2 u ∂ x 2 or equivalently, the first order matrix equation ∂ ∂ t u v = 0 I c 2 ∂ 2 ∂ x 2 0 u v

Eigenvalues, eigenfunctions, and their significance We are especially interested in the eigenvalues λ and associated eigenfunctions of the wave equation, such that 0 I ∂ 2 ∂ x 2 0 u v = λ u v , ∂ 2 u ∂ x 2 = λ 2 u Since only trigonometric functions satisfy both our equation and our boundary conditions, our eigenfunctions take the form u(x)=Asin(λx)+Bcos(λx). Applying our boundary condition at x=0 to u(x) reveals that B=0. Since we can then set A as an arbitrary scaling factor, our eigenfunction u(x) is simply sinλx. By applying our second boundary condition at x=ℓ, we can see that λ is of the form iπnℓ for any nonzero integer n. We then get the eigenpairs λ n = i π n ℓ , u n ( x ) = sin ( λ n x ) These eigenfunctions constitute an infinite-dimensional basis for any solution to the wave equation, with ui(x) orthogonal to uj(x) for i≠j with respect to the inner product 〈 u i , u j 〉 ≡ ∫ 0 ℓ u i ( x , t ) u j ( x , t ) d x Intuitively, these correspond to the fundamental modes of a string - any vibration of the string can be decomposed into a linear combination of the fundamentals. The magnitude of each eigenvalue, likewise, is related to the frequency at which the corresponding fundamental mode vibrates - in other words, each eigenvalue is tied to a note in the progression of the Western scale. As we will see, this linear progression of the eigenvalues is lost when a single string is replaced by a network of strings, leading to more of a dissonant sound when a network is plucked.

Finite element solution method In this report, we use the finite element method to numerically solve for solutions to the wave equation. The idea behind this method is based on picking a finite-dimensional set of N basis functions φi(x) that span the space on which the solution is defined. We then calculate the best approximation u ( x , t ) = ∑ j = 1 N c j ( t ) φ j ( x ) to the solution from the span of these basis functions via the solution to a matrix equation Ac=f. Recall the definition of our inner product 〈ui,uj〉≡∫0ℓui(x,t)uj(x,t)dx. Then, A is A = φ 1 , φ 1 φ 1 , φ 2 ... φ 1 , φ N φ 2 , φ 1 φ 2 , φ 2 ... φ 2 , φ N ⋮ ⋱ φ N , φ 1 φ n , φ 2 ... φ n , φ N , f = f , φ 1 f , φ 2 ⋮ f , φ N A is called the Gramian matrix - a matrix whose ijth entry is the inner product between the i and jth basis functions. After solving for the vector c=[c1,c2,...,cN]T, we can reconstruct our best approximation to the solution. We first rearrange our PDE into a more flexible form. Given a function v(x) obeying the same boundary conditions as u, multiply both sides of our wave equation by this function and integrate over the interval [0,ℓ] ∫ 0 ℓ u t t v d x = ∫ 0 ℓ u x x v d x If we integrate the right hand side by parts and apply Dirichlet boundary conditions, we get ∫ 0 ℓ u t t v d x = - ∫ 0 ℓ u x v x d x This form of the wave equation is called the equation's “weak form". Notice there is only one derivative with respect to x on u(x,t) now. We now expand u(x,t) in the space spanned by our basis functions u N ( x , t ) = ∑ j = 1 N c j ( t ) φ j ( x ) Let v(x)=φi(x) for i∈{1,2,...,N}. Plugging this into the wave equation's weak form, we get the relation ∑ j = 1 N c j ' ' ( t ) ∫ 0 ℓ φ i ( x ) φ j ( x ) d x = ∑ j = 1 N c j ( t ) ∫ 0 ℓ φ i ' ( x ) φ j ' ( x ) d x Note that if we define a new “energy" inner product au,v≡ux,vx, we can then rewrite our whole relation as ∑ j = 1 N c i ' ' ( t ) φ i , φ j = ∑ j = 1 N c i ( t ) a φ i , φ j for i=1,2,...,N. Thus, we have N unknowns along with N linear equations; we can now formulate our problem as the matrix equation M c ' ' = K c where M is the Gramian matrix created using regular inner products, and K is the Gramian matrix resulting from energy inner products. Using the finite element method, we choose our basis functions to be piecewise linear “hat" functions. If we partition the space [0,ℓ] into n segments of the form [xk-1,xk], with x1<x2<...<xN, we can define these hat functions as φ k ( x ) = x - x k - 1 x k - x k - 1 if x ∈ [ x k - 1 , x k ] , x k + 1 - x x k + 1 - x k if x ∈ [ x k , x k + 1 ] , 0 otherwise for k=1,...,N.

A hat function centered at x=.6 with a step size h=.2. Since the support of φi and φj overlap only if |i-j|≤1, most of the entries of M and K are automatically zero. For the rest of the terms, the inner products are easy to compute. If we take a uniform discretization of [0,1] into these n segments, with h=1/(N+1) and xk=kh, then for |i-j|=1, φi,φi=2h/3, φi,φj=h/6, aφi,φj=1/h, and aφi,φi=-2/h. M and K are just M = h 6 4 1 1 4 ⋱ ⋱ ⋱ 1 1 4 , K = 1 h - 2 1 1 - 2 ⋱ ⋱ ⋱ 1 1 - 2 We can solve for our coefficients c by rewriting Mc''=Kc as a system of equations c ' = d d ' = M - 1 K c ∂ ∂ t c d = 0 I M - 1 K 0 c d We can see the relation to the continuous system, ∂ ∂ t u v = 0 I ∂ 2 ∂ x 2 0 u v where ∂2∂x2 is approximated by M-1K. With this discretization, we can numerically calculate the time solution of the wave equation given some initial condition, as well as approximate the eigenvalues λ.

Damping A closely related equation is the wave equation with viscous damping (resulting from a viscous medium in which the string vibrates, i.e. air). To simulate this effect, a velocity-dependent damping function a(x) is added to the equation ∂ 2 u ∂ t 2 = ∂ 2 u ∂ x 2 - 2 a ( x ) ∂ u ∂ t For the cases we consider here, we shall take a(x)=a, some constant. Thankfully, the finite element discretization of this equation doesn't involve much new work; all we do is reuse some of our calculations. If we make the substitution for u u N = ∑ j = 1 N c j ( t ) φ j ( x ) we get ∑ j = 1 N c j ' ' ( t ) φ j ( x ) = ∑ j = 1 N c j ( t ) φ j ' ' ( x ) - 2 a ∑ j = 1 N c j ' ( t ) φ j ( x ) Taking an inner product with φk for k=1,2,...,N leads us to the following discretization M c ' ' = K c - 2 a M c ' We usually refer to the matrix -2aM as the damping matrix G. Again, we can solve this by writing it out as a system of ordinary differential equations c ' = d d ' = M - 1 K c - M - 1 G ∂ ∂ t c d = 0 I M - 1 K M - 1 G c d

Eigenvalues computed from a finite element discretization of a simple string. The progression of the eigenvalues as a grows is towards the left and towards the real axis As the damping factor grows from 0 to 6π, the eigenvalues shift further left in the complex plane. At factors of π, the two smallest magnitude eigenvalues become completely real, with one moving left and one moving right in the complex plane. The physical significance of this lies in the fact that the real part of the eigenvalue furthest in the right half plane is determines proportionally how quickly the displacement u decays. For large values of a, the string becomes overdamped, floating in midair, while for smaller values of a, the system oscillates before coming to rest. At a=π, the damping is optimal for bringing the string to rest most quickly. This behavior is shown in .

Displacement of the midpoint of a string for different damping terms. Notice for damping factor a=π, the displacement reaches a steady state fastest.

Networks of strings Unlike our simple one dimensional case, it is much more difficult to determine the closed form eigenvalues and eigenfunctions of a network of strings. To this end, we apply the finite element method to numerically simulate the behavior of a network wave equation.

Network wave equation Let the ith string in a network of strings be defined on an interval from [0,ℓi], where ℓi is the length of that particular string. To generalize the wave equation to a network of strings in three dimensions, we reference Schmidt's system of equations for the planar displacement ui(xi,t) of the ith string, where xi∈[0,ℓi] . We define the "stensor" matrix P i = k i [ ( s i - 1 ) I + v i v i T ] where ki is stiffness, si>1 is prestress (string tension), and vi is a unit vector specifying 3-dimensional orientation of the ith string. We characterize network movement by ρ i I ∂ 2 u i ∂ t 2 = P i ∂ 2 u i ∂ x i 2 where ρi is the ith strings density. I is the 3-by-3 identity matrix. Our boundary conditions are Dirichlet at endpoints (displacement is fixed at 0) and a condition enforcing force balance laws and connectivity of each leg at the joint. We define an end of the first string to have position 0, and for the other endpoints, we consider them to be at position ℓk on their respective kth string. Our Dirichlet conditions can be written as u 1 ( 0 , t ) = 0 , u k ( ℓ k , t ) = 0 If we define the set Si to be the set of integer indices of all strings incident to a joint at the end of the ith string, the force-balance joint conditions connecting strings in the set {i,Si} can be described by P i ∂ u i ∂ x i ( ℓ i , t ) = ∑ j ∈ S i P j ∂ u j ∂ x j ( 0 , t ) This network wave equation matrix Pi can also be mathematically derived from the nonlinear model of Antman; the linear, one dimensional wave equation is derived by taking the orientation vector v to be a standard basis vector.

An example of the notation for the simple tritar case. The network wave equation is much more tractable for a concrete example. We begin by covering the network wave equation for the simplest net - a Y-shaped net called a “tritar", in honor of the guitar with Y-shaped strings (see http://www.tritare.com). For our simple case, then, we have the boundary conditions u 1 ( 0 , t ) = 0 , u 2 ( ℓ 2 , t ) = 0 , u 3 ( ℓ 3 , t ) = 0 with the force balance equation P 1 ∂ u 1 ∂ x 1 ( ℓ 1 , t ) = P 2 ∂ u 2 ∂ x 2 ( 0 , t ) + P 3 ∂ u 3 ∂ x 3 ( 0 , t ) We will investigate this example further using a discretization of the network.

Finite element discretization of the network wave equation To model behavior and structure of a continuous network, we discretize and solve our equations using the finite element method. For the most part, applying FEM to our network model is the same as applying it to a simple string - the hat functions overlap and form a basis for the structure of each leg. The exception is at a joint, which has a new type of hat function, with its support spanning a small section of each string connected at that joint.

Finite element discretization of a tritar, with a pyramidal hat function φn1 at the joint. r1,r2 and r3 denote the first, second, and third strings, respectively, which are discretized into n1,n2, and n3 parts, respectively.

The tritar example case Let us write out the discretization for the example net in . If we take a uniform discretization of each string into n1,n2, and n3 pieces (with N=n1+n2+n3), respectively, we can again derive a system of differential equations to describe the evolution of the coefficients ck(t) over time. Define the N basis hat functions as being . Consider first the kth hat function on string i, where k≠n1. We multiply each side of the network wave equation by the non-joint hat functions φk and integrate over the support of that function. After integration by parts, we have the relation ρ i I ∫ 0 ℓ i ∂ 2 u i ∂ t 2 ( x i , t ) φ k ( x i ) d x i = - P i ∫ 0 ℓ i ∂ u i ∂ x i ∂ φ k ∂ x i d x i analagous to the one dimensional finite element discretization of a string. If we substitute in our approximation from the basis of hat functions u N = ∑ j = 1 N c j ( t ) φ j ( x ) we arrive at the relation ρ i I ∑ j = 1 N ∂ 2 c j ( t ) ∂ t 2 ∫ 0 ℓ i φ j ( x i ) φ k ( x i ) d x i = - P i ∑ j = 1 N c j ( t ) ∫ 0 ℓ i ∂ φ j ∂ x i ∂ φ k ∂ x i d x i Let L be the number of connections in our web; L=3 for our tritar. Defining our inner products ·,· and a·,· as u , v = ∑ i = 1 L ∫ 0 ℓ i u ( x i ) v ( x i ) d x i , a u , v = ∑ i = 1 L ∫ 0 ℓ i ∂ u ( x i ) ∂ x i ∂ v ( x i ) ∂ x i d x i we see these inner products behave much like the simple string inner products on the topology our network. This gives the relation ρ i I ∑ j = 1 N ∂ 2 c j ( t ) ∂ t 2 φ j , φ k = - P i ∑ j = 1 N c j ( t ) a φ j , φ k The joint is a different case. Let us our joint hat function be φn1(x). Then, since integration by parts moves a derivative from one function to another with the addition of a boundary value, we get ρ 1 I ∫ 0 ℓ 1 ∂ 2 u 1 ∂ t 2 ( x 1 , t ) φ n 1 ( x 1 ) d x 1 = P 1 ∂ u 1 ∂ x 1 ( ℓ 1 , t ) - P 1 ∫ 0 ℓ 1 ∂ u 1 ∂ x 1 ∂ φ n 1 ∂ x 1 d x 1 ρ 2 I ∫ 0 ℓ 2 ∂ 2 u 2 ∂ t 2 ( x 2 , t ) φ n 1 ( x 2 ) d x 2 = - P 2 ∂ u 2 ∂ x 2 ( 0 , t ) - P 2 ∫ 0 ℓ 2 ∂ u 2 ∂ x 2 ∂ φ n 1 ∂ x 2 d x 2 ρ 3 I ∫ 0 ℓ 3 ∂ 2 u 3 ∂ t 2 ( x 3 , t ) φ n 1 ( x 3 ) d x 3 = - P 3 ∂ u 3 ∂ x 3 ( 0 , t ) - P 3 ∫ 0 ℓ 3 ∂ u 3 ∂ x 3 ∂ φ n 1 ∂ x 3 d x 3 after integrating over each string where the joint hat function is nonzero. If we recall that our force balance equation was P 1 ∂ u 1 ∂ x 1 ( ℓ 1 , t ) - P 2 ∂ u 2 ∂ x 2 ( 0 , t ) - P 3 ∂ u 3 ∂ x 3 ( 0 , t ) = 0 , however, we can sum these equations together to achieve the relation ∑ i = 1 3 ρ i I ∫ 0 ℓ i ∂ 2 u i ∂ t 2 ( x i , t ) φ n 1 ( x i ) d x i = - ∑ i = 1 3 P i ∫ 0 ℓ i ∂ u i ∂ x i ∂ φ n 1 ∂ x i d x i Conveniently, the force balance equation allows us to generalize this condition to joints with multiple legs as well. Next, substituting in uN=∑j=1Ncj(t)φj(x), we get ∑ i = 1 3 ρ i I ∑ j = 1 N ∂ 2 c j ( t ) ∂ t 2 ∫ 0 ℓ i φ j ( x i ) φ n 1 ( x i ) d x i = - ∑ i = 1 3 P i ∑ j = 1 N c j ( t ) ∫ 0 ℓ i ∂ φ j ∂ x i ∂ φ n 1 ∂ x i d x i If we define ρ¯=ρ1+ρ2+ρ3 and P¯=P1+P2+P3, we are then left with the relation ρ ¯ I ∑ j = 1 N ∂ 2 c j ( t ) ∂ t 2 φ j , φ n 1 = - P ¯ ∑ j = 1 N c j ( t ) a φ j , φ n 1 Together, equations () and () provide us with a system of equations Mc''=Kc from which to determine our coefficients ck(t), where M and K are M = ρ 11 I φ 1 , φ 1 ρ 12 I φ 1 , φ 2 ... ρ 1 N I φ 1 , φ N ρ 21 I φ 2 , φ 1 ρ 22 I φ 2 , φ 2 ... ρ 2 N I φ 2 , φ N ⋮ ⋱ ρ N 1 I φ N , φ 1 ρ N 2 I φ N , φ 2 ... ρ N N I φ N , φ N K = P 11 a φ 1 , φ 1 P 12 a φ 1 , φ 2 ... P 1 N a φ 1 , φ N P 21 a φ 2 , φ 1 P 22 a φ 2 , φ 2 ... P 2 N a φ 2 , φ N ⋮ ⋱ P N 1 a φ N , φ 1 P N 2 a φ N , φ 2 ... P N N a φ N , φ N , where ρjk and Pjk are linear combinations of ρi and Pi that are determined by the geometry of the network. If we assume ℓi=1 for i=1,2,3, then the inner product of two non-joint hat functions is exactly the same as in the one-dimensional case, where φ i , φ j = 2 h / 3 , if i = j ; h / 3 , if | i - j | = 1 ; 0 , otherwise. and a φ i , φ j = - 2 / h , if i = j ; 1 / h , if | i - j | = 1 ; 0 , otherwise. Let us take n1=4 and n2=n3=3. For this network , if we assume all the legs of the tritar lie at equal angles from each other, we can define the orientation v1=[1,0,0], v2=[.5,.5,0], v3=[.5,-.5,0]. Suppose si=2, ki=1, ρi=1. Then, P 1 = 2 0 0 0 1 0 0 0 1 , P 2 = 5 / 4 1 / 4 0 1 / 4 5 / 4 0 0 0 1 , P 3 = 5 / 4 - 1 / 4 0 - 1 / 4 5 / 4 0 0 0 1 and we can assemble M and K as follows M = h 6 4 ρ 1 I ρ 1 I 0 0 0 0 0 0 0 0 ρ 1 I 4 ρ 1 I ρ 1 I 0 0 0 0 0 0 0 0 ρ 1 I 4 ρ 1 I ρ 2 I 0 0 0 0 0 0 0 0 ρ 2 I 2 ρ ¯ I ρ 2 I 0 0 ρ 3 I 0 0 0 0 0 ρ 2 I 4 ρ 2 I ρ 2 I 0 0 0 0 0 0 0 0 ρ 2 I 4 ρ 2 I ρ 2 I 0 0 0 0 0 0 0 0 ρ 2 I 4 ρ 2 I 0 0 0 0 0 0 ρ 3 I 0 0 0 4 ρ 3 I ρ 3 I 0 0 0 0 0 0 0 0 ρ 3 I 4 ρ 3 I ρ 3 I 0 0 0 0 0 0 0 0 ρ 3 I 4 ρ 3 I K = 1 h - 2 P 1 P 1 0 0 0 0 0 0 0 0 P 1 - 2 P 1 P 1 0 0 0 0 0 0 0 0 P 1 - 2 P 1 P 2 0 0 0 0 0 0 0 0 P 2 - P ¯ P 2 0 0 P 3 0 0 0 0 0 P 2 - 2 P 2 - P 2 0 0 0 0 0 0 0 0 P 2 - 2 P 2 P 2 0 0 0 0 0 0 0 0 P 2 - 2 P 2 0 0 0 0 0 0 P 3 0 0 0 - 2 P 3 P 3 0 0 0 0 0 0 0 0 P 3 - 2 P 3 P 3 0 0 0 0 0 0 0 0 P 3 - 2 P 3 We can reverse engineer some of the geometry of our network from examination of these matrices - notice that each leg has 3 blocks assigned to it, corresponding to the 3 non-joint hat functions on each string. The far off-diagonal terms capture the connection of the first string to the third string, and the presence of ρ¯ and P¯ on the diagonal stems from the inner product of the joint hat function with the hat functions on each of the strings.

Generalized numbering scheme and -adaptivity Unfortunately, for larger and more complex webs, writing the system out by hand becomes far too tedious. We seek a more systematic and flexible way of producing our finite element discretizations. We should note two things about finite element discretizations. First, if we stay consistent, a reordering of the nodes does not affect our discretization, though it may change the structure of our matrix. Secondly, our hat functions are not required to be either uniform or symmetric - they can vary in width depending on index, and one side can have a different width than another. This idea is known as h-adaptivity; advanced finite element methods tend to adapt their discretizations by using error estimates from iteration to iteration to pinpoint areas where a coarse discretization should be refined to allow for greater accuracy. Knowing this, it is possible to produce a generalized finite element discretization of a web given only physical constants, a set of nodal points and each point's neighbors. Given this, we can calculate the step size h and orientation v from node to node, and thus reconstruct our Pi matrices. Knowing the neighbors of each node would allow us to reconstruct the structure of our M and K matrices as well. If we examine M and K, we can see they are formed out of an N by N block grid, where each block is a 3 by 3 matrix. Previously, the index i differentiated between constants and different legs/connections in our network. In this generalized scheme, we allow i and j to reference different nodes in our discretization instead. Thus, a connection between a node i and j implies a nonzero entry in the ijth block, and kij, and sij refers to the value of the physical constants on the shared support of the hat functions φi and φj. Utilizing this generalized scheme allows for much more flexibility in terms of our physical constants as well; for example, if the stiffness k varied as function of xi, we could capture this by varying our stiffness k from node to node. We go into more detail on this in the next section. Many of the concepts from the single-string case carry over to networks.

Assembling the and matrices We begin by describing the notation of the information represented by our data structures. We denote the ith node to have xyz position vector pi. For each node at pi, there is an associated set of the indices of connected neighbor nodes Ni and a set Ci containing the physical constants kij,sij pertaining to the connection between nodes i and j∈Ni. Since we can divide through by ρi on both sides of the network wave equation, we can assume without loss of generality that the constant ρi=1, and that any data carried by the density ρi is now contained in kij. Assuming we are given a set of N nodes, along with the xy positions of each node (the z positions are assumed to be 0, such that the web is planar in the xy plane at rest), our first goal is to compute our step sizes hij and orientation vectors vij for connections between two nodes i and j∈Ni. To account for Dirichlet boundary conditions, we also create an anchored node for each endpoint. In this implementation, if a node has only one neighbor, we assume it is connected to a pinned endpoint whose position is in the opposite direction but the same distance away as the only neighbor (this is required to calculate an inner product). For a node connected to an endpoint, we append the index 0 to Ni. Given pi, Ni, Ci, we proceed as follows for each i from 1,2,...,Nfor all j∈Ni, calculate hij=∥pi-pj∥2if the number of elements in any set Ni=1create an “endpoint" node at position p0=2pi-pj, with step size hi0=∥pi-p0∥2set Ni={Ni,0} to indicate the ith node is connected to an anchor In practice, we normalize the positions of our nodes such that the web lies within a box of a desired arbitrary size s. We do this by calculating the maximum distance dmax between the anchored endpoints of a web, then scaling the positions pi of every node by the factor s/dmax. Since the absolute positions of the nodes don't affect our discretization, we don't need to worry about subtracting off the centroid of all the node positions. With all our variables now in place, we can now proceed to the actual construction of our discretization matrices. This requires knowing φi,φj, aφi,φj, and Pij. We deal first with constructing the M matrix, which requires only knowledge of φi,φj. Just as in the case of finite element on the single string, most of the basis functions φi and φj don't share support and their inner products are zero. However, in addition to calculating inner products of regular hat functions, we need to compute the inner products with joint, generalized and nonsymmetric hat functions as well. Starting with our M matrix, we only need to calculate the inner product φi,φj for j=i and j∈Ni. For the diagonal case j=i, we note that our inner product u,v=∑i=1L∫0ℓiu(xi)v(xi)dxi needs only be calculated on the support of φi. For given i, φ i , φ i = ∑ j ∈ N i ∫ p i p i + h j φ i 2 ( x ) d x i = ∑ j ∈ N i h i j 3 The last part is a generalization of our inner product for a uniform grid on a single string. For the off-diagonal case j∈Ni, j≠i, the inner product is analogous to the single-string case. φ i , φ j = h i j 6 since two different hat function can overlap on at most one leg (otherwise two legs of a hat function could cover the same support). Next, we can create our K matrix, which requires knowledge of both aφi,φj and our Pi matrix between nodes i and j, Pij. Starting with Pij, if we have done our bookkeeping correctly, we have all the variables we need to compute P i j = k i j [ ( s i j - 1 ) I - v i j v i j T ] , after which we only need aφi,φj for |i-j|≤1. For i=j on our main diagonal, the energy inner product is just the sum of the integrals ∫(φi''(xi))2dxi evaluated for each leg of the hat function. If each leg lives on a support of length hij, the energy inner product is a φ i , φ i = ∑ j ∈ N i 1 h i j . For |i-j|=1, two hat functions can share support on at most one leg, so our energy inner product is a φ i , φ j = - 1 h i j which again is analogous to our single-string case.

Damping The case of the damped network wave equation is worth examining as well, especially in the mathematical modeling of a spider's web. The material properties of spiderwebs also make it ideal for simulation via the second order wave equation. These include minimal torsion (twisting) in vibrations, low stiffness, no hysteresis under small strains, and a loss of energy primarily through aerodynamic damping. The wave equation assumes negligible torsion and low stiffness, is meant to model string movement specifically under small strains, and is easy to add a constant aerodynamic/viscous damping term to. Since the structure of our damping matrix G is built from the same inner products as our M matrix; the only difference is that we now have to keep track of one more constant, the damping coefficient on a connection between two nodal points aij. The ijth block of G is then just the ijth block of M scaled by aij. This allows us to again vary damping from connection to connection, which proves useful in the simulation of spider webs, since the radial and axial fibers of a spiderweb are often subject to different levels of damping.

Matlab GUI With this last bit of information, we know each block entry of our N blocks by N blocks discretization matrices, and can construct a finite element discretization for a web given only a list of nodes, their positions, and their connectivity. To implement this in an accessible way, a Matlab “point-and-click" GUI was developed to allow users to trace and experiment with their own webs through numerical simulations of web motion and analysis of the eigenvalues and fundamental modes.

A screenshot of the GUI. With the web outline drawn, we can continue to refine our grid until we achieve a desired size.

Setting up the web Using a GUI to wrap around our framework which allows the user to point and click to place nodes down, then to click from one node to another to specify the connection pattern. Endpoints (where the nodes are pinned down, enforcing Dirichlet boundary conditions) are assumed to be nodes with only one neighbor (i.e., not a link in a chain). Once the initial pattern is set, the user can change the discretization fineness as desired, as well as rescale the size of the web to a larger or small grid. When the user is done, the positions and connection pattern of the nodes can be used to create a finite element discretization of the network of strings. While code behind the discretization of the web remains the main engine driving the mathematics of this model, the GUI has probably been where most of the work has gone, making the creation of webs accessible to anyone, although it's still possible to create a web simply by loading a data file.

Solving for the solution with FEM Along with calculating eigenvalues, we can solve a system of differential equations in order to solve our second order PDE once given our M,K, and G matrices. Since the wave equation with damping can be written as an ODE system ∂ ∂ t c d = 0 I M - 1 K M - 1 G c d Since M is positive definite, we can compute the matrix-vector product M-1Kc and M-1Gd using the Cholesky factorization M=CTC, where C is an upper triangular matrix. We can then compute d ' = C - 1 C - T K c - G d and use a numerical solver. We feed into Matlab's ode45 to compute our solution over time given any initial condition. To simulate a smooth initial ripple, I coded in a 3-point Gaussian as the initial displacement to a single node when the user chooses to pluck the web at that specific node. Plucking multiple nodes sums the displacement up over each node, so that any overlap of the Gaussian initial condition between two nodes is accounted for.

Solving the eigenvalue problem with FEM Once we're given the mass and stiffness matrices, it's easy to numerically solve for the eigenvalues. Writing our second order system with damping (Mc''=Kc-Gc') as a linear system, we can solve for our eigenvalues by solving the system 0 I - K - G u v = λ I 0 0 M u v This is typically done using eig(A,B) in Matlab to solve this generalized eigenvalue problem, where A,B are, respectively, the left and right hand matrices in the above equation. However, a comparison of the FEM eigenvalues of a simple string in to the closed form solutions yields a large discrepancy for larger magnitude eigenvalues.

Comparison of FEM eigenvalue approximations to closed form solutions. This behavior is typical of finite element discretizations; given the non-smooth nature of our basis functions, it is difficult to accurately approximate high-frequency eigenfunctions, and typically, only about 10-15% of the smallest approximated eigenvalues from FEM tend to be accurate. In order to compensate for this inaccuracy as much as possible, we try to find the few smallest eigenvalues of our system by finding the largest eigenvalues of the inverse of the system (which is arrived at from algebraic manipulation and similarity transforms) and then inverting those values. The eigs package in Matlab is suited perfectly for that; eigs can quickly compute the few smallest or largest eigenvalues of a large matrix where traditional eigensolvers typically fail due to high algorithmic complexity costs.

Other approaches Since the finite element method does a poor job of representing higher frequency modes, work has been done by Dr. Embree and Jeremy Morell in implementing a similar generalized spectral discretization using Chebyshev polynomials, which proved much more accurate in describing eigenmodes corresponding to larger magnitude eigenvalues.

Interpretation of results In this section, we examine the eigenvalues of networks computed from our finite element discretization. While analysis of these values turns out to be difficult, we examine closed form solutions for a similar eigenvalue problem on networks from a paper by Joachim Von Below .

Numerical results Unlike our original wave equation uxx=utt, our networks of strings are allowed three-dimensional freedom of motion. We can apply our network wave equation to a single string u x x = P u t t where v is the unit vector specifying orientation of our string and P=k[(s-1)I-vvT]. Assume without loss of generality that v=k→=[0,0,1]T and that k=1, s=2. Then P is simply a diagonal matrix and our equation uxx=Putt becomes u i u j u k t t = 1 0 0 0 1 0 0 0 2 u i u j u k x x where ui,uj,uk are the the displacements of our string in the i,j,k directions. Since each of the equations is independent of the others, we can solve for the eigenvalues and eigenfunctions of each one-dimensional wave equation separately λ i 2 u i = ∂ 2 u i ∂ x 2 λ j 2 u j = ∂ 2 u j ∂ x 2 2 λ k 2 u k = ∂ 2 u k ∂ x 2 Then, if λ is an eigenvalue of our one-dimensional wave equation, the eigenvalues λi,λj and λk for the three dimensional wave equation are λ i = λ λ j = λ λ k = 2 λ We can see this captured in Figures and - the eigenvalues of our discretization are the interleaved eigenvalues of three one dimensional wave equations. For general orientations, the result is the same.

Computed eigenvalues of the three-dimensional wave equation.

Closed-form eigenvalues of the three-dimensional wave equation. We can still trace out the linear progression of the eigenvalues here. However, the eigenvalues of a network of strings turn out to be far more interesting and unpredictable. presents the first few eigenvalues of a Y-shaped network of strings, similar to our tritar mentioned previously. Even among simple webs such as the tritar, the pattern of the progression of eigenvalues is not easily deduced.

Eigenvalues of a tritar

Eigenvalues of a more complex network. We can observe a few parts at which the eigenvalue behavior mimics the three dimensional single string. At values around 1.5i and 3i, there are double eigenvalues reminiscent of our double eigenvalues in , but the pattern of the rest of the eigenvalues is much less coherent.

Mode 5

Mode 10

Mode 11

Mode 12 Figures to are FEM calculations of several eigenmodes of a more complex network. Note that as the number of legs and connections increase, the number of degrees of freedom for the movement of each leg (and thus the number of possible eigenmodes of the network) should increase as well. The eigenvalues for the more complex network exhibit a similarly nonlinear pattern in the progression of the eigenvalues of the tritar.

Mode 5

Mode 9

Mode 10

Mode 12

Closed form solutions To better understand the nonlinear progression of our eigenvalues, we seek out closed form solutions for the eigenvalues of networks. Joachim Von Below's provides one such solution in his examination of networks of strings in “A Characteristic Equation Associated to an Eigenvalue Problem on c2-Networks", Linear Algebra and its Applications, Volume 71 (1985), p309-325. While his paper focuses on the heat equation on networks, given Neumann boundary conditions instead of Dirichlet, the eigenvalues of the heat equation are given to be constant μ such that ∂ 2 u ∂ x 2 = μ u Since μ has to be less than or equal to zero in order for a solution to the heat equation, the eigenvalues of the heat equation are the squares of the eigenvalues of the wave equation ∂ 2 u ∂ x 2 = λ 2 u Therefore, taking the square roots of the eigenvalues μ, we can recover the eigenvalues of the wave equation on our network. For the rest of our stated examples, the eigenvalues mentioned are assumed to be eigenvalues of the wave equation. Von Below studies the network heat equation with Neumann conditions at the endpoints ∂ u j ∂ t = a j ∂ 2 u j ∂ x j 2 where uj is the temperature of the jth edge, aj is the positive diffusion coefficient of the jth edge, and xj∈[0,ℓj] is position on the jth edge, which has length ℓj. The notation is the same as the notation with which the network wave equation was described. By defining simply a set of nodes and connections between them, Von Below calculates the eigenvalues and eigenmodes of a network with the help of spectral graph theory (a study of the connections between a graph and the eigenvalues of its adjacency matrix). The eigenvalues of the heat equation on networks of strings progress according to several different patterns, some of which are common to every network, some of which are network-specific and dependent on the eigenvalues of the network graph's row-normalized adjacency matrix. For nonspecific parameter values, the best we can hope to do is to solve a transcendental equation for our eigenvalues. For specific parameter values, such as when aj=ℓj2, we can solve analytically for both the eigenfunctions and eigenvalues of the heat equation, breaking both down into cases. There are, in general, 5 total classes into which eigenpairs of the heat equation on a network can be organized, with each network having at most 4 classes applicable to it.

The tritar's eigenvalues sorted by cases

The tritar's eigenvalues sorted by magnitude. These eigenvalues of the tritar are calculated via Von Below's methods. For each eigenvalue case, there exist an infinite number of eigenvalues For each one of these classes mentioned, there is an infinite number of eigenvalues associated with it. Organizing the progression by k instead of by magnitude gives more insight into the nonlinear progression of our network eigenvalues - there is a clear bifurcating/splitting patten of the eigenvalues as k increases. As the progression of eigenvalues splits into two separate patterns, each pattern is linear with different rates. The fundamental difference here is that the eigenvalues for the one dimensional transverse network wave equation progress with two different rates. While there are two distinct sets of eigenvalues for our three dimensional string, to have eigenvalues progressing at different rates for one dimensional transverse string motion is atypical. This isn't a complete explanation for the irregular progression of the eigenvalues for the three dimensional wave equation; however, it does offer some clues as to the nonlinear ordering of those eigenvalues.

Acknowledgments Thanks to Jeremy Morell, whose work I built upon, as well as Robert Likamwa, who helped nail down a lot of details concerning visualization and graphics. Thanks to Dr. Cox and Dr. Embree for their guidance their trust that I could actually do something with the project they gave me. This work was partially supported by NSF DMS Grant 0240058