The network wave equation

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

Since $M$ is positive definite, we can compute the matrix-vector product $M^{-1}Kc$ and $M^{-1}Gd$ using the Cholesky factorization $M=C^{T}C$ , where $C$ is an upper triangular matrix. We can then compute

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 ( $M{c}^{\text{'}\text{'}}=Kc-G{c}^{\text{'}}$ ) as a linear system, we can solve for our eigenvalues by solving the system

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.

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; eigscan 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.