Skip to content Skip to navigation

Connexions

You are here: Home » Content » The network wave equation

Navigation

Content Actions
  • Download module PDF
  • Add to ...
    Add the module to:
    • My Favorites
    • A lens
    • An external social bookmarking service
    • My Favorites (What is 'My Favorites'?)
      'My Favorites' is a special kind of lens which you can use to bookmark modules and collections directly in Connexions. 'My Favorites' can only be seen by you, and collections saved in 'My Favorites' can remember the last module you were on. You need a Connexions account to use 'My Favorites'.
    • A lens (What is a lens?)

      Definition of a lens

      Lenses

      A lens is a custom view of Connexions content. You can think of it as a fancy kind of list that will let you see Connexions through the eyes of organizations and people you trust.

      What is in a lens?

      Lens makers point to Connexions materials (modules and collections), creating a guide that includes their own comments and descriptive tags about the content.

      Who can create a lens?

      Any individual Connexions member, a community, or a respected organization.

      What are tags? tag icon

      Tags are descriptors added by lens makers to help label content, attaching a vocabulary that is meaningful in the context of the lens.

    • External bookmarks
  • E-mail the author
  • Rate this module (How does the rating system work?)

    Rating system

    Ratings

    Ratings allow you to judge the quality of modules. If other users have ranked the module then its average rating is displayed below. Ratings are calculated on a scale from one star (Poor) to five stars (Excellent).

    How to rate a module

    Hover over the star that corresponds to the rating you wish to assign. Click on the star to add your rating. Your rating should be based on the quality of the content. You must have an account and be logged in to rate content.

    		(0 ratings)

Recently Viewed

This feature requires Javascript to be enabled.

The network wave equation

Module by: Jesse Chan

Summary: This module introduces an overview of the three-dimensional network wave equation, and discusses numerical solutions and eigenvalue approximations using the finite element method. A Matlab GUI for drawing webs is presented, and eigenvalues from FEM are compared to closed form solutions to the eigenvalues of the one-dimensional network wave equation. As of present, this module contains a rough draft of the material. Links and external files will be added soon.

Note: Your browser may not currently support MathML. See our browser support page for additional details. You can always view the correct math in the PDF version.

Jesse Chan Web Report (1st draft)

The motion of most musical instrument strings can be described by the one dimensional wave equation on an interval x∈[0,ℓ]x∈[0,ℓ], with u(t,0)=u(t,ℓ)=0u(t,0)=u(t,ℓ)=0, where uu 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 c2c2-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 tell 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 to reverse-engineer 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 ∂ 2 u ∂ t 2 = c 2 ∂ 2 u ∂ x 2 , u ( 0 , t ) = u ( ℓ , t ) = 0 (1)

where cc is a constant describing wave speed and ℓℓ is the length of the string. We take Dirichlet boundary conditions at the ends and assume without loss of generality c=1c=1. 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 ∂ u ∂ t = v ∂ v ∂ t = c 2 ∂ 2 u ∂ x 2 (2)

or equivalently, the first order matrix equation

∂ ∂ t u v = 0 I c 2 ∂ 2 ∂ x 2 0 u v ∂ ∂ t u v = 0 I c 2 ∂ 2 ∂ x 2 0 u v (3)

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 0 I ∂ 2 ∂ x 2 0 u v = λ u v , ∂ 2 u ∂ x 2 = λ 2 u (4)

Since only trigonometric functions satisfy both our equation and our boundary conditions, our eigenfunction to take the form u(x)=Asin(λx)+Bcos(λx)u(x)=Asin(λx)+Bcos(λx). Applying our boundary condition at x=0x=0 to u(x)u(x) reveals that B=0B=0. Since we can then set AA as an arbitrary scaling factor, our eigenfunction u(x)u(x) is simply sinλxsinλx. By applying our second boundary condition at x=ℓx=ℓ, we can see that λλ is of the form iπnℓiπnℓ for any nonzero integer nn. We then get the eigenpairs

λ n = i π n ℓ , u n ( x ) = sin ( λ n x ) λ n = i π n ℓ , u n ( x ) = sin ( λ n x ) (5)

