Linear Methods of Applied Mathematics
Calculations of some Fourier series

Solving the Wave Equation

(c) Copyright 1994-1997 by Evans M. Harrell II and James V. Herod.     All rights reserved.

Notes for the instructor.

This contains calculations and examples which correlate with chapter 6 of the WWW text by Harrell and Herod.  Students can be encouraged to cut and paste from this notebook to do homework.

Instructions

This notebook uses Mathematica  to perform calculations for Harrell and Herod's hypertext book, Linear Methods of Applied Mathematics.   The student needs only a basic knowledge of Mathematica  to use the notebook, which is designed both to show how to work problems in the text and to provide a template for using Mathematica  to work other problems of the student's own design.   Calculations will be performed when the reader presses enter in a given calculation cell (bold print).  It is best to activate the cells in order, so that Mathematica  will be able to call on operators and functions defined in earlier cells.  Red color coding is used to warn the reader when a given calculation relies on an earlier one.

A useful substitution, which we shall often make, is:

TrigId = {Cos[Pi n_] -> (-1)^n, Sin[Pi n_] -> 0};

In Mathematica  notation, the wave equation could be written 
					
					D[u, {t,2}] = c^2 D[u, {x,2}] 
					
We classify this PDE as second order, linear, and homogeneous.

Definitions

"Second order" means that derivatives up to second order occur in the equation, but no higher-order derivatives.  Linear means that the differential equation is formed from a linear combination of the solution u and its derivatives, and that u comes into the equation only in this linear combination.  Homogeneous  

This will be the case for the great majority of PDE's in this class.

Discussion of the wave operator and the method of separation of variables

We shall solve the wave equation in two entirely different ways.  The method in chapter VI is known as separation of variables, and later we shall solve the same equation with d'Alembert's method of characteristics.  The first method is often particularly good for linear equations.  "Linear" refers to the superposition property, i.e., if two functions u1 and u2 are solutions, then so is any linear combination, a1 u1 + a2 u2. 

In fact, (WE) is a statement that the function u is in the kernel of a linear operator. 

A convenient way to set up a linear differential operator in Mathematica is as follows - the operator acts on functions, so identify a function f with an understroke to show that it is the independent variable in the 
wave operator:

WaveOp[f_,c_] := (1/c^2) D[f, {t,2}] -  D[f, {x,2}]

We'll often set c = 1 for convenience , so let's also define

WaveOp[f_] := WaveOp[f,1]

