• Packages
• Introduction
• Initialization
• Example 1: "Pick"
• Example 2: "Hammer Blow"
• Example 3: Left Wall (Dirichlet BC)
• Example 4: Left Wall (Neumann BC)
• Example 5: Two Walls (Mixed)
• Example 6: Two Walls (Dirichlet)
• References
• Disclaimer

Partial Differential Equations

The vibrating string: d'Alembert's formula.

Anton Dzhamay

Department of Mathematics

The University of Michigan

Ann Arbor, MI 48109

wPage: http://www.math.lsa.umich.edu/~adzham

email: adzham@umich.edu

Copyright  2004  by Anton Dzhamay

All rights reserved

Packages

Some packages that we use in this worksheet:

>	restart: with(plottools): with(plots):

Warning, the names arrow and changecoords have been redefined

Introduction

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   . Since the equation is second order in   , we need to specify the initial position    and the initial velocity   . 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.

>

Initialization

First we initialize some of the variables. Below    is the time interval we are interested in,    is the number of characteristic curves that we want to plot per unit interval, and    is the number of steps in the    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    for half-line or about      and    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:

Example 1: "Pick"

In this example we take    and take    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]);

Example 2: "Hammer Blow"

In this example we take    and take    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 clearly

>	display([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]);

Example 3: Left Wall (Dirichlet BC)

In order to satisfy the Dirichlet BC at    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    for all   . 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]);

Example 4: Left Wall (Neumann BC)

In order to satisfy the Neumann BC at    we need to extend our functions from the right half-line to the whole line by taking the even extension about   .

>	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]);

Example 5: Two Walls (Mixed)

In order to satisfy the Neumann BC at    we take an even  extension about    and then to satisfy the Dirichlet BC at    we take an odd  extension about the origin. This gives us the correct extension on the interval   , and then we extend it periodically to the whole real line. (In the example below   .) 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]);

Example 6: Two Walls (Dirichlet)

This time we need to take the odd  extension of    about both    and   .  

>	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]);

>

References

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

Disclaimer

"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."