Orthogonal series and boundary value problems Calculations of some Fourier series

This is an evaluated Mathematica notebook. If you have Mathematica or MathReader (which is available free from WRI, you may download the notebook file. This notebook contains some Mathematica calculations for chapter VI, Notes on a vibrating string, of Harrell's WWW textbook.

A string with variable density and Airy's equation

A realistic string or optical fiber may not be uniform, so some of the simplifying
assumptions in our derivation fo the wave equation are not valid. For example, the
constant c^2 = (spring constant)/(mass density) for the spring may be different at different
positions, leading to an equation of the form

           2            2         2        2             2
   \partial u/\partial t  = (c[x]) \partial  u/\partial x
                                           (ModWE1)

or, as a thorough examination of the derivation of the wave equation reveals, we could
more generally have:

         2            2      
 \partial u/\partial t  

     = p(x) (\partial /\partial x) s(x) \partial u/\partial x 
                                           (ModWE2)

for some potentially complicated positive functions p(x) and s(x). How well does the method of separation of variables do for problems like this? Rather well, actually, although we may have to encounter some new functions.

Model Problem VI.3.. Suppose that the wave speed depends on position, so that
c^2 = 1/(1 + x), 0 < x < 1, DBC at 0 and 1.
Find the normal modes of vibration.

Solution.
When we attempt to solve the equation with the ansatz u[t,x] = T[t] X[x], we find the
following eigenvalue problem.

- X''[x] = (1 + x) mu X[x]
There are actually some special functions, called Airy functions, which solve the ODE
y''[x] = x y. Two independent solutions are called Ai(x) and Bi(x), or, according to
Mathematica , AiryAi[x] and AiryBi[x].

In[1]:=

  Plot[AiryAi[x], {x, -10,10}]

In[2]:=

  Plot[AiryAi[x], {x, -5,3}];
  Plot[AiryBi[x], {x, -5,3}]

Out[2]=

  -Graphics-

The function Bi explodes exponentially to the right, while Ai decays exponentially. They
both oscillate to the left (why?). By changing variables we can get these functions to solve
our eigenvalue equation:

In[3]:=

  D[AiryAi[-mu^(1/3) (x+1)],{x,2}]

Out[3]=

                          1/3
  -(mu (1 + x) AiryAi[-(mu    (1 + x))])

We need a linear combination Ai(- mu^(1/3) (1+x)) + C Bi ( - mu^(1/3) (1+x)) which is 0
at x=0 and 0 at x=1. I avoid the cube root at this stage by letting p = mu^1/3:

In[4]:=

  FindRoot[{AiryAi[- p] + C AiryBi[- p] == 0, \
       AiryAi[- 2 p] + C AiryBi[- 2 p] == 0}, \
          {p,1.5},{C,2}]

Out[4]=

  {p -> 1.87088, C -> 0.819688}

In[5]:=

  Plot[AiryAi[-1.87088 (1+x)] +.819688 AiryBi[-1.87088 (1+x)], \
        {x,0,1}]

Out[5]=

  -Graphics-

The fact that this function has no nodes between 0 and 1, and therefore resembles sin(9 x),
indicates that this is the spatial part of the fundamental (lowest-frequency) mode. The
eigenvalue and normal mode are:

In[6]:=

  mu0 = p /. %%;
  Mode0[t_,x_] = (A0 Cos[mu0 t] + B0 Cos[mu0 t]) *\
          (AiryAi[-1.87088 (1+x)] +.819688 AiryBi[-1.87088 (1+x)])

Out[6]=

  (AiryAi[-1.87088 (1 + x)] + 0.819688 AiryBi[-1.87088 (1 + x)]) 
   
    (A0 Cos[1.87088 t] + B0 Cos[1.87088 t])

In[7]:=

  FindRoot[{AiryAi[- p] + C AiryBi[- p] == 0, \
       AiryAi[- 2 p] + C AiryBi[- 2 p] == 0}, \
          {p,3},{C,-2}]