(Equivalently, we could choose a clock with time variable t' = c t.)  Let's
test out the wave operator

WaveOp[t Sin[x]^2]

WaveOp[Sin[t - x]]

We next try the ansatz that u[t,x] is of the special form T[t] X[x], known as a product solution.

u[t_,x_] := T[t] X[x]

WaveOp[u[t,x]]

This simplifies if we divide through by T[t] X[x]:

%/(T[t] X[x])

Simplify[%]

If u solves the wave equation, this expression must be zero. (We could quibble about what happens when T or X = 0, but that would either mean that u is the uninteresting "trivial" solution or else the problem occurs only at isolated values of x and t.)  Thus 
                        T''/T = X''/X.
Since the left side is independent of x and the right side is independent of t, both sides must in fact equal a constant. We don't know its value yet, so let us call it - mu. (It turns out to be negative.)
Thus X''[x] = - mu X[x]:

DSolve[- X''[x]== mu X[x],X[x],x]

We now remind Mathematica about Euler's interesting formula for complex exponentials:

% /. {E^(I Sqrt[mu] x) -> Cos[Sqrt[mu] x] + 
               I Sin[Sqrt[mu] x],
         E^(- I Sqrt[mu] x) -> Cos[Sqrt[mu] x] -
                  I Sin[Sqrt[mu] x]}

Evidently, we could reexpress this as

Clear[X,D1,D2]

X[x_] := D1 Cos[Sqrt[mu] x] + 
        D1 Sin[Sqrt[mu] x]

X[0]

This means that in order to satisfy the BC at x=0, we must have:

X[x_] := Sin[Sqrt[mu] x]

Solve[X[L] == 0, mu]

Clear[X]

X[x_,n_Integer] := Sin[n Pi x/L]

Since the second derivative of this is - (n Pi/L)^2 Sin[n Pi x/L], we see that

mu[n_Integer] := -(n Pi/L)^2

Now let's look at the other part of the PDE which we "separated off":
T''[t] = - mu T[t]

As we saw, there are many possible values of \mu which are consistent with the boundary condition for the function X, and the same possibilities occur in equation (6.3) (multiply \mu by c^2 here if c is not set to 1). Therefore
         T[t_,n_] = A[n] Cos[n \pi t] + B[n] Sin[n \pi t ]. 
         
We cannot go further with the function T without being given the initial condition, that is, information about what happens at t=0.

The general solution we come up with is a linear combination of the particular solutions we get by separating the equation. It is not yet obvious that this is a completely general solution, but it is. The most general linear combination is of the form:

u[t_,x_] = Sum[(A[n] Cos[n Pi t ] + B[n] Sin[n Pi t ] )           Sin[n Pi x],  {n,1,infinity}]

The coefficients A and B have to be determined from the initial conditions.

Suppose that 

										u[0,x] = f[x]

and 

                    ut[0,x] =   g[x].
                    
Then we get:

A[n_,f_] := 2 Integrate[f[x] Sin[n Pi x], {x,0,1}]
B[n_,g_] := (2/(n Pi)) Integrate[g[x] Sin[n Pi x], {x,0,1}]

Illustration of normal mode solutions.  (Thanks to Alfred D. Andrew for this section.)

The normal modes are the solutions which have a single frequency (pure frequency) in time.  For the one dimensional wave equation, for each frequency there is basically a single product solution  T[n,t] X[n,t] which is a normal mode for that frequency.  All others are simply phase-shifted versions of this one.  You can sketch a few of them by executing the following cells, and the results can be animated to show how the string vibrates.

T[n_,t_] = Cos[ n Pi t]

X[n_,x_] = Sin[n Pi x]

U[n_,x_,t_] = T[n,t] X[n,x]

Table[Plot[U[1,x,t],{x,0,1},
	PlotRange->{-1,1},
	PlotStyle->
		Thickness[.007]],
	{t,0,.9,.1}]

Table[Plot[U[3,x,t],{x,0,1},
	PlotRange->{-1,1},
	PlotStyle->
		Thickness[.007]],
	{t,0,.3,.1/3}]

Model Problem.  
Solve the problem of the vibrating string for 0 < x < 1, c=1, 
         u[t,0] = u[t,1] = 0
using Fourier sine series, and initial conditions 
              u(x,0) = x - x^3, 
              u_t(x,0) = 0.

TrigId = {Cos[Pi n_] -> (-1)^n, Sin[Pi n_] -> 0}

f[x_] := x - x^3

A[n,f]

(2*(-6*n*Pi*Cos[n*Pi] + 6*Sin[n*Pi] - 2*n^2*Pi^2*Sin[n*Pi]))/
 
  (n^4*Pi^4)

% /. TrigId

(-12*(-1)^n)/(n^3*Pi^3)

The B's are clearly zero, even without the aid of software.  Thus the solution for times t>0 is:

U[t_,x_,infinity_] :=        - Sum[Cos[n Pi t] Sin[n Pi x] 12 (-1)^n /(n^3 Pi^3),            {n,1,infinity}]

Table[Plot[U[t/4,x,7], {x, 0, 1}, 
        PlotRange -> {-.5,.5}], {t,0,8}]

If you are reading this notebook with <I>Mathematica</I> or MathReader you may
select the <BR>
graphics just calculated and animate them to see the string vibrate.

Model Problem VI.1.  
Solve the problem of the vibrating string for 0 < x < 1, c=1, 
         u[t,0] = u[t,1] = 0
using Fourier sine series, and initial conditions 
              u(x,0) = 0, 
              u_t(x,0) = 1

TrigId = {Cos[Pi n_] -> (-1)^n, Sin[Pi n_] -> 0}

The A coefficients are clearly zero this time, even without the aid of software.  As for the B coefficients,

g[x_] := 1

B[n,g]

% /. TrigId

U[t_,x_,infinity_] :=         Sum[Sin[n Pi t] Sin[n Pi x] 4 /(n^2 Pi^2),            {n,1,infinity, 2}]

The final 2 here ensures that only odd terms are included.  The even coefficients are 0

Table[Plot[U[t/4,x,7], {x, 0, 1}, 
        PlotRange -> {-.5,.5}], {t,0,8}]

If you are reading this notebook with Mathematica or MathReader you may select the graphics just calculated and animate them to see the string vibrate.

Model Problem VI.2.  
Set up the general solution of the vibrating string with Neumann BC., with the simplifications c = L = 1.

X[x_,mu_] := D1 Cos[Sqrt[mu] x] + D1 Sin[Sqrt[mu] x]

D[X[x,mu],x]

Possibility 1. mu=0, which means that the output just given is not really the solution to (6.1), which has become X''[x] = 0.  The solutions of this are of the form A + Bx, but B has to be 0 so that X'[0] = 0. Thus with this possibility, X[x] = 1 (times any constant). Eq. (6.3) now also becomes 
          T''(t) = 0, 
          which leads to the product solution of (WE):
	(A0 + B0 t) 1 = A0 + B0 t.

Since the boundary conditions do not apply to the t variable, we keep the
general linear combination.  A0 and B0 will be determined by the initial conditions.

Possibility 2. mu > 0. The boundary condition at x=0 now forces 
    X[x] =  Cos[Sqrt[mu] x] 
(again, up to a constant multiple), and the other boundary condition forces one of the 
values mun = (n Pi)^2.  The product solutions become
	(An Cos[n Pi t] + Bn Sin[n Pi t] ) Cos[n Pi t],
and the general solution will be

u[t_,x_] = A[0] + B[0] t + Sum[(A[n] Cos[n Pi t ] + B[n] Sin[n Pi t ] )                  Cos[n Pi x], {n,1,infinity}]