The Vibrating String
Calculating Eigenvalues and Eigenfunctions
with Mathematica

Initialization

Needs["Graphics`Colors`"];
Off[General::spell];
Off[General::spell1];

Consider the wave equation as a model for a vibrating string:

D[p[x] D[u[x,t], x], x] == rho[x] D[u[x,t], t, t]

Where u[x,t] is the wave functions, p[x] is the modulus of elasticity times the cross sectional area, and rho is the density per unit length.

Suppose that coordinates are normalized so that the length of the string is 1, and suppose that the ends are clamped.  Then the boundary conditions are given by

u[0,t] == u[1,t] == 0

The problem to we solve here is to find the solution with a motionless initial profile (string plucked and then released)

u[x,0] == f[x],   Derivative[0,1][u][x,0] == 0.

This problem can be solved by separation of variables.  Suppose that u[x,t] = v[x]w[t], then, the wave equation  is equivalent to

D[p[x] D[v[x], x], x]/(rho[x] v[x]) == D[w[t], t, t]/w[t]

If both sides are equal to a constant, say -L^2, then the variables separate and the solution form v[x]w[t] is valid.  The values L will be referred to as the eigenvalues (or the spectrum).  Because of the nature of the problem, the set of functions v[k][x] for all of the infinite number of values L[k] is complete, so any reasonable function f[x] which satisfies the boundary conditions (f[0] == f[1] == 0) can be expressed by

f[x] = Sum[a[k] v[k][x],{k,1,Infinity}]

And thus, the solution can be expressed by

u[x,t] = Sum[a[k] v[k][x] w[k][t],{k,1,Infinity}]

The form of the functions w[k][t] is simple, since w[k][x] satisfies

D[w[k][t], t, t] + L[k]^2 w[L[k]] == 0,  w[k][0] == 1,  w[k]'[0] == 0]

So   w[k][t] = Cos[L[k] t]

The problem of finding the eigenvalues L[k] and the eigenfunctions v[k][x] is more difficult, but provides and excellent example of using a variety of Mathematica's numerical functions.

First start with a familar problem which can be solved symbolically as well; the unifrom string

The uniform string.

Suppose that physical units are scaled so that p[x] == rho[x] == 1.  Then you need to solve the equation with the boundary conditions.

DSolve[{
	D[vv[x],x,x] + LL^2 vv[x] == 0,
	vv[0] == 0,vv[1] == 0},
	vv[x],x]

DSolve does not necessarily find all possible solutions, so we will have to reformulate the problem.  Since the equation is linear, the size of the first derivative at zero does not matter.  So to look for non-zero solutions, it makes sense to try:

