The wave equation

This worksheet shows how to solve the wave equation governing the motion of an elastic string of length L (assumed to be elastic) with given initial displacement and initial velocity.

sm311o_wave1.mws,wdj,2-9-98

Series method, fixed ends

First, we load the program which computes the sine series coefficients.

> sine_series_coeff:=proc(f,L,n::integer)
#f is a function of x, 0<x<L
#n is the index of b_n
#L is half the period
local x,b;
b:=2*int(f(x)*sin(n*Pi*x/L),x=0..L)/L;
RETURN(b);
end;

The function below will represent the initial displacement of the string of length L=Pi.

> f:=x->piecewise(x<Pi/2,-1,x>Pi/2,2);

> plot(f(x),x=0..Pi);

Below we plot f(x) against several partial Sine series to illustrate how they approximate the function.

> S2:=add(sine_series_coeff(f,Pi,n)*sin(n*x),n=1..2);
S5:=add(sine_series_coeff(f,Pi,n)*sin(n*x),n=1..5):
S10:=add(sine_series_coeff(f,Pi,n)*sin(n*x),n=1..10):
S15:=add(sine_series_coeff(f,Pi,n)*sin(n*x),n=1..15):
plot([f(x),S2,S5,S10,S15],x=0..Pi);

Let u(x,t) denote the displacement from equilibrium (with downward positive) of the point x on the string at time t. Let f(x) denote the initial displacement, u(x,0)=f(x), let g(x) denote the initial velocity, . The program below solves the wave equation

> wave:=proc(f,g,L,a,m::integer,x,t)
local n,s;
s:=add((sine_series_coeff(f,L,n)*cos(a*n*Pi*t/L)+
(L/(a*n*Pi))*sine_series_coeff(g,L,n)*sin(a*n*Pi*t/L))*sin(n*x*Pi/L),n=1..m):
RETURN(s);
end;

Our initial velocity:

> g:=x->1;

Below we plot g(x) against several partial Sine series to illustrate how they approximate the function.

> S2:=add(sine_series_coeff(g,Pi,n)*sin(n*x),n=1..2);
S10:=add(sine_series_coeff(g,Pi,n)*sin(n*x),n=1..10):
plot([g(x),S2,S10],x=0..Pi);

The first few tems of the series solution of the wave equation. Unlike the heat equation, this series solution may not converge quickly so the approximation a few terms gives to the true solution may be poor.

> wave(f,g,Pi,1,3,x,t);

> with(plots):

> animate(wave(f,g,Pi,1,50,x,t),x=0..Pi,t=0..(2*Pi),frames=50,title=`wave equation`);

D'Alembert's method

> restart;

> u:=(x,t)->U(x+a*t)+V(x-a*t);

> a^2*Diff(u(x,t),x$2)-Diff(u(x,t),t$2)=0;
expand(a^2*diff(u(x,t),x$2)-diff(u(x,t),t$2));
MAPLE can check that this satisfies the wave equation:

Let's solve the wave equation with a=1, initial displacement f(x)=x^2 and initial velocity g(x)=1. (For this method to work, both f(x) and g(x) must be twice differentiable.)

> f:=x->x^2; g:=x->1;
a:=1;
u:=(x,t)->(f(x+a*t)+f(x-a*t))/2+(2*a)^(-1)*int(g(z),z=(x-a*t)..(x+a*t)):
u(x,t);

> with(plots):

> plot([u(x,0),u(x,.5),u(x,1),u(x,1.5),u(x,2),u(x,2.5),u(x,3)],x=0..3);
Bottom red curve is y=u(x,0).

> plot3d(u(x,t),x=0..3,t=0..3,axes=FRAME);

Another example. This time, let's try the same functions as in the series example above. Note: f(x) in this case is discontinuous so the theory does not apply here.

> f:=x->piecewise(x<Pi/2,-1,x>Pi/2,2); g:=x->1;
a:=1;
u:=(x,t)->(f(x+a*t)+f(x-a*t))/2+(2*a)^(-1)*int(g(z),z=(x-a*t)..(x+a*t)):
u(x,t);
u(x,0);

> plot([seq(u(x,n*0.2),n=0..15)],x=0..Pi);
The bottom red line is y=u(x,0)

> plot3d(u(x,t),x=0..3,t=0..3,axes=FRAME);

>