Unit 5: Boundary Value Problems for Partial Differential Equations

Chapter 24: Wave Equation

Sections24.1: the plucked string

Copyright

Copyright * 2001 by Addison Wesley Longman, Inc.

All rights reserved. No part of this publication may be reproduced, stored in a retrieval system, or transmitted, in any form or by any means, electronic, mechanical, photocopying, recording, or otherwise, without the prior written permission of the publisher. Printed in the United States of America.

Initializations

> restart;

> interface(showassumed=0);
assume(n,integer);

> with(plots):
with(plottools):

Warning, the name changecoords has been redefined

>

The following Maple code defines the PX command which creates a periodic extension of a function. See lesson 2-4-1 for details of installation and behavior.

> PX := proc(h, g)
local f, var, d;
if type(h, procedure) then RETURN(subs(
{'D' = rhs(g) - lhs(g), 'F' = h, 'L' = lhs(g)}, proc(
x::algebraic)
local y;
y := floor((x - L)/D); F(x - y*D)
end))
fi;
if type(g, equation) then var := lhs(g); d := rhs(g) fi;
f := unapply(h, var);
subs({'L' = lhs(d), 'D' = rhs(d) - lhs(d), 'F' = f}, proc(
x::algebraic)
local y;
y := floor((x - L)/D); F(x - y*D)
end)
end:

>

Wave Equation on the Finite String

An idealized string, uniform and without internal friction, is stretched along an -axis from to . The ends are fixed, the string is taut, and it is constrained to move in a vertical plane. The equilibrium position for the string coincides with the interval on the -axis, and displacements in the string are measured vertically by . At time , and at location in , the displacement of the string from equilibrium is . For fixed , the shape of the string is given by .

Figure 24.1, below, shows the shape of the string at a particular instant .

> p1 := plot(x*(1-x),x=0..1,color=black):
p2 := textplot([1,-.05,`L`], font=[TIMES,ITALIC,12]):
p3 := arrow([.5,.09],[.5,0],.01,.03,.2,color=black):
p4 := arrow([.5,.16],[.5,.25],.01,.03,.2,color=black):
p5 := textplot([.5,.125,`u(x,t)`]):
### WARNING: note that `I` is no longer of type `^`
p6 := textplot([.525,.136,`^`]):
display([p||(1..6)],scaling=constrained, labels=[x,u], xtickmarks=[-1,2], ytickmarks=[-1,1], labelfont=[TIMES,ITALIC,12]);

>

Wave Equation (1) for

fixed left end (2)

fixed right end (3)

initial shape (4)

initial velocity (5)

__________________________________________________

Table 24.1

As we will derive in Sections24.4, displacements in the string are governed by , the wave equation (1) in Table 24.1. This is a partial differential equation (PDE) in the two independent variables and . The parameter will turn out to be the wave speed , the speed at which a disturbance can propagate horizontally along the string. The boundary conditions (2) and (3) in Table 24.1 assure that the endpoints of the string remain fized. Because they proscribe values of at the endpoints, they are called Dirichlet conditions. Since the prescribed values are zero, these conditions are homogeneous Dirichlet conditions. The initial conditions (4) and (5) in Table 24.1 give the initial shape and velocity of the string. At first, we will take , so that the string is given a (small) initial displacement , and released. This corresponds to a gentle plucking of the string, much like the action in a harpsichord.

Together, the five equations in Table 24.1 constitute an initial-boundary value problem (or BVP) for the finite string. Our goal is to obtain a solution to this BVP, and to understand how the function so determined describes the physical motion of the string.

>

Example 24.1

We will show that is a solution to the BVP in Table 24.1 when , and . A calculation shows

and

To implement this calculation in Maple, set so in the equation

> q := diff(u(x,t),t,t) = 9*diff(u(x,t),x,x);

>

Write as

> U := sin(x)*cos(3*t);

>

and make the substitution

> simplify(subs(u(x,t)=U,q));

>

Moreover, = 0, as we see by

> eval(subs(x=0,U));

>

and

> eval(subs(x=Pi,U));

>

Finally, = , as we see by

> eval(subs(t=0,U));

>

Figure 24.2 consisting of Figures 24.2a and 24.2 b, shows the shape of the string at two different times.

Figure 24.2a is

> U := sin(x)*cos(3*t):
p1 := plot(subs(t=.3,U),x=0..Pi,color=black, scaling=constrained, xtickmarks=[0], ytickmarks=[-1,1], labels=[x,u], labelfont=[TIMES,ITALIC,12]):
p2 := textplot([Pi,-.1,`p`],font=[SYMBOL,10]):
p3 := textplot([2.7,.56,`t = .3`]):
display([p||(1..3)]);