dsol = DSolve[{
	D[vv[x],x,x] + LL^2 vv[x] == 0,
	vv[0] == 0,vv'[0] == 1},
	vv,x][[1]]

Now you need to find the values of LL for which vv[1] == 0.  The Solve[] function will not do this for you because it uses the ArcSin[] function which only preduces the zero solution for Sin[LL] == 0.  Of course, you know that the solutions lie at L[k] = k Pi, but for less simple equations, this would not be known.

Lets continue with a numerical-graphical procedure to find the roots.

Define a function of LL which gives the value of vv[1]

r[LL_] = vv[1] /. dsol

This is a function that can be used in FindRoot.

FindRoot[r[LL],{LL,3}]

r[LL /. %]

But now you are faced with finding a sequence of the values.  How can you choose which values to give to FindRoot for starting values?  One trick is to use a plot to choose points where the graph crosses the axis.

plot = Plot[Chop[r[LL]],{LL,.01,20}]

One way to preceed is to type in values which appear on the graph to bracket the root of interest.  However, the plot contains the information in a useable form if you know how to extract it.

Position[plot,Line[___]]

pvals = plot[[1,1,1,1]];
svals = {};
Do[If[Sign[pvals[[j,2]]] != Sign[pvals[[j+1,2]]],
		svals = Join[svals,{{pvals[[j,1]],pvals[[j+1,1]]}}]],
	{j,1,Length[pvals]-1}];
{n = Length[svals],svals}

Now that we have good starting values, we can apply FindRoot[] to compute the first 6 nonzero eigenvalues.

Table[Chop[LL /. FindRoot[r[LL],{LL,svals[[k]]}]],
	{k,1,n}]

The values are reasonably accurate.

Table[N[(L[k] = k Pi) - %[[k]]],{k,1,n}]

Now that you have found the eigenvalues, it is possible to find the eigenfunctions.

Clear[v];
Table[v[k][x_] = (vv[x] /. dsol /. LL->L[k]),{k,1,n}]

These are eigenfunctions, but for the purposes of finding the expansion to solve the equation, it would be nice to have them as an orthonormal set such that

Integrate[v[k][x] v[j][x],{x,0,1}] == If[k == j,1,0]

Do[v[k][x_] = v[k][x]/Sqrt[Integrate[(v[k][x])^2,{x,0,1}]],
	{k,1,n}]

Now we have an orthonormal basis

TableForm[
	Table[Integrate[v[k][x] v[j][x],{x,0,1}],{k,1,n},{j,1,n}]]

The functions in time are

Table[w[k_][t_] := Cos[L[k] t];w[k][t],{k,1,n}]

Lets find the solution for the equation for the function f[x] = 1-Abs[2(x-.5)].

f[x_] := 1-Abs[2(x-.5)]

Table[a[k] = Integrate[2 x v[k][x],{x,0,1/2}] +
  			 Integrate[(2-2x) v[k][x],{x,1/2,1}],
	{k,1,n}]

The solution to the wave equation (truncated to six modes) is

usol[x_,t_] = Sum[a[k] v[k][x] w[k][t],{k,1,n}]

You can easily produce an animation of the solution to the wave equation.

Do[Plot[Chop[usol[x,t]],{x,0,1},PlotRange->{{0,1},{-1,1}}]
	,{t,0,2,.1}]

A nonuniform string

Now consider a string which has thickened ends.   Assume that the modulus of elsticity and the density are uniform, so that the only dependence is on the cross sectional area.

b = 1/5;
s = 1/4;
c = 25;
p[x_] := 1 + s (Tanh[c(b-x)]+1)/2+ s (Tanh[c(x-(1-b))]+1)/2;
Plot[{p[x],-p[x]},
{x,0,1},PlotRange->All]

Remember, that these are scaled units, not necessarily equal to the physical units.  In the physical units, the string might look more like this.

Plot[{p[x],-p[x]},
{x,0,1},PlotRange->All,AspectRatio->.1]

Before, you were able to use DSolve to find the solution.  With this function, that is no longer possible, but here is where Mathematica's numerical capabilities are very useful.

Clear[vv];
vsol = NDSolve[
	{D[p[x] D[vv[x],x],x] + Pi^2 p[x] vv[x] == 0,
	vv[0] == 0,vv'[0] == 1},
	vv,
	{x,0,1}]

Plot[Evaluate[vv[x] /. vsol],{x,0,1}]

Rather than having to put each value in to the right spot in the NDSolve command, it is better to define a function which gives the solution as a function of LL.

vsol[LL_] := (vv /. NDSolve[
	{D[p[x] D[vv[x],x],x] + LL^2 p[x] vv[x] == 0,
	vv[0] == 0,vv'[0] == 1},
	vv,
	{x,0,1}])[[1]]

vsol[3.3]

Plot[%[x],{x,0,1}]

Now it is easy to define a function which gives the value at the right boundary,

r[LL_] := vsol[LL][1]

r[3.4]

Which can be used in FindRoot

FindRoot[r[LL],{LL,3.3,3.4}]

Here is a comparison of the first eigenfunction for the thick and uniform strings.

Plot[Evaluate[{vsol[LL /. %][x],Sin[Pi x]/Pi}],{x,0,1},
	PlotStyle->{Black,Red}]

You want to get a sequence of the first several eigenvalues.  You could use the trick shown above of making a plot and extracting bracketing starting values, but it could take a long time to produce the plot because each evaluation of the function r[LL] requires solving a differential equation.  Since we know that this problem is a perturbation of the previous, it should suffice to use the previous sequence of starting values.

And the same step (minus the Chop[]) to obtain the eigenvalues: (I put in the "Print[k]") so that I could see where it was.

Do[L[k] = LL /. FindRoot[r[LL],{LL,svals[[k]]}];
	Print["L[",k,] = ",L[k]];
	{k,1,n}]

And the unnormalized eigenfunctions.

Clear[v];
Do[v[k] = vsol[L[k]],{k,1,n}]

Now how to find the orthonormalized eigenfunctions.  We can use the same procedure as before, but will get better accuracy by producing a new InterpolatingFunction object which is the normalized eigenfunction.

First extend the definition of vsol:

vsol[LL_,ux_] := (vv /. NDSolve[
	{D[p[x] D[vv[x],x],x] + LL^2 p[x] vv[x] == 0,
	vv[0] == 0,vv'[0] == ux},
	vv,
	{x,0,1}])[[1]]

Then we can use

Do[Print[k," ",iv = NIntegrate[p[x] (v[k][x])^2,{x,0,1}]];
	v[k] = vsol[L[k],(v[k])'[0]/Sqrt[iv]],
	{k,1,n}]

Here is a plot of the normalized eigenfunctions

Plot[Evaluate[Table[v[k][x],{k,1,n}]],{x,0,1},
	PlotStyle->Table[RGBColor[1-j/(n-1),0,j/n],{j,0,n-1}]]

This shows that we do have an approximately orthonormal set.

TableForm[
	Table[NIntegrate[p[x] v[k][x] v[j][x],{x,0,1}],{k,1,n},{j,1,n}]]

Now we are in a position to find the coefficients for the function f[x] = 1-Abs[2 (x-.5)] in this case.

Do[Print[a[k] = NIntegrate[p[x] f[x] v[k][x],{x,0,.5,1}]],{k,1,n}]

Here is the solution to the wave equation.

thick[x_,t_] = Sum[a[k] v[k][x] w[k][t],{k,1,n,2}];

Do[Plot[thick[x,t],{x,0,1},
PlotRange->{{0,1},{-1,1}}],{t,0,2}]

Automating the process

One of the great powers of Mathematica is taht it alow you to automate a process such as this.  Once you have it automated, you can try new functions or parameters very easily.

Redefine vsol and r to include in its arguments the relevant functions.

vsol[LL_,ux_,p_,rho_,{x_,xmin_,xmax_}] := (vv /. NDSolve[
	{D[p D[vv[x],x],x] + LL^2 rho vv[x] == 0,
	vv[0] == 0,vv'[0] == ux},
	vv,
	{x,xmin,xmax}])[[1]]

r[LL_,ux_,p_,rho_,{x_,xmin_,xmax_}] := 
	vsol[LL,ux,p,rho,{x,xmin,xmax}][1]

Now you can define a function which will do the entire process.   It arguments are the relevant functions, the interval, and a set of starting values for FindRoot[] to use.

Eigenexpansion[f_,p_,rho_,{x_,xmin_,xmax_},svals_] :=
Module[{k,n = Length[svals],iv},
	Do[
		L[k] = LL /. FindRoot[
				r[LL,1/LL,p,rho,{x,xmin,xmax}],
				{LL,svals[[k]]}];
		Print["L[",k,"] = ",L[k]];
		v[k] = vsol[L[k],1/L[k],p,rho,{x,xmin,xmax}];
		iv = NIntegrate[rho (v[k][x])^2,{x,xmin,xmax}];
		v[k] = vsol[L[k],1/(L[k]*Sqrt[iv]),p,rho,{x,xmin,xmax}];
		a[k] = NIntegrate[rho f v[k][x],{x,xmin,xmax}];
		Print["a[",k,"] = ",a[k]],
		{k,1,n}];
	Function[
		Evaluate[{x,t}],
		Evaluate[Sum[a[k] v[k][x] Cos[L[k] t],{k,1,n}]]]];

Generate enough starting values so that the first 11 eigenvalues are found.

svals = Table[N[k Pi + {-.5,.5}],{k,1,11}]

This then calulates the solution for the string with thick ends.

b = .2;
s = .25;
c = 25;
p[x_] := 1 + s (Tanh[c(b-x)]+1)/2+ s (Tanh[c(x-(1-b))]+1)/2;
f[x] = 1-Abs[2(x-.5)];

thick = Eigenexpansion[f[x],p[x],p[x],{x,0,1},svals];

Since the process is so simple, you might want to find out what happens if only one end of the string is thickened.

b = .2;
s = .25;
c = 25;
p[x_] := 1 + s (Tanh[c(b-x)]+1)/2;
f[x] = 1-Abs[2(x-.5)];

half = Eigenexpansion[f[x],p[x],p[x],{x,0,1},svals];

For comparison, find the solution for the uniform string.

usol = Eigenexpansion[f[x],1,1,{x,0,1},svals];

Now lets do some plots with all three together.

Do[Plot[{usol[x,t],thick[x,t],half[x,t]},{x,0,1},
	PlotRange->{{0,1},{-1,1}},
	PlotStyle->{Black,Red,Green}],{t,0,2,.1}]

Do[Plot[{usol[x,t],thick[x,t],half[x,t]},{x,0,1},
	PlotRange->{{0,1},{-1,1}},
	PlotStyle->{Black,Red,Green}],{t,0,20,2}]

Finite Differences

An alternative way of approachine this is to use finite differences:

Let vv[j] represent v[h*j], where h is the constant step size between points x[j] == j*h, and
let pp[i] represent p[h*i] and rrho[i] represent rho[h*i]

Then the operator D[p[x] D[v[x], x], x]/rho[x] can be approximated by

fdiff[j_,h_,vv_,pp_,rrho_]  := 
	(pp[j+1/2]*(vv[j+1]-vv[j])-
	 pp[j-1/2]*(vv[j]-vv[j-1]))/(h^2 rrho[j])

This produces a matrix with coefficients corresponding to the finite difference formula.

fdmat[n_,h_,pp_,rrho_] := Block[{vv,j,k},
Table[
	fd = Expand[fdiff[j,h,vv,pp,rrho]];
	Table[Coefficient[fd,vv[k]],{k,1,n-1}],
	{j,1,n-1}]]

f[x_] := 1-Abs[2(x-.5)]

n = 100;
h = 1./n;
rrho = 1&;
pp = 1&;
fvec = Table[f[h*j],{j,0,n}];

matrix = fdmat[n,h,pp,rrho];

Compute teh eigenvalues and eigenvectors.

{evals,evecs} = Eigensystem[matrix];

Do[Lvec[k] = Sqrt[-evals[[n-k]]];
	Print["L[",k,"] = ",Lvec[k]];
	vvec[k] = Join[{0},evecs[[n-k]],{0}];
	vvec[k] = vvec[k]*Sign[vvec[k][[2]]];
	vvec[k] = vvec[k]/
		Sqrt[h*Sum[rrho[j] (vvec[k][[j+1]])^2,{j,0,n}]];
	a[k] = h*(fvec.vvec[k]);
	Print["a[",k,"] = ",a[k]],
	{k,1,11}]

Compare the approximate eigenfunctions to the exact ones.

Do[Print[Max[Abs[N[vvec[k]-Table[Sqrt[2] Sin[Pi h  k j],{j,0,n}]]]]],
{k,1,11}]

The solution can be expressed in semi-discrete form. For each (continuum) value of t, you have a discrete solution

solvec[t_] = Transpose[
{h*Range[0,n],Sum[a[k]*Cos[Lvec[k] t]*vvec[k],{k,1,11}]}];

This can be plotted using the ListPlot[] command.

Do[ListPlot[solvec[t]],{t,0,2}]

As for the other process, this can be automated.  The arguments are the relevant functions, the number of points to use in the discretizations, the number of eigenfunctions to use, and the interval.

DiscreteEigenexpansion[f_,p_,rho_,n_,nk_,{x_,xmin_,xmax_}] :=
Module[{h,fvec,pp,rrho,evals,evecs},
	h = N[(xmax-xmin)/n];
	fvec = Function[{x},f] /@ (xmin + h*Range[0,n]);
	pp = Function[{j},Block[{x = xmin+h*j},p]];
	rrho = Function[{j},Block[{x = xmin+h*j},rho]];
	{evals,evecs} = Eigensystem[fdmat[n,h,pp,rrho]];
	Do[
		Lvec[k] = Sqrt[-evals[[n-k]]];
		Print["L[",k,"] = ",Lvec[k]];
		vvec[k] = Join[{0},evecs[[n-k]],{0}];
		vvec[k] = vvec[k]*Sign[vvec[k][[2]]];
		vvec[k] = vvec[k]/
			Sqrt[h*Sum[rrho[j] (vvec[k][[j+1]])^2,{j,0,n}]];
		a[k] = h*(fvec.vvec[k]);
		Print["a[",k,"] = ",a[k]],
		{k,1,nk}];
	Function[{t},Evaluate[
		Transpose[{
			xmin+h*Range[0,n],
			Sum[a[k]*Cos[Lvec[k] t]*vvec[k],{k,1,nk}]}]]]]

f[x] = 1-Abs[2(x-.5)];

This can now be used to compute different cases as before.

dusol = DiscreteEigenexpansion[f[x],1,1,100,11,{x,0,1}];

b = .2;
s = .25;
c = 25;
p[x_] := 1 + s (Tanh[c(b-x)]+1)/2+ s (Tanh[c(x-(1-b))]+1)/2;
f[x] = 1-Abs[2(x-.5)];

dthick = DiscreteEigenexpansion[f[x],p[x],p[x],100,11,{x,0,1}];