(*^
::[	Information =

	"This is a Mathematica Notebook file.  It contains ASCII text, and can be
	transferred by email, ftp, or other text-file transfer utility.  It should
	be read or edited using a copy of Mathematica or MathReader.  If you 
	received this as email, use your mail application or copy/paste to save 
	everything from the line containing (*^ down to the line containing ^*)
	into a plain text file.  On some systems you may have to give the file a 
	name ending with ".ma" to allow Mathematica to recognize it as a Notebook.
	The line below identifies what version of Mathematica created this file,
	but it can be opened using any other version as well.";

	FrontEndVersion = "Macintosh Mathematica Notebook Front End Version 2.2";

	MacintoshStandardFontEncoding; 
	
	fontset = title, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeTitle, center, M7, bold, L3, e8,  24, "New York"; 
	fontset = subtitle, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeTitle, center, M7, bold, L2, e6,  18, "New York"; 
	fontset = subsubtitle, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeTitle, center, M7, bold, L2, e6,  14, "New York"; 
	fontset = section, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, grayBox, M22, bold, L2, a20,  14, "New York"; 
	fontset = subsection, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, blackBox, M19, bold, L2, a15,  12, "New York"; 
	fontset = subsubsection, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, whiteBox, M18, bold, L2, a12,  10, "New York"; 
	fontset = text, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, L2,  12, "New York"; 
	fontset = smalltext, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, L2,  10, "New York"; 
	fontset = input, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeInput, M42, N23, bold, L2,  12, "Courier"; 
	fontset = output, output, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, L2,  12, "Courier"; 
	fontset = message, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, R65535, L2,  12, "Courier"; 
	fontset = print, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, L2,  12, "Courier"; 
	fontset = info, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, L2,  12, "Courier"; 
	fontset = postscript, PostScript, formatAsPostScript, output, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeGraphics, M7, l34, w282, h287, L2,  12, "Courier"; 
	fontset = name, inactive, nowordwrap, nohscroll, noKeepOnOnePage, preserveAspect, M7, italic, B65535, L2,  10, "Geneva"; 
	fontset = header, inactive, noKeepOnOnePage, preserveAspect, M7, L2,  10, "Times"; 
	fontset = leftheader, inactive, L2,  10, "Times"; 
	fontset = footer, inactive, noKeepOnOnePage, preserveAspect, center, M7, L2,  12, "Times"; 
	fontset = leftfooter, inactive, center, L2,  12, "Times"; 
	fontset = help, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, L2,  10, "Geneva"; 
	fontset = clipboard, inactive, noKeepOnOnePage, preserveAspect, M7, L2,  12, "New York"; 
	fontset = completions, inactive, nowordwrap, noKeepOnOnePage, preserveAspect, M7, L2,  12, "New York"; 
	fontset = special1, inactive, nowordwrap, noKeepOnOnePage, preserveAspect, M7, L2,  12, "New York"; 
	fontset = special2, inactive, nowordwrap, noKeepOnOnePage, preserveAspect, center, M7, L2,  12, "New York"; 
	fontset = special3, inactive, nowordwrap, noKeepOnOnePage, preserveAspect, right, M7, L2,  12, "New York"; 
	fontset = special4, inactive, nowordwrap, noKeepOnOnePage, preserveAspect, M7, L2,  12, "New York"; 
	fontset = special5, inactive, nowordwrap, noKeepOnOnePage, preserveAspect, M7, L2,  12, "New York"; 
	paletteColors = 128; currentKernel; 
]
:[font = title; inactive; preserveAspect; fontSize = 19; fontName = "Times"; startGroup]
Linear Methods of Applied Mathematics
Calculations of some Fourier series
:[font = subtitle; inactive; dontPreserveAspect]
The Method of d'Alembert
:[font = text; inactive; preserveAspect; plain; bold; fontName = "Times"]
(c) Copyright 1994-1997 by Evans M. Harrell II and James V. Herod.     All rights reserved.
:[font = text; inactive; preserveAspect; plain; bold; fontName = "Times"]
Acknowledgment.  We are grateful to Alfred D. Andrew for contributions which have been incorporated into this notebook.
:[font = text; inactive; Cclosed; preserveAspect; plain; bold; fontSize = 13; fontName = "Times"; startGroup]
Notes for the instructor.
:[font = text; inactive; preserveAspect; fontSize = 13; fontName = "Times"; endGroup]
This contains calculations and examples which correlate with chapter 7 of the WWW text by Harrell and Herod.  Students can be encouraged to cut and paste from this notebook to do homework.
:[font = subsubsection; inactive; Cclosed; preserveAspect; fontSize = 11; fontName = "Times"; startGroup]
Instructions
:[font = text; inactive; preserveAspect; fontSize = 13; fontName = "Times"; endGroup; endGroup]
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.
;[s]
11:0,0;19,1;30,0;96,1;133,0;181,1;192,0;315,1;326,0;537,1;548,0;720,-1;
2:6,15,10,Times,0,13,0,0,0;5,15,10,Times,2,13,0,0,0;
:[font = text; inactive; preserveAspect]
As in the preceding chapter, we define the wave operator by:
;[s]
2:0,0;60,1;61,-1;
2:1,17,12,New York,0,12,0,0,0;1,14,9,Times,0,12,0,0,0;
:[font = input; preserveAspect]
WaveOp[f_,c_] := (1/c^2) D[f, {t,2}] -  D[f, {x,2}]
:[font = text; inactive; preserveAspect]
The wave equation is the statement that the unknown, u, is in the null space of this operator.  We'll often set c = 1 for convenience , so let's also define
:[font = input; dontPreserveAspect]
 WaveOp[f_] := WaveOp[f,1]
:[font = section; inactive; preserveAspect; startGroup]
d'Alembert's solution of the wave equation.  

Derivation and discussion.
:[font = text; inactive; preserveAspect; fontName = "Times"]
Let's test out the wave operator on a couple of instructive functions.  We set c=1 fot this discussion.
;[s]
1:0,1;104,-1;
2:0,14,9,Times,0,12,0,0,0;1,17,12,New York,0,12,0,0,0;
:[font = input; preserveAspect; startGroup]
 WaveOp[Sin[t - x]]
:[font = output; output; inactive; dontPreserveAspect; endGroup]
0
;[o]
0
:[font = input; preserveAspect; startGroup]
WaveOp[Sin[Exp[Cos[x-t]]]]
:[font = output; output; inactive; preserveAspect; endGroup]
0
;[o]
0
:[font = text; inactive; preserveAspect]
This can't be a coincidence!  In fact, early in the nineteenth century, d'Alembert realized that if f is any  function depending on the combination x-t, WaveOp[f] is automatically zero:
;[s]
7:0,0;100,2;101,0;105,1;109,0;160,2;161,0;186,-1;
3:4,17,12,New York,0,12,0,0,0;1,17,12,New York,2,12,0,0,0;2,19,13,Symbol,0,12,0,0,0;
:[font = input; dontPreserveAspect]
WaveOp[Phi[x + t]]
:[font = text; inactive; dontPreserveAspect]
The same is true for any function of x+t:
:[font = input; dontPreserveAspect; startGroup]
WaveOp[Psi[x - t]]
:[font = output; output; inactive; dontPreserveAspect; endGroup]
0
;[o]
0
:[font = text; inactive; preserveAspect]
Even better:
:[font = input; dontPreserveAspect; startGroup]
WaveOp[Phi[x + t] + Psi[x - t]]
:[font = output; output; inactive; dontPreserveAspect; endGroup]
0
;[o]
0
:[font = text; inactive; dontPreserveAspect; fontName = "Times"]
Since the wave equation is linear, given any two twice differentiable functions of a single variable, the expression f[x+ct]+ y[x-ct] is also a solution.   D'Alembert had the insight that not only is f[x+ct]+ y[x-ct] a solution, it is the general solution.  To understand this we need to think about how f and y could be determined in specific cases.

Suppose now that we have an infinitely long string, and initial conditions of the usual form:
                    u[0,x] = f[x]
                    ut[0,x] = g[x].
                    
Here there are no boundary conditions to impose, since we let x run to infnity in both directions.


How do we find a solution of the form u[t,x] = Phi[x+ct]+ Psi[x-ct] which solves this problem?  

Suppose first that g[x] := 0.  Then, as in the text, we substitute and find that:
                  f[x] = y[x] = f[x]/2.

;[s]
19:0,1;117,3;118,1;126,3;127,1;200,3;201,1;209,3;210,1;304,3;305,1;310,3;311,1;501,2;502,1;836,3;837,1;843,3;844,1;859,-1;
4:0,14,9,Times,0,12,0,0,0;10,17,12,New York,0,12,0,0,0;1,25,16,New York,64,12,0,0,0;8,19,13,Symbol,0,12,0,0,0;
:[font = subsection; inactive; Cclosed; preserveAspect; startGroup]
Example
:[font = text; inactive; preserveAspect]
Suppose that f[x] is a function with one bump near x=0 such as 
Exp[-x^2], and there is no initial velocity.  e would see the bump break into the superposition of two bumps of half the height, one moving to the right and one to the left.
:[font = input; dontPreserveAspect]
Table[Plot[Exp[-(x+n*.3)^2] /2 + Exp[-(x-n*.3)^2]  /2, 
{x, -4, 4}, PlotRange -> {-.1,1}], {n,0,10}]
:[font = text; inactive; preserveAspect; endGroup]
(You can animate these plots and watch the wave move.) 
;[s]
1:0,1;56,-1;
2:0,17,12,New York,0,12,0,0,0;1,14,9,Times,0,12,0,0,0;
:[font = text; inactive; dontPreserveAspect; fontName = "Times"]
Next, let's turn the tables and suppose that f[x] = 0  but g[x]  0.  As in the text, we learn that: Psi[x] = - Phi[x], and 
       Phi[x] = Integrate[g[x], x] 
Consequently, the solution to (WE) with these IC is
:[font = input; dontPreserveAspect]
u2[g_,t_,x_] := (1/2) Integrate[(g /. x -> z), {z,x - t, x + t}]
:[font = text; inactive; dontPreserveAspect]
If, for example,
;[s]
1:0,0;17,-1;
1:1,14,9,Times,0,12,0,0,0;
:[font = input; dontPreserveAspect]
g[x_] := x/(1 + x + x^2)
:[font = input; preserveAspect; startGroup]
u2[g[x],t,x]
:[font = output; output; inactive; preserveAspect; endGroup]
(ArcTan[(1 - 2*t + 2*x)/3^(1/2)]/3^(1/2) - 
 
    ArcTan[(1 + 2*t + 2*x)/3^(1/2)]/3^(1/2) - 
 
    Log[1 - t + x + (-t + x)^2]/2 + Log[1 + t + x + (t + x)^2]/2)/2
;[o]
        1 - 2 t + 2 x           1 + 2 t + 2 x
 ArcTan[-------------]   ArcTan[-------------]
           Sqrt[3]                 Sqrt[3]
(--------------------- - --------------------- - 
        Sqrt[3]                 Sqrt[3]
 
                            2                           2
    Log[1 - t + x + (-t + x) ]   Log[1 + t + x + (t + x) ]
    -------------------------- + -------------------------) / 2
                2                            2
:[font = input; dontPreserveAspect; endGroup; animationSpeed = 0]
Table[Plot[%, 
{x, -10, 10}, PlotRange -> {-5,5}], {t,1,6}]
:[font = section; inactive; preserveAspect; startGroup]
Automating the general solution
:[font = text; inactive; dontPreserveAspect]
We can try to write the general solution as follows. 
;[s]
2:0,1;53,0;54,-1;
2:1,14,9,Times,0,12,0,0,0;1,17,12,New York,0,12,0,0,0;
:[font = input; dontPreserveAspect]
Clear[u]
Clear[f]
Clear[g]
:[font = text; inactive; dontPreserveAspect]
The general solution (with arbitrary c) is:
;[s]
2:0,1;43,0;44,-1;
2:1,14,9,Times,0,12,0,0,0;1,17,12,New York,0,12,0,0,0;
:[font = input; dontPreserveAspect; endGroup]
u[c_,f_,g_,t_,x_] := \
          (1/2) (f[x + c t] + f[x - c t]) \
          + (1/(2 c)) Integrate[g[x1], {x1,x - c t, x + c t}]
:[font = section; inactive; preserveAspect]
Model Problem VII.2
:[font = text; inactive; dontPreserveAspect]
We use the general solution to solve for the wave equation with the following initial conditions, and to plot the solution:
:[font = input; preserveAspect]
Clear[f,g]
:[font = input; dontPreserveAspect]
f[x_] := 4 Cos[4 x]
g[x_] := x^2 - 1/3
:[font = input; Cclosed; dontPreserveAspect; startGroup]
u[1,f,g,t,x]
:[font = output; output; inactive; preserveAspect; endGroup]
((-t - x)/3 + (-t + x)/3 - (-t + x)^3/3 + (t + x)^3/3)/2 + 
 
  (4*Cos[4*(-t + x)] + 4*Cos[4*(t + x)])/2
;[o]
                          3          3
-t - x   -t + x   (-t + x)    (t + x)
------ + ------ - --------- + --------
  3        3          3          3
-------------------------------------- + 
                  2
 
  4 Cos[4 (-t + x)] + 4 Cos[4 (t + x)]
  ------------------------------------
                   2
:[font = text; inactive; preserveAspect; fontName = "Times"]
Although the formula for the solution is fairly simple, the wave forms it describes can be complicated.  For example, at unit time steps beginning at t=0, this solution looks as shown below.

First we add an optional step which allows Mathematica  to sketch the graphs more efficiently:
;[s]
4:0,1;235,2;246,1;286,0;287,-1;
3:1,14,9,Times,0,12,0,0,0;2,17,12,New York,0,12,0,0,0;1,17,12,New York,2,12,0,0,0;
:[font = input; preserveAspect; startGroup]
u[t_,x_] = %
;[s]
3:0,0;11,1;12,0;13,-1;
2:2,14,10,Courier,1,12,0,0,0;1,14,10,Courier,1,12,65535,0,0;
:[font = output; output; inactive; preserveAspect; endGroup]
((-t - x)/3 + (-t + x)/3 - (-t + x)^3/3 + (t + x)^3/3)/2 + 
 
  (4*Cos[4*(-t + x)] + 4*Cos[4*(t + x)])/2
;[o]
                          3          3
-t - x   -t + x   (-t + x)    (t + x)
------ + ------ - --------- + --------
  3        3          3          3
-------------------------------------- + 
                  2
 
  4 Cos[4 (-t + x)] + 4 Cos[4 (t + x)]
  ------------------------------------
                   2
:[font = input; dontPreserveAspect]
Table[Plot[u[t,x], 
{x, -3, 3}], {t,0,4}]
:[font = text; inactive; dontPreserveAspect]
We next consider how d'Alembert can incorporate boundary conditions

Model Problem VII.3.  As explained in the text, we extend the function f[x_] := x (1-x) from the interval  0 ² x ² 1 to the full range of x by reflecting in the odd manner about both end points:
;[s]
3:0,0;69,1;88,0;264,-1;
2:2,14,9,Times,0,12,0,0,0;1,14,9,Times,1,12,0,0,0;
:[font = input; dontPreserveAspect]
Clear[f]
:[font = input; dontPreserveAspect]
f[x_] := x (1 - x) /; 0 <= x <= 1
f[x_] := x (1 + x) /; -1 <= x < 0
:[font = input; dontPreserveAspect]
Plot[f[x], {x,-1,1}]
:[font = input; dontPreserveAspect]
f[x_] := - (2 - x) (1 - (2 - x)) /; 1 < x <= 2
:[font = input; dontPreserveAspect; startGroup]
Plot[f[x], {x,-1,2}]
:[font = text; inactive; dontPreserveAspect]
This new territory can then be flipped over to the interval [-2,-1], etc:
;[s]
1:0,0;74,-1;
1:1,14,9,Times,0,12,0,0,0;
:[font = input; dontPreserveAspect; endGroup]
f[x_] := x (1 - x) /; 0 <= x <= 1
f[x_] := x (1 + x) /; -1 <= x < 0
f[x_] := - (2 - x) (1 - (2 - x)) /; 1 < x <= 2
f[x_] := + (2 + x) (1 - (2 + x)) /; -2 <= x < -1
f[x_] := + (x - 2) (1 - (x - 2)) /; 2 < x <= 3
:[font = input; dontPreserveAspect; startGroup]
Plot[f[x], {x,-2,3}]
:[font = output; output; inactive; preserveAspect]
Graphics["<<>>"]
;[o]
-Graphics-
:[font = text; inactive; dontPreserveAspect]
We could go on like this all day, or we could use modular arithmetic (the int and frac commands in programming languages, or, in Mathematica, Floor[x] and x - Floor[x] ) to automate it better, as follows.  Int and Frac give the integer and fractional parts of a number:
;[s]
3:0,0;129,1;140,0;270,-1;
2:2,14,9,Times,0,12,0,0,0;1,14,9,Times,2,12,0,0,0;
:[font = input; dontPreserveAspect]
Int[x_] := Floor[x]
Frac[x_] := x - Floor[x]
:[font = input; dontPreserveAspect; startGroup]
{Int[3.89], Frac[3.89], Int[-3.89], Frac[-3.89]}
:[font = output; output; inactive; dontPreserveAspect; endGroup; endGroup]
{3, 0.8900000000000000002, -4, 0.1099999999999999999}
;[o]
{3, 0.89, -4, 0.11}
:[font = text; inactive; dontPreserveAspect]
Notice that when we piece the function together, it is periodic with period 2.  Here is an alternative way to write the extended function f(x):
;[s]
1:0,0;144,-1;
1:1,14,9,Times,0,12,0,0,0;
:[font = input; dontPreserveAspect]
f[x_] := If[0 <= Frac[x/2] < 1/2, Frac[x] (1 - Frac[x]), \
              Frac[x] (Frac[x] - 1)]
           
:[font = input; dontPreserveAspect; startGroup]
Plot[f[x], {x,-2,3}]
:[font = output; output; inactive; dontPreserveAspect; endGroup]
Graphics["<<>>"]
;[o]
-Graphics-
:[font = text; inactive; dontPreserveAspect]
Finally, we solve the IBVP:
;[s]
1:0,0;28,-1;
1:1,14,9,Times,0,12,0,0,0;
:[font = input; dontPreserveAspect]
g[x_] := 0
:[font = input; dontPreserveAspect]
Table[Plot[u[1,f,g,t,x], 
{x, -1, 2}, PlotRange -> {-.3,.3}],{t,0,1.6,.2}]
:[font = text; inactive; dontPreserveAspect]
It might be interesting to animate these pictures.  Notice that the wave does not simply jiggle up and down; its shape is changing, too.
;[s]
1:0,0;137,-1;
1:1,14,9,Times,0,12,0,0,0;
^*)