Partial Differential EquationsThe vibrating string: d'Alembert's formula.Anton DzhamayDepartment of MathematicsThe University of MichiganAnn Arbor, MI 48109wPage: http://www.math.lsa.umich.edu/~adzhamemail: adzham@umich.eduCopyright \251 2004 by Anton DzhamayAll rights reserved

Some packages that we use in this worksheet:restart: with(plottools): with(plots):

The goal of this worksheet is to illustrate the use of the method of characteristics and d'Alembert's formula for the one-dimentional wave equation (i.e., the vibrating string). The equation is NiMvLSUlZGlmZkc2JC0lInVHNiQlInhHJSJ0Ry0lIiRHNiRGKyIiIyomJSJjR0YvLUYlNiRGJy1GLTYkRipGLyIiIg== . Since the equation is second order in NiMlInRH , we need to specify the initial position NiMvLSUidUc2JCUieEciIiEtJSJmRzYjRic= and the initial velocity NiMvLSUlZGlmZkc2JC0lInVHNiQlInhHIiIhJSJ0Ry0lImdHNiNGKg== . If we consider this equation on the whole real line, then we do not need to worry about boundary conditions and the solution is given by the d'Alembert's formula:u:=(x,t)->(f(x-c*t)+f(x+c*t))/2+int(1/(2*c)*g(v),v=x-c*t..x+c*t);For semi-infinite and finite intervals with homogeneous Dirichlet or Neumann boundary conditions the strategy is to extend the initial data to the whole real line in such a way that the boundary conditions are satisfied. It turns out that one has to take the odd extension about the boundary point for the Dirichlet boundary condition and the even extension for the Neumann boundary condition. In the examples below we illustrate this technique.

First we initialize some of the variables. Below NiMlIlRH is the time interval we are interested in, NiMlJG5jaEc= is the number of characteristic curves that we want to plot per unit interval, and NiMlKG50c3RlcHNH is the number of steps in the NiMlInRH variable.c:=2: nch:=3: T:=3: ntsteps:=25:The shift command below allows us to combine the graph of the characteristics (which is two-dimensional) with the graph of the solution (three-dimensional).shift:=transform((x,y)->[x,y,-0.1]):In order to satisfy different boundary conditions we need to take even and odd extensions of functions about NiMvJSJ4RyIiIQ== for half-line or about NiMvJSJ4RyIiIQ== and NiMvJSJ4RyUiTEc= for a finite interval.We define the oddext0 (evenext0) procedures that takes a function or an expression expr in a variable var and returns a functional operator in var that is the odd (even) extension of expr about the origin. oddext0:=proc(expr,var) unapply(signum(var)*unapply(expr,var)(abs(var)),var); end proc:evenext0:=proc(expr,var) unapply(unapply(expr,var)(abs(var)),var); end proc:The procedures for even, odd, and periodc extensions of fucntions over the finite interval are defined similarly.pext:=proc(expr,var,L) unapply(unapply(expr,var)(var - floor((var+L)/(2*L))*2*L),var) end proc: oddext:=proc(expr,var,L) local x; x:=var - floor((var+L)/(2*L))*2*L; unapply(signum(x)*unapply(expr,var)(abs(x)),var) end proc: evenext:=proc(expr,var,L) local x; x:=var - floor((var+L)/(2*L))*2*L; unapply(unapply(expr,var)(abs(x)),var) end proc:

In this example we take NiMvLSUiZ0c2IyUieEciIiE= and take NiMtJSJmRzYjJSJ4Rw== to look like a "guitar pick" (note that this function has discontinuous derivatives, and so the corresponding solution is the weak solution):f:=x->piecewise(x < 1,0,x < 2,x-1,x < 3,3-x,0):g:=x->0: plot([f(x),g(x)],x=-1..5,color=[red,green],thickness=2, scaling=constrained,title="Initial Conditions",legend=["f(x)","g(x)"]);Now we can see how the initial profile changes in time:display([seq(plot(u(x,t/ntsteps),x=-1..5,color=red,thickness=2), t=0..T*ntsteps)],scaling=constrained,insequence=true,title="Pick");The plot below helps to see how the initial profile splits into the left and right moving waves.display([seq(plot([u(x,t/ntsteps),1/2*f(x-c*t/ntsteps),1/2*f(x+c*t/ntsteps)],x=-1..5, color=[red,green,blue],thickness=2), t=0..T*ntsteps)],scaling=constrained,insequence=true);And this is how the graph of the solution looks like (from the top). We color it according to the height of the solution, and this makes it easy to see the left and right waves. The characteristics are plotted below, rotate the graph to see them.char_r:=display([seq(plot([c*t+1/nch*x,t,t=0..T],color=red,thickness=1),x=-1*nch..5*nch)]): char_l:=display([seq(plot([-c*t+1/nch*x,t,t=0..T],color=black,thickness=1),x=-1*nch..5*nch)]):display([plot3d(u(x,t),x=-1..5,t=0..T,axes=boxed,scaling=constrained,numpoints=2000, style=patchnogrid,shading=zhue), shift(char_r),shift(char_l)],orientation=[-90,0],view=[-1..5,0..T,-0.1..1]);