These eigenfunctions constitute an infinite-dimensional basis for any solution to the wave equation, with ui(x)ui(x) orthogonal to uj(x)uj(x) for i≠ji≠j with respect to the inner product 〈ui,uj〉≡∫0ℓui(x,t)uj(x,t)dx〈ui,uj〉≡∫0ℓui(x,t)uj(x,t)dx . 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 NN basis functions φi(x)φi(x) that span the space on which the solution is defined. We then calculate the best approximation uN=∑j=1NcjφjuN=∑j=1Ncjφj to the solution from the span of these basis functions via the solution to a matrix equation Ac=fAc=f, where

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 = φ 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 (6)

AA is called the Gramian matrix - a matrix whose ijijth entry is the inner product between the ii and jjth basis functions. After solving for the vector c=[c1,c2,...,cN]Tc=[c1,c2,...,cN]T, we can reconstruct our best approximation to the solution.

We first derive what is called the “weak form" of our PDE. Given a function v(x)v(x) obeying the same boundary conditions as uu, multiply both sides of our wave equation by this function and integrate over the interval [0,ℓ][0,ℓ]

∫ 0 ℓ u t t v d x = ∫ 0 ℓ u x x v d x ∫ 0 ℓ u t t v d x = ∫ 0 ℓ u x x v d x (7)

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 ∫ 0 ℓ u t t v d x = - ∫ 0 ℓ u x v x d x (8)

This form of the wave equation is called the weak form. We now expand uu in the space spanned by our basis functions

u N ( x , t ) = ∑ j = 1 N c j ( t ) φ j ( x ) u N ( x , t ) = ∑ j = 1 N c j ( t ) φ j ( x ) (9)

Let v(x)=φi(x)v(x)=φi(x) for i∈{1,2,...,N}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 ∑ 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 (10)

Note that if we define a new “energy" inner product au,v≡ux,vxau,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 ∑ j = 1 N c i ' ' ( t ) φ i , φ j = ∑ j = 1 N c i ( t ) a φ i , φ j (11)

for i=1,2,...,Ni=1,2,...,N. Thus, we have NN unknowns along with NN linear equations; we can now formulate our problem as the matrix equation

M c ' ' = K c M c ' ' = K c (12)

where MM is the Gramian matrix created using regular inner products, and KK is the Gramian matrix resulting from energy inner products.

Using the finite element method, we choose these basis functions to be piecewise linear “hat" functions. If we partition the space [0,ℓ][0,ℓ] into nn segments of the form [xk-1,xk][xk-1,xk], with x1<x2<...<xNx1<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 φ 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 (13)

for k=1,...,Nk=1,...,N. Since the support of φiφi and φjφj overlap only if |i-j|≤1|i-j|≤1, most of the entries of MM and KK 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][0,1] into these nn segments, with h=1/(N+1)h=1/(N+1) and xk=khxk=kh, then for |i-j|=1|i-j|=1, φi,φi=2h/3φi,φi=2h/3, φi,φj=h/6φi,φj=h/6, aφi,φj=1/haφi,φj=1/h, and aφi,φi=-2/haφi,φi=-2/h. MM and KK are just

M = h 6 4 1 1 4 ⋱ ⋱ ⋱ 1 1 4 , K = 1 h - 2 1 1 - 2 ⋱ ⋱ ⋱ 1 1 - 2 M = h 6 4 1 1 4 ⋱ ⋱ ⋱ 1 1 4 , K = 1 h - 2 1 1 - 2 ⋱ ⋱ ⋱ 1 1 - 2 (14)

We can solve for our coefficients cc by rewriting Mc''=KcMc''=Kc as a system of equations

c ' = d d ' = M - 1 K c c ' = d d ' = M - 1 K c (15)
∂ ∂ t c d = 0 I M - 1 K 0 c d ∂ ∂ t c d = 0 I M - 1 K 0 c d (16)

We can see the relation to the continuous system,

∂ ∂ t u v = 0 I ∂ 2 ∂ x 2 0 u v ∂ ∂ t u v = 0 I ∂ 2 ∂ x 2 0 u v (17)

where ∂2∂x2∂2∂x2 is approximated by M-1KM-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)a(x) is added to the equation

∂ 2 u ∂ t 2 = ∂ 2 u ∂ x 2 - 2 a ( x ) ∂ u ∂ t ∂ 2 u ∂ t 2 = ∂ 2 u ∂ x 2 - 2 a ( x ) ∂ u ∂ t (18)