>

and Figure 24.2b is

> p1 := plot(subs(t=1.5,U),x=0..Pi,color=black, scaling=constrained, xtickmarks=[0], ytickmarks=[-1,1], labels=[x,u], labelfont=[TIMES,ITALIC,12]):
p2 := textplot([Pi,-.1,`p`],font=[SYMBOL,10]):
p3 := textplot([2.7,-.25,`t = 1.5`]):
display([p||(1..3)]);

>

Figure 24.3 shows the solution surface whereby is graphed as a surface in the -plane.

> p1 := plot3d(U,x=0..Pi,t=0..4*Pi/3, axes=frame, labels=[` x`,` t`,`u `], labelfont=[TIMES,ITALIC,12], color=red, style=hidden, tickmarks=[3,4,3], orientation=[45,60]):
p1;

>

Figure 24.4 shows several snapshots of the string superimposed on the solution surface.

> F := z -> spacecurve([x,z,sin(x)*cos(3*z),x=0..Pi], color=black, thickness=3):
p2 := display3d([seq(F(2*Pi/30*k),k=0..20)]):
display3d([p1,p2],orientation=[25,80], scaling=constrained, axes=frame);

>

The solution surface is traced out in time by a succession of moving images of the dynamic string, or equivalently, the plane sections of the solution surface are snapshots in time of the physical motion of the string. The obvious animation of the string tracing out the solution surface follows.

> p2 := display3d([seq(F(2*Pi/30*k),k=0..20)], insequence=true):
display3d([p1,p2],orientation=[25,80], scaling=constrained, axes=frame);

>

Figures 24.3 and 24.4 show that the motion of the string is an example of a standing wave because the points of zero displacement do not translate along the string.

In fact, a complete animation of the motion of the string is given by

> animate(U,x=0..Pi,t=0..2*Pi/3, frames=20, color=black, labels=[x,u], labelfont=[TIMES,ITALIC,12], xtickmarks=3, ytickmarks=3);

>

Example 24.2

A wave free to move along the string is called a traveling wave. Suppose the initial disturbance in the string has the shape shown in figure 24.5, below.

> g := x -> piecewise(x<3*Pi/8,0, x<Pi/2,x-3*Pi/8, x<5*Pi/8,Pi/8, x<3*Pi/4,3*Pi/4-x,0):
plot(g,0..Pi, scaling=constrained, labels=[x,u], labelfont=[TIMES,ITALIC,12], xtickmarks=3, ytickmarks=[0,.4]);

>

Under the action of the wave equation, that initial disturbance would travel along the string, reflecting at the endpoints, and bouncing back and forth between between them. Figure 24.6 (consisting of Figures 24.6a and 24.6b) shows the traveling wave at two successive times after .

Figure 24.6a is

> g1 := piecewise(x<0,-g(-x), x<Pi,g(x)):
G := PX(g1,x=-Pi..Pi):
U1 := (G(x-t)+G(x+t))/2:
pp1:=plot(subs(t=.7,U1),x=0..Pi,color=red,xtickmarks=3, ytickmarks=[-.4,0,.4], scaling=constrained, view=[0..Pi,-.4.. .4]):
pp2 := textplot([2.9,.3,`t = .7`]):
display([pp1,pp2], labels=[x,u], labelfont=[TIMES,ITALIC,12]);

>

while Figure 24.6b is

> pp3 := plot(subs(t=1.47,U1),x=0..Pi, color=red, xtickmarks=3, ytickmarks=[-.4,0,.4], scaling=constrained, view=[0..Pi,-.4.. .4]):
pp4:=textplot([2.9,.3,`t = 1.47`]):
display([pp3,pp4], labels=[x,u], labelfont=[TIMES,ITALIC,12]);

>

The complete animation showing how the initial energy in the string is distributed under the action of the wave equation is given by

> g1 := piecewise(x<0,-g(-x), x<Pi,g(x)):
G := PX(g1,x=-Pi..Pi):
U1 := (G(x-t)+G(x+t))/2:
animate(U1,x=0..Pi,t=0..2*Pi, frames=63, color=red, scaling=constrained, labels=[x,u], labelfont=[TIMES,ITALIC,12], xtickmarks=3, ytickmarks=[-.4,0,.4]);

>

The initial energy splits into two equal parts, with each new wave traveling in opposite directions. When the wave reaches the boundary, it is absorbed, and then re-transmitted, only now, with a reflection. Figure 24.7 shows the solutoin surface for the traveling wave.

