The wave equation

The vibration of a string in one dimension can be understood through the standard wave equation, given by

where uu describes string displacement, cc 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=1c=1. This equation is derived in much more detail in [1]. This second order partial differential equation can likewise be rewritten as a system of two ordinary differential equations in time

or equivalently, the first order matrix equation

Eigenvalues, eigenfunctions, and their significance

We are especially interested in the eigenvalues λλ and associated eigenfunctions of the wave equation, such that

Since only trigonometric functions satisfy both our equation and our boundary conditions, our eigenfunctions 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

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

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

to the solution from the span of these basis functions via the solution to a matrix equation Ac=fAc=f. Recall the definition of our inner product 〈ui,uj〉≡∫0ℓui(x,t)uj(x,t)dx〈ui,uj〉≡∫0ℓui(x,t)uj(x,t)dx. Then, AA is

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 rearrange our PDE into a more flexible form. 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,ℓ]

If we integrate the right hand side by parts and apply Dirichlet boundary conditions, we get

This form of the wave equation is called the equation's “weak form". Notice there is only one derivative with respect to xx on u(x,t)u(x,t) now. We now expand u(x,t)u(x,t) in the space spanned by our basis functions

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

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

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

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 our 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

for k=1,...,Nk=1,...,N.

Figure 2: A hat function centered at x=.6x=.6 with a step size h=.2h=.2.

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

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

We can see the relation to the continuous system,

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

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 ) (21)

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 ) (22)

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 ' (23)

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 (24)

∂ ∂ 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 (25)

Figure 3: Eigenvalues computed from a finite element discretization of a simple string. The progression of the eigenvalues as aa grows is towards the left and towards the real axis

As the damping factor grows from 0 to 6π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 uu decays. For large values of aa, the string becomes overdamped, floating in midair, while for smaller values of aa, the system oscillates before coming to rest. At a=πa=π, the damping is optimal for bringing the string to rest most quickly. This behavior is shown in Figure 4.

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

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 [3] 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

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

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

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

This network wave equation matrix PiPi 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 5: 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

with the force balance equation

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 6: 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.