For the cases we consider here, we shall take a(x)=aa(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 uu

u N = ∑ j = 1 N c j ( t ) φ j ( x ) u N = ∑ j = 1 N c j ( t ) φ j ( x ) (19)

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 ) ∑ 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 ) (20)

Taking an inner product with φkφk for k=1,2,...,Nk=1,2,...,N leads us to the following discretization

M c ' ' = K c - 2 a M c ' M c ' ' = K c - 2 a M c ' (21)

We usually refer to the matrix -2aM-2aM as the damping matrix GG. 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 c ' = d d ' = M - 1 K c - M - 1 G (22)
∂ ∂ t c d = 0 I M - 1 K M - 1 G c d ∂ ∂ t c d = 0 I M - 1 K M - 1 G c d (23)

Should I write about eigenvalues here?

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 iith string in a network of strings be defined on an interval from [0,ℓi][0,ℓi], where ℓiℓ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 [Find REFERENCE ] system of equations for the planar displacement ui(xi,t)ui(xi,t) of the iith string, where xi∈[0,ℓi]xi∈[0,ℓi] . We define the stensor matrix

P i = k i [ ( s i - 1 ) I + v i v i T ] P i = k i [ ( s i - 1 ) I + v i v i T ] (24)

where kiki is stiffness, si>1si>1 is prestress (string tension), and vivi is a unit vector specifying 3-dimensional orientation of the iith string. We characterize network movement by

ρ i I ∂ 2 u i ∂ t 2 = P i ∂ 2 u i ∂ x i 2 ρ i I ∂ 2 u i ∂ t 2 = P i ∂ 2 u i ∂ x i 2 (25)

where ρiρi is the iith strings density. II is the 3-by-3 identity matrix. Our boundary conditions are Dirichlet at endpoints 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ℓk on their respective kkth string. Our Dirichlet conditions can be written as

u 1 ( 0 , t ) = 0 , u k ( ℓ k , t ) = 0 u 1 ( 0 , t ) = 0 , u k ( ℓ k , t ) = 0 (26)

If we define the set SiSi to be the set of integer indices of all strings incident to a joint at the end of the iith string, the force-balance joint conditions connecting strings in the set {i,Si}{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 ) P i ∂ u i ∂ x i ( ℓ i , t ) = ∑ j ∈ S i P j ∂ u j ∂ x j ( 0 , t ) (27)

This network wave equation stensor matrix can also be mathematically derived from the nonlinear model of Antman; the linear, one dimensional wave equation is derived by taking the orientation vector vv to be a standard basis vector.

Figure 1: An example of the notation for the simple tritar case. Here, xixi denotes position on the iith string, and S1={2,3}S1={2,3}.
Figure 1 (tritarpure.png)

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 http://www.tritare.com guitar with Y-shaped strings. 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 u 1 ( 0 , t ) = 0 , u 2 ( ℓ 2 , t ) = 0 , u 3 ( ℓ 3 , t ) = 0 (28)

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 ) P 1 ∂ u 1 ∂ x 1 ( ℓ 1 , t ) = P 2 ∂ u 2 ∂ x 2 ( 0 , t ) + P 3 ∂ u 3 ∂ x 3 ( 0 , t ) (29)

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.

Figure 2: Finite element discretization of a tritar, with a pyramidal hat function φn1φn1 at the joint. r1,r2r1,r2 and r3r3 denote the first, second, and third strings, respectively, which are discretized into n1n1,n2n2, and n3n3 parts, respectively.
Figure 2 (tritarHats.png)

The tritar example case

Let us write out the discretization for the example net in Figure 2. If we take a uniform discretization of each string into n1,n2n1,n2, and n3n3 pieces (with N=n1+n2+n3N=n1+n2+n3), respectively, we can again derive a system of differential equations to describe the evolution of the coefficients ck(t)ck(t) over time. Define the NN basis hat functions as being . Consider first the kkth hat function on string ii, where k≠n1k≠n1. We multiply each side of the network wave equation by the non-joint hat functions φkφ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 ρ 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 (30)

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 ) u N = ∑ j = 1 N c j ( t ) φ j ( x ) (31)

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 ρ 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 (32)

Let LL be the number of connections in our web; L=3L=3 for our tritar. Defining our inner products ·,··,· and a·,·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 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 (33)

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 ρ i I ∑ j = 1 N ∂ 2 c j ( t ) ∂ t 2 φ j , φ k = - P i ∑ j = 1 N c j ( t ) a φ j , φ k (34)