> plot3d(U1,x=0..Pi,t=0..2*Pi,axes=boxed, labels=[x,t,u], labelfont=[TIMES,ITALIC,12], color=black, style=hidden, shading=zgrayscale, tickmarks=[[0,1,2,3],5,[0]], grid=[25,50], scaling=constrained, orientation=[50,60]);

>

>

The initial disturbance splits into two waves of half the initial height, one traveling left, and one traveling right. Without changing shape, these waves move in the -plane along characteristics , the lines + = constant, which are therefore said to "carry the initial information."

A more dynamic view of the propagation of the initial shape is given by the following animation.

> ff := z -> spacecurve([x,z,subs(t=z,U1),x=0..Pi], color=black, thickness=3, tickmarks=[2,2,[0]]):
display3d([seq(ff(2*Pi/50*k),k=0..50)], insequence=true, axes=boxed, labels=[x,t,u], labelfont=[TIMES,ITALIC,14], scaling=constrained, orientation=[55,65]);

>

Solution by Separation of Variables

The classical method of solution for such BVPs is the method of separation of variables . Assume a solution of the separated form

a product of one function just containing x and one function just containing t . Recall the initial example

from above.

Now, generalize this to the product , that is, to

> U := X(x)*T(t);

>

to which we next apply the fixed endpoint conditions, obtaining

and

that is,

> subs(x=0,U) = 0;
subs(x=L,U) = 0;

>

Ruling out the choice because that implies = 0 (identically), we conclude = 0.

Under the assumption of a separated solution, the initial conditions become

=

ans

T'(0) =

Unless either or is zero, no conclusion can be drawn from these two equations. Assuming for simplicity that the initial velocity is zero, so that , we conclude

= X(x) T'(0) = 0 => T'(0) = 0

The initial condition can be satisfied with the aid of a Fourier series. Before we can see that, we must apply the separation assumption to the PDE itself, namely

> q := diff(u(x,t),t,t) = c^2*diff(u(x,t),x,x);

>

giving

that is,

> q1 := eval(subs(u(x,t)=U,q));

>

We have just deliberately committed a notational no-no. The derivatives on the left are with respect to , while the derivatives on the right are with respect to . Yet, we have used the same symbol for each differentiation, a mathematical imprecision and ambiguity that is too convenient to avoid.

Next, divide by , with in the form

This general strategy, which seems to work with linear second-order PDEs, gives

that is,

> q2 := q1/(c^2*U);

>

the separated form of the PDE wherein each side of the equation is a function of only one of the variables. Since on the left there is a function purely of t , and on the right, a function purely of x , and these two functions must be equal for all values of x and t , each must be equal to the same constant.

This argument is generally not transparent to the novice. Consider the following. Student A is asked to state a "favorite function of ." Student B is asked to state a "favorite function of ." Typically, the functions now written on the board are something like and . Student C is asked to select a "favorite value of ," and student D is asked for a "favorite value of ." The functions and are now computed at values such as and . The students then realize that because = 0, the only way for a function of to equal, consistently, a function of is for the two functions to be the same constant function.

This constant is typically taken as and is called the Bernoulli separation constant . Hence, we have the two ordinary differential equations

and

which, together with the relevant data, form the two problems

= 0

and

The first problem is a Sturm-Liouville eigenvalue problem, while the second is a modified initial value problem.

Obtaining the two separated ordinary differential equations in Maple is not a simple task. For example, we begin with

> q3 := rhs(q2) = lambda;
q4 := lhs(q2) = lambda;

>

Then, re-arranging, we get

> q5 := lhs(X(x)*q3) - rhs(X(x)*q3) = 0;
q6 := lhs(c^2*T(t)*q4) - rhs(c^2*T(t)*q4) = 0;

>

In effect, we did all the work. Maple really only played the role of bookkeeper, because it does not have any built-in code for performing manipulations of this type.

The solution of this Sturm-Liouville eigenvalue problem, obtained in Sections16.1, is

with eigenvalues

Again, Maple presently does not have the tools with which to construct these solutions, so we enter them as

> Xn := B[n]*sin(n*Pi*x/L);

>

with eigenvalues given by

> lambda[n] := -(n*Pi/L)^2;

>

for , ....

For each value of , the function T( ) is the solution of

so we write

a solution we obtain in Maple by writing the differential equation as

> q7 := subs(lambda=lambda[n],q6);

>

then solving with

> q8 := dsolve({q7,D(T)(0)=0},T(t));

>

Emphasizing the dependence on , we write