In this example we take NiMvLSUiZkc2IyUieEciIiE= and take NiMtJSJnRzYjJSJ4Rw== to represent kicking the string by a hammer:f:=x->0:g:=x->piecewise(x < 2,0,x < 3,1,0): plot([f(x),g(x)],x=-1..5,color=[red,green],thickness=2, scaling=constrained,title="Initial Conditions",legend=["f(x)","g(x)"]);Now we can see how the initial profile changes in time (this calculation may take a while):display([seq(plot(u(x,t/ntsteps),x=-1..5,color=red,thickness=2), t=0..T*ntsteps)],scaling=constrained,insequence=true,title="Hammer Blow");View from the top and the characteristics:char_r:=display([seq(plot([c*t+1/nch*x,t,t=0..T],color=red,thickness=1),x=-1*nch..5*nch)]): char_l:=display([seq(plot([-c*t+1/nch*x,t,t=0..T],color=black,thickness=1),x=-1*nch..5*nch)]):display([plot3d(u(x,t),x=-1..5,t=0..T,axes=boxed,scaling=constrained,numpoints=2000, style=patchnogrid,shading=zhue), shift(char_r),shift(char_l)],orientation=[-90,0],view=[-1..5,0..T,-0.1..0.3]);Note that there are different (color) zones, in each of these zones the solution is given by a different analytic formula.The following "close-up" picture helps to see the zones more clearlydisplay([plot3d(u(x,t),x=1.5..3.5,t=0..1,axes=boxed,numpoints=2000, style=patchcontour,shading=zhue), shift(char_r),shift(char_l)],orientation=[-90,0],view=[1.5..3.5,0..1,-0.1..0.3]);

In order to satisfy the Dirichlet BC at NiMvJSJ4RyIiIQ== we need to extend our functions from the right half-line to the whole line by taking the odd extension about 0. ff:=piecewise(x < 1,0,x < 2,x-1,x < 3,3-x,0):g:=x->0: f:=oddext0(ff,x): plot([ff(x),f(x)],x=-5..5,color=[red,green],thickness=2, scaling=constrained,title="f(x) and its odd extension");f(x);Now we can see how the initial profile changes in time:T:=4:display([seq(plot(u(x,t/ntsteps),x=0..5,color=red,thickness=2), t=0..T*ntsteps)],scaling=constrained,insequence=true,title="Left Wall (Dirichlet)");It looks like the left-moving wave hits the wall and gets flipped and reflected. Note that NiMvLSUidUc2JCIiISUidEdGJw== for all NiMlInRH . Peeking "behind" the wall we see the "image" right-moving wave emerging from the wall simultaneously with the left-moving wave disappearing into the wall:movie:=display([seq(plot([u(x,t/ntsteps),1/2*f(x-c*t/ntsteps),1/2*f(x+c*t/ntsteps)],x=-5..5, color=[red,green,blue],thickness=2), t=0..T*ntsteps)],scaling=constrained,insequence=true):display([movie,rectangle([-5,5],[0,-5],color=grey)],view=[-3..5,-1..1],title="Left Wall (Dirichlet)");And this is how the graph of the solution looks like (from the top). plot3d(u(x,t),x=0..5,t=0..T,axes=boxed,scaling=constrained,numpoints=2000, style=patchnogrid,shading=zhue,orientation=[-90,0]);