The joint is a different case. Let us our joint hat function be φn1(x)φ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 ( ℓ 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 ( 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 ( 0 , t ) - P 3 ∫ 0 ℓ 3 ∂ u 3 ∂ x 3 ∂ φ n 1 ∂ x 3 d x 3 ρ 1 I ∫ 0 ℓ 1 ∂ 2 u 1 ∂ t 2 ( x 1 , t ) φ n 1 ( x 1 ) d x 1 = P 1 u 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 ( 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 ( 0 , t ) - P 3 ∫ 0 ℓ 3 ∂ u 3 ∂ x 3 ∂ φ n 1 ∂ x 3 d x 3 (35)

after integrating over each string where the joint hat function is nonzero. If we recall that our force balance equation was P1∂x1u1(ℓ1,t)-P2∂x2u2(0,t)-P3∂x3u3(0,t)=0P1∂x1u1(ℓ1,t)-P2∂x2u2(0,t)-P3∂x3u3(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 ∑ 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 (36)

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)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 ∑ 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 (37)

If we define ρ¯=ρ1+ρ2+ρ3ρ¯=ρ1+ρ2+ρ3 and P¯=P1+P2+P3P¯=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 ρ ¯ I ∑ j = 1 N ∂ 2 c j ( t ) ∂ t 2 φ j , φ n 1 = - P ¯ ∑ j = 1 N c j ( t ) a φ j , φ n 1 (38)

Together, equations (Equation 34) and (Equation 38) provide us with a system of equations Mc''=KcMc''=Kc from which to determine our coefficients ck(t)ck(t), where MM and KK 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 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 (39)
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 , 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 , (40)

where ρjkρjk and PjkPjk are linear combinations of ρiρi and PiPi that are determined by the geometry of the network.

If we assume ℓi=1ℓi=1 for i=1,2,3i=1,2,3, then the inner product of two non-joint hat functions is exactly the same as in the one-dimensional case, where, for |j-k|>1|j-k|>1, φj,φj=2h/3φj,φj=2h/3, φj,φk=h/6φj,φk=h/6, aφj,φk=1/haφj,φk=1/h, and aφj,φj=-2/haφj,φj=-2/h. Let us take n1=4n1=4 and n2=n3=3n2=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]v1=[1,0,0], v2=[.5,.5,0]v2=[.5,.5,0], v3=[.5,-.5,0]v3=[.5,-.5,0]. Suppose si=2si=2, ki=1ki=1, ρi=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 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 (41)

and we can assemble MM and KK 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 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 (42)
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 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 (43)

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¯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 hh-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 hh and orientation vv from node to node, and thus reconstruct our PiPi stensors. Knowing the neighbors of each node would allow us to reconstruct the structure of our MM and KK matrices as well. If we examine MM and KK, we can see they are formed out of an NN by NN block grid, where each block is a 3 by 3 matrix. Previously, the index ii differentiated between constants and different legs/connections in our network. In this generalized scheme, we allow ii and jj to reference different nodes in our discretization instead. Thus, a connection between a node ii and jj implies a nonzero entry in the ijijth block, and kijkij, and sijsij refers to the value of the physical constants on the shared support of the hat functions φiφi and φjφj. Utilizing this generalized scheme allows for much more flexibility in terms of our physical constants as well; for example, if the stiffness kk varied as function of xixi, we could capture this by varying our stiffness kk 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 iith node to have xyzxyz position vector pipi. For each node at pipi, there is an associated set of the indices of connected neighbor nodes NiNi and a set CiCi containing the physical constants kij,sijkij,sij pertaining to the connection between nodes ii and j∈Nij∈Ni. Since we can divide through by ρiρi on both sides of the network wave equation, we can assume without loss of generality that the constant ρi=1ρi=1, and that any data carried by the density ρiρi is now contained in kijkij.

Assuming we are given a set of NN nodes, along with the xyxy positions of each node (the zz positions are assumed to be 0, such that the web is planar in the xyxy plane at rest), our first goal is to compute our step sizes hijhij and orientation vectors vijvij for connections between two nodes ii and j∈Nij∈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 NiNi.

Given pipi, NiNi, CiCi, we proceed as follows

  1. for each ii from 1,2,...,N1,2,...,N
    1. for all j∈Nij∈Ni, calculate hij=∥pi-pj∥2hij=∥pi-pj∥2
    2. if the number of elements in any set Ni=1Ni=1
      1. create an “endpoint" node at position p0=2pi-pjp0=2pi-pj, with step size hi0=∥pi-p0∥2hi0=∥pi-p0∥2
      2. set Ni={Ni,0}Ni={Ni,0} to indicate the iith 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 ss. We do this by calculating the maximum distance dmaxdmax between the anchored endpoints of a web, then scaling the positions pipi of every node by the factor s/dmaxs/dmax. Since the absolute positions of the nodes doesn'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φi,φj, aφi,φjaφi,φj, and PijPij. We deal first with constructing the MM matrix, which requires only knowledge of φi,φjφi,φj. Just as in the case of finite element on the single string, most of the basis functions φiφi and φjφ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 MM matrix, we only need to calculate the inner product φi,φjφi,φj for j=ij=i and j∈Nij∈Ni. For the diagonal case j=ij=i, we note that our inner product u,v=∑i=1L∫0ℓiu(xi)v(xi)dxiu,v=∑i=1L∫0ℓiu(xi)v(xi)dxi needs only be calculated on the support of φiφi. For given ii,

φ i , φ i = ∑ j ∈ N i ∫ p i p i + h j φ i 2 ( x ) d x i = ∑ j ∈ N i h i j 3 φ i , φ i = ∑ j ∈ N i ∫ p i p i + h j φ i 2 ( x ) d x i = ∑ j ∈ N i h i j 3 (44)

The last part is a generalization of our inner product for a uniform grid on a single string. For the off-diagonal case j∈Nij∈Ni, j≠ij≠i, the inner product is analogous to the single-string case.

φ i , φ j = h i j 6 φ i , φ j = h i j 6 (45)

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 KK matrix, which requires knowledge of both aφi,φjaφi,φj and our stensor between nodes ii and jj, PijPij. Starting with PijPij, 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 ] , P i j = k i j [ ( s i j - 1 ) I - v i j v i j T ] , (46)

after which we only need aφi,φjaφi,φj for |i-j|≤1|i-j|≤1. For i=ji=j on our main diagonal, the energy inner product is just the sum of the integrals ∫(∂xiφi)2dxi∫(∂xiφi)2dxi evaluated for each leg of the hat function. If each leg lives on a support of length hijhij, the energy inner product is

a φ i , φ i = ∑ j ∈ N i 1 h i j . a φ i , φ i = ∑ j ∈ N i 1 h i j . (47)

For |i-j|=1|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 a φ i , φ j = - 1 h i j (48)

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 proper ties 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 GG is built from the same inner products as our MM 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 aijaij. The ijijth block of GG is then just the ijijth block of MM scaled by aijaij. 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 NN blocks by NN 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.

Thanks goes to Robert Likamwa, who coded much of the 3-dimensional visualization code, and to Jeremy Morell, who was the first to work on larger networks.

Figure 3: A screenshot of the GUI. With the web outline drawn, we can continue to refine our grid until we achieve a desired size.
Figure 3 (GUIscreenshot.png)

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 MM,KK, and GG 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 ∂ ∂ t c d = 0 I M - 1 K M - 1 G c d (49)

Since MM is positive definite, we can compute the matrix-vector product M-1KcM-1Kc and M-1GdM-1Gd using the Cholesky factorization M=CTCM=CTC, where CC is an upper triangular matrix. We can then compute

d ' = C - 1 C - T K c - G d d ' = C - 1 C - T K c - G d (50)

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'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 0 I - K - G u v = λ I 0 0 M u v (51)

This is typically done using eig(A,B) in Matlab to solve this generalized eigenvalue problem, where A,BA,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 Figure 4 to the closed form solutions yields a large discrepancy for larger magnitude eigenvalues.

Figure 4: Comparison of FEM eigenvalue approximations to closed form solutions.
Figure 4 (FEMeigs.png)

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%10-15% of the smallest approximated eigenvalues from FEM tend to be accurate. In order to capitalize as much as possible on this inaccuracy, we try to find the few smallest eigenvalues of our system by finding the largest eigenvalues of the inverse of the system 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. All that's left is to find the inverse of the system.

Suppose we take a Cholesky factorization of M=LL-1M=LL-1. We can then write the system using a similarity transform as

I 0 0 L - 1 0 I - K - G I 0 0 L - T L T