> Tn := subs(_C2=C[n],rhs(q8));

>

( Caution : If Maple uses _C1 instead of _C2 for the constant of integration, please adjust the worksheet accordingly.)

A single eigenfunction, is therefore the product

that is,

> Un := Xn*Tn;

>

which contains two constants where only one is needed. Hence, replace the product with the single constant via

> un := algsubs(B[n]*C[n]=b[n],Un);

>

The set of eigensolutions can be thought of as the analog of a fundamental set for a linear ODE. The general solution of a linear ODE is a linear combination of all the members of the fundamental set. Similarly, the most general solution of the linear PDE is a linear combination of all the members of its "fundamental set," namely, the set of eigensolutions .

Since there are an infinite number of eigensolutions, our general solution is the infinite sum

=

that is,

> U := Sum(un,n=1..infinity);

>

The final condition to be applied in this BVP is the initial shape requirement,

=

that is,

> simplify(subs(t=0,U)) = f(x);

>

The coefficients must therefore be the Fourier sine series coefficients for the function . There is no other choice for these constants if this final condition is to be satisfied. If is to be equal to a sum of sine terms, then the coefficients in that sum must be the Fourier sine series coefficients. Therefore, the 's are determined by the Fourier sine series formulas

and the solution is given by the infinite series

This infinite sum, along with the defining integrals for the , constitutes the formal solution to the given BVP for the finite plucked string.

>

Example 24.3

A string of length has its ends fixed at , and on an x -axis. The string is under tension, so that when it is given the initial shape , as shown in Figure 24.8 (below), and released, it oscillates. We solve for the motion of the string at any time .

The associated BVP that models the motion of this string is

, t > 0

=

The solution, given by

and

is implemented in Maple as follows. The initial shape function is

> f := x*(Pi-x)/10;

>

and Figure 24.8 (which shows its graph) is

> plot(f,x=0..Pi,color=black, scaling=constrained, xtickmarks=3, ytickmarks=2, labels=[x,u], labelfont=[TIMES,ITALIC,12]);

>

The Fourier coefficients are given by the integral

> q := (2/Pi)*Int(f*sin(n*x),x=0..Pi);

>

whose value is

> b := simplify(value(q));

>

A peek at the first ten coefficients

> seq(subs(n=k,b),k=1..10);

>

reveals that the odd-indexed terms "survive," the even-indexed terms are zero. Moreover, the coefficients rapidly converge to zero for this function. In fact, they go to zero as , a rapid decrease, as the following conversion to floating point numbers shows.

> seq(evalf(subs(n=k,b)),k=1..10);

>

Hence, a very few terms will be needed to get a good approximation to . We can test this by examining how many terms it takes to get the Fourier series of to approximate well. The first two distinct partial sums in the Fourier series for are

> p1 := sum(b*sin(n*x),n=1..1);
p3 := sum(b*sin(n*x),n=1..3);

>

with (thin black) and (thick red) shown in Figure 24.9, below.

> plot([f,p1],x=0..Pi,color=[black,red], thickness=[1,3], scaling=constrained, xtickmarks=3, ytickmarks=2, labels=[x,u], labelfont=[TIMES,ITALIC,12]);

>

This experiment suggests that two terms of the series for might yield a reasonably good approximation, so write

=

that is,

> U := sum(b*sin(n*x)*cos(c*n*t),n=1..3);

>

as an approximation to the solution. Set , that is,

> U1 := subs(c=1,U);

>

to obtain the graph of the solution surface seen in Figure 24.10, below.

> plot3d(U1, x=0..Pi, t=0..4*Pi, axes=frame, labels=[x,`t `,`u `], labelfont=[TIMES,ITALIC,12], style=hidden, color=red, tickmarks=[3,6,3], orientation=[-45,60]);

>

The plane sections are snapshots in time, freezing the physical motion of the string for each value of . Exhibiting these shapshots in succession forms a movie, or animation, of the physical motion of the string. This animation is given by the following Maple command.

> animate(U1,x=0..Pi,t=0..2*Pi, frames=60, color=black, scaling=constrained, xtickmarks=3, ytickmarks=2, labels=[x,`u `], labelfont=[TIMES,ITALIC,12]);

>

As the string moves in space-time, its history generates a surface, a small part of which we now animate.

>

> F := z -> plot3d(U1, x=0..Pi, t=0..z):
display3d([seq(F(Pi/10*k),k=1..30)],insequence=true, axes=boxed, labels=[x,t,u], labelfont=[TIMES,ITALIC,12], style=hidden, color=red, tickmarks=[3,5,3]);

>