In order to satisfy the Neumann BC at NiMvJSJ4RyIiIQ== we need to extend our functions from the right half-line to the whole line by taking the even extension about NiMvJSJ4RyIiIQ== .ff:=x->piecewise(x < 1,0,x < 2,x-1,x < 3,3-x,0):g:=x->0: f:=evenext0(ff(x),x): plot([ff(x),f(x)],x=-5..5,color=[red,green],thickness=2,numpoints=2000, scaling=constrained,title="f(x) and its even extension");T:=4:display([seq(plot(u(x,t/ntsteps),x=0..5,color=red,thickness=2), t=0..T*ntsteps)],scaling=constrained,insequence=true,title="Left Wall (Neumann)");Again, let's "peek" behind the wall:movie:=display([seq(plot([u(x,t/ntsteps),1/2*f(x-c*t/ntsteps),1/2*f(x+c*t/ntsteps)],x=-5..5, color=[red,green,blue],thickness=2), t=0..T*ntsteps)],scaling=constrained,insequence=true):display([movie,rectangle([-5,5],[0,-5],color=grey)],view=[-3..5,-1..1],title="Left Wall (Neumann BC)");plot3d(u(x,t),x=0..5,t=0..T,axes=boxed,scaling=constrained,numpoints=2000, style=patchnogrid,shading=zhue,orientation=[-90,0]);

In order to satisfy the Neumann BC at NiMvJSJ4RyUiTEc= we take an even extension about NiMvJSJ4RyUiTEc= and then to satisfy the Dirichlet BC at NiMvJSJ4RyIiIQ== we take an odd extension about the origin. This gives us the correct extension on the interval NiM3JCwkKiYiIiMiIiIlIkxHRichIiJGJQ== , and then we extend it periodically to the whole real line. (In the example below NiMvJSJMRyIiJQ== .) This is a very non-optimal way to define such an extension, but it is enough for our purposes.ff:=x->piecewise(x < 1,0,x < 2,x-1,x < 3,3-x,0):g:=x->0: f:=pext(oddext0(evenext(ff(x),x,4)(x),x)(x),x,8): plot([ff(x),f(x)],x=-10..10,color=[red,green],thickness=2, scaling=constrained,title="f(x) and its extension");Now we can see how the initial profile changes in time:T:=8: ntsteps:=15:display([seq(plot(u(x,t/ntsteps),x=0..4,color=red,thickness=2), t=0..T*ntsteps)],scaling=constrained,insequence=true,title="Left Dirichlet and Right Neumann BC");"Peeking behind the walls" shows us the left and right moving waves:movie:=display([seq(plot([u(x,t/ntsteps),1/2*f(x-c*t/ntsteps),1/2*f(x+c*t/ntsteps)],x=-5..9, color=[red,green,blue],thickness=2), t=0..T*ntsteps)],scaling=constrained,insequence=true):display([movie,rectangle([-5,5],[0,-5],color=grey), rectangle([4,5],[9,-5],color=grey)],view=[-5..9,-1..1]);Looking at the graph of the solution (from the top) we can see the characteristics "bounching off the walls":plot3d(u(x,t),x=0..4,t=0..T,axes=boxed,scaling=constrained,numpoints=2000, style=patchnogrid,shading=zhue,orientation=[-90,0]);

This time we need to take the odd extension of NiMtJSJmRzYjJSJ4Rw== about both NiMvJSJ4RyIiIQ== and NiMvJSJ4RyUiTEc= . ff:=x->x^2*(1-x):g:=x->0: f:=oddext(ff(x),x,1): plot([ff(x),f(x)],x=-3..3,color=[red,green],thickness=2, scaling=constrained,title="f(x) and its odd extension",view=-1/2..1/2);T:=4:display([seq(plot(u(x,t/ntsteps),x=0..1,color=red,thickness=2), t=0..T*ntsteps)],scaling=unconstrained,insequence=true);And here are the left and right-waves (set animation to cycle for a better movie):movie:=display([seq(plot([u(x,t/ntsteps),1/2*f(x-c*t/ntsteps),1/2*f(x+c*t/ntsteps)],x=-2..3, color=[red,green,blue],thickness=2), t=0..T*ntsteps)],scaling=constrained,insequence=true):display([movie,rectangle([-5,0.2],[0,-0.2],color=grey),rectangle([1,0.2],[3,-0.2],color=grey)], scaling=unconstrained);plot3d(u(x,t),x=0..1,t=0..3,axes=boxed,scaling=constrained,numpoints=2000, style=patchnogrid,shading=zhue,orientation=[-90,0]);

Walter A. Strauss, Partial Differential Equations: An Introduction, Wiley, 1992Richard Haberman, Elementary Applied Partial Differential Equations, 3rd edition, Prentice Hall

"While every effort has been made to validate the solutions in this worksheet, Waterloo Maple Inc. and the contributors are not responsible for any errors contained and are not liable for any damages resulting from the use of this material."