37WvEqu.mws

Partial Differential Equations PowerTool

by Dr. Jim Herod

Section 8.3 Numerical Solutions for the One Dimensional Wave Equation

Maple Packages used in Section 8.3

>    restart;

>    with(plots): 

Warning, the name changecoords has been redefined

>    with(LinearAlgebra):

>   

The problem we consider is

         diff(u,`$`(x,2))  = diff(u,`$`(t,2))  - F(t,x)

with boundary conditions      u(t, 0) = u(t,1) = 0,

and with initial conditions      u(0, x) = f(x) and diff(u,t) (0,x) = g(x).

By now, we expect to create a sequence of vectors U[1], U[2] , ... U[M]  along the time axis. Each U[m]  has its values created from what has come before. These U 's are found from the preceding ones in a manner based on equations found by replacing the partial derivatives with appropriate approximations.

If the vector U[m]  has N  terms, then the first and last term are determined by the boundary conditions. That is, for each m, U[m] (0) and U[m] (N) are zero for this problem. In a similar manner, we expect that the vector U[0] , corresponding to t  = 0, is determined by the initial condition, u(0,x) = f(x). The implication is that, if x[n]  is the mesh along the x  axes, then

            U[0] (n) = f( x[n] ) for n = 1, 2, ... , N.

How the remaining initial condition enters the approximation scheme will be developed shortly.

We state the discrete replacement terms: Let h = 1/N and k = 1/M.

     diff(u,`$`(x,2))    is replaced by     (U[m](n+1)-2*U[m](n)+U[m](n-1))/(h^2)

and

     diff(u,`$`(t,2))  - F(t,x)  is replaced by    (U[m+1](n)-2*U[m](n)+U[m-1](n))/(k^2)   - F[m](n)

From the statement of the problem, these should be equal. Let delta  = h/k  .  We let Maple solve this equality for us.

>    eq:=(U[m](n+1)-2*U[m](n)+U[m](n-1))/h^2 =
      (U[m+1](n)-2*U[m](n)+U[m-1](n))/k^2-F[m](n);

eq := (U[m](n+1)-2*U[m](n)+U[m](n-1))/h^2 = (U[m+1](n)-2*U[m](n)+U[m-1](n))/k^2-F[m](n)

>    solve(eq,U[m+1](n));

(k^2*U[m](n+1)-2*k^2*U[m](n)+k^2*U[m](n-1)+2*h^2*U[m](n)-h^2*U[m-1](n)+F[m](n)*h^2*k^2)/h^2

>   

 Here is the Summary Equation:

   U[m+1](n)  =     delta^2*U[m](n-1)+2*(1-delta^2)*U[m](n)+delta^2*U[m](n+1)-U[m-1](n)+k^2*F[m](n) .

This equation is appropriate for m > 1 and for 0 < n < N. We have already discussed what happens for U[0] . Here is the story on U[1] . There are two considerations. From the Summary Equation we see that what should happen is that

    U[1](n)  =     delta^2*U[0](n-1)+2*(1-delta^2)*U[0](n)+delta^2*U[0](n+1)-U[-1](n)+k^2*F[0](n) .

What must be done with that U[-1]  term comes from the remaining initial condition. The discrete replacement term we use is

            (U[1](n)-U[-1](n))/(2*k)   =  g(n / N).

We now know the replacement for the minus one term in the equation defining U[1] :

     U[1](n)  =     delta^2*U[0](n-1)+2*(1-delta^2)*U[0](n)+delta^2*U[0](n+1)-(U[1](n)-2*k*g(n/N))+k^2*F[0](n)  .

Or,

     U[1](n)  = delta^2/2*U[0](n-1)+(1-delta^2)*U[0](n)+delta^2/2*U[0](n+1)+k*g(n/N)+k^2*F[0](n)/2 .

We write these replacement equations in the matrix form. For exposition, suppose that N = 4 = M. Here is the matrix A  that will be used in this analysis.

>    N:=4; M:=4;
delta:=(1/N^2)/(1/M^2);

N := 4

M := 4

delta := 1

>    A:=Matrix([[0,0,0,0,0],[delta^2,1-delta^2,delta^2,0,0],
     [0,delta^2,1-delta^2,delta^2,0],
     [0,0,delta^2,1-delta^2,delta^2],[0,0,0,0,0]]);

A := Matrix(%id = 19660964)

>   

The vector U[0]  is defined using only one of the initial conditions:

>    U[0]:=Vector([f(0),f(1/4),f(2/4),f(3/4),f(1)]);

U[0] := Vector(%id = 19504448)

The vector U[1]  is defined using the above analysis. We first take care of the non-homogeneous term, F(t,x),  with a sequence of vectors V.

>    for m from 1 to M do
  V[m]:=Vector([seq(F(m/M,n/N),n=0..N)]):
od:

Here is U[1] .

>    U[1]:=A.U[0]+k^2*V[1]+2*k*Vector([seq(g(n/N),n=0..N)]);

U[1] := k^2*Vector(%id = 18927984)+2*k*Vector(%id = 19779916)+Vector(%id = 18501596)

The remaining U 's are determined by the iteration

                 U[m+1]  = A*U[m]+V[m+1] .

A word should be said about the size of delta  for which the discrete replacement scheme can be expected to be appropriate. Reference is made to the Courant-Friedrichs-Lewy Convergence Criterion. In this problem, the criteria is that

                                 0 < delta <= 1 .

Problem 1

We work the following problem

         diff(u,`$`(x,2))  = diff(u,`$`(t,2))  - 16 sin( Pi  t),

with boundary conditions      u(t, 0) = u(t,1) = 0,

and with initial conditions      u(0, x) = 0 and diff(u,t) (0,x) = 0.

First, we draw the graph with Maple' s built-in numerical techniques.

>    pde1:=diff(u(t,x),x,x)=diff(u(t,x),t,t)-16*sin(Pi*t);

pde1 := diff(u(t,x),`$`(x,2)) = diff(u(t,x),`$`(t,2))-16*sin(Pi*t)

Here are the boundary conditions and initial condition.

>    bc1:={u(t,0)=0,u(t,1)=0,u(0,x)=0,D[1](u)(0,x)=0};

bc1 := {u(t,0) = 0, u(t,1) = 0, u(0,x) = 0, D[1](u)(0,x) = 0}

Now, we call the numerical solver for this pde. The output is a module. We recall how to make a graph of the solution.

>    u1:=pdsolve(pde1,bc1,type=numeric,time=t,range=0..1,timestep=1/100,spacestep=1/6);

u1 := module () local INFO; export plot, plot3d, animate, value, settings; option `Copyright (c) 2001 by Waterloo Maple Inc. All rights reserved.`; end module

>    u1:-plot3d(u(t,x),t=0..2,axes=normal,orientation=[25,65],
      style=wireframe);

[Maple Plot]

Now, we construct the discrete replacement scheme and compare with the above.

>    f:=x->0:
g:=x->0:

>    F:=(t,x)->16*sin(t*Pi);

F := proc (t, x) options operator, arrow; 16*sin(t*Pi) end proc

>    N:=4;
M:=4;
k:=1/M;
delta:=(1/N)/(1/M);

N := 4

M := 4

k := 1/4

delta := 1

We check that delta  is small enough.

>    is(delta > 0); is(delta<=1);

true

true

Here is the matrix that is used to advance the time steps.

>    A:=Matrix([[0,0,0,0,0],[delta^2,1-delta^2,delta^2,0,0],
     [0,delta^2,1-delta^2,delta^2,0],
     [0,0,delta^2,1-delta^2,delta^2],[0,0,0,0,0]]);

A := Matrix(%id = 3267436)

We need to input the initial conditions.

>    U[0]:=Vector([seq(f(n/N),n=0..N)]);
GINT:=Vector([seq(g(n/N),n=0..N)]);

U[0] := Vector(%id = 18940088)

GINT := Vector(%id = 18930928)

Here is the non-homogeneous term.

>    for m from 0 to 2*M do
  VF[m]:=Vector([0,seq(F(m/M,n/N),n=0..N-2),0]):
od:

>    U[1]:=A.U[0]+k^2*VF[0]/2+k*GINT;

U[1] := Vector(%id = 19139972)

>    for m from 2 to 2*M do
  U[m]:=A.U[m-1]-U[m-2]+k^2*VF[m-1]:
od:

>    for m from 1 to 2*M do J[m]:=pointplot3d([seq([(n-1)/N,m/(M),U[m][n]],n=1..N+1)],axes=normal,connect=true,thickness=3,color=BLACK,orientation=[25,65]):
od:

>    Jlim:=u1:-plot3d(u(t,x),t=0..2,axes=normal,orientation=[25,65],
      style=wireframe):

>    display3d([seq(J[p],p=1..2*M),Jlim]);

[Maple Plot]

This makes a good fit to the graph Maple found with the built-in numerical methods.

>   

>   

Problem 2

The next problem represents a string at rest and tied down at one end. The action is that we move the left end up and down and watch as the wave runs through the string.

>    pde1:=diff(u(t,x),x,x)=diff(u(t,x),t,t);

pde1 := diff(u(t,x),`$`(x,2)) = diff(u(t,x),`$`(t,2))

Here are the boundary conditions and initial condition.

>    bc1:={u(t,0)=sin(Pi*t),u(t,1)=0,u(0,x)=0,D[1](u)(0,x)=0};

bc1 := {u(t,0) = sin(Pi*t), u(t,1) = 0, u(0,x) = 0, D[1](u)(0,x) = 0}

Now, we call the numerical solver for this pde. The output is a module. We recall how to make a plot.

>    u1:=pdsolve(pde1,bc1,type=numeric,time=t,range=0..1,timestep=1/100,spacestep=1/6);

u1 := module () local INFO; export plot, plot3d, animate, value, settings; option `Copyright (c) 2001 by Waterloo Maple Inc. All rights reserved.`; INFO := INFO; plot := proc () if nargs = 0 then error...

>    u1:-plot3d(u(t,x),t=0..4,axes=normal,orientation=[25,65],
      style=wireframe);

[Maple Plot]

An animation really shows what is happening.

>    u1:-animate(u(t,x),t=0..5,frames=30,
              labels=["x","u(x,t)"], labelfont=[TIMES,ROMAN,14]);

[Maple Plot]

Now, we follow methods such as we have seen to make a numerical approximation.

>    f:=x->0:
g:=x->0:

>   

>    F:=(t,x)->0;

F := 0

>    N:=4;
M:=4;
k:=1/M;
delta:=(1/N)/(1/M);

N := 4

M := 4

k := 1/4

delta := 1

>    is(delta > 0); is(delta<=1);

true

true

>    A:=Matrix([[0,0,0,0,0],[delta^2,1-delta^2,delta^2,0,0],
     [0,delta^2,1-delta^2,delta^2,0],
     [0,0,delta^2,1-delta^2,delta^2],[0,0,0,0,0]]);

A := Matrix(%id = 16466356)

>    U[0]:=Vector([seq(f(n/N),n=0..N)]);
GINT:=Vector([seq(g(n/N),n=0..N)]);

U[0] := Vector(%id = 3276620)

GINT := Vector(%id = 19731192)

>    for m from 0 to 2*M do
  VF[m]:=Vector([0,seq(F(m/M,n/N),n=0..N-2),0]):
od:

>    Vector([sin(1/4*Pi),seq(0,i=1..N)]);

Vector(%id = 19258728)

>    U[1]:=A.U[0]+k^2*VF[0]/2+k*GINT+
           Vector([sin(1/4*Pi),seq(0,i=1..N)]);

U[1] := Vector(%id = 16724940)

>    for m from 2 to 2*M do
  U[m]:=A.U[m-1]-U[m-2]+k^2*VF[m-1]+
           Vector([sin(m/M*Pi)+U[m-2][1],seq(0,i=1..N)]):
od:

>    for m from 1 to 2*M do J[m]:=pointplot3d([seq([(n-1)/N,m/(M),U[m][n]],n=1..N+1)],axes=normal,connect=true,thickness=3,color=BLACK,orientation=[25,65]):
od:
m:='m':

>    Jlim:=u1:-plot3d(u(t,x),t=0..2,axes=normal,orientation=[25,65],
      style=wireframe):

>    display3d([seq(J[p],p=1..2*M),Jlim]);

[Maple Plot]

>   

Again, this seems to be a good fit to the graph found by Maple's numerical methods.

>   

>   

Problem 3

Consider this third partial differential equation.

   

         diff(u,`$`(x,2))  = diff(u,`$`(t,2))  + u(t,x)

with boundary conditions      u(t, 0) = u(t,1) = 0,

and with initial conditions      u(0, x) = 1 - 2 | x - 1/2 |  and diff(u,t) (0,x) = 0.

 

This is a problem for which we can draw a graph of the solution with three different methods:

   1. Fourier Series methods will provide an analytic solution for this problem.

   2. The Discrete Replacement Methods will provide an approximation for the graph of the solution for this problem.

   3. Maple's PDSolve routine will not only provide an approximation for the graph, it can also provide an animation for the solution.

Concerning the first method, we have computed the Fourier solution. The first five terms of  the infinite series are

>    w:=(t,x)->8/Pi^2*sum((-1)^(n+1)/(2*n-1)^2*
   cos(sqrt(1+(2*n-1)^2*Pi^2)*t)*sin((2*n-1)*Pi*x),n=1..5);

w := proc (t, x) options operator, arrow; 8/Pi^2*sum((-1)^(n+1)/(2*n-1)^2*cos(sqrt(1+(2*n-1)^2*Pi^2)*t)*sin((2*n-1)*Pi*x),n = 1 .. 5) end proc

Concerning the second method, note that the replacement equations take only a slightly different form.

>    eq:=(U[m](n+1)-2*U[m](n)+U[m](n-1))/H^2 =
      (U[m+1](n)-2*U[m](n)+U[m-1](n))/K^2+U[m](n);

eq := (U[m](n+1)-2*U[m](n)+U[m](n-1))/H^2 = (U[m+1](n)-2*U[m](n)+U[m-1](n))/K^2+U[m](n)

>    solve(eq,U[m+1](n));

(K^2*U[m](n+1)-2*K^2*U[m](n)+K^2*U[m](n-1)+2*H^2*U[m](n)-H^2*U[m-1](n)-U[m](n)*H^2*K^2)/H^2

It is unlikely that programming this iteration scheme will add significantly to the ideas of this section.

However, it does seem well to see an animation of how the string vibrates. We use Maple's built in PDE solver to draw this. Already, we know how to get a graph for the solution.

>    pde1:=diff(u(t,x),x,x)=diff(u(t,x),t,t)+u(t,x);

Here are the boundary conditions and initial condition.

pde1 := diff(u(t,x),`$`(x,2)) = diff(u(t,x),`$`(t,2))+u(t,x)

>    bc1:={u(t,0)=0,u(t,1)=0,u(0,x)=1-2*abs(x-1/2),D[1](u)(0,x)=0};

bc1 := {u(0,x) = 1-2*abs(x-1/2), u(t,0) = 0, u(t,1) = 0, D[1](u)(0,x) = 0}

Now, we call the numerical solver for this pde. The output is a module. We illustrate how to make a plot.

>    u1:=pdsolve(pde1,bc1,type=numeric,time=t,range=0..1,timestep=1/100,spacestep=1/20);

u1 := module () local INFO; export plot, plot3d, animate, value, settings; option `Copyright (c) 2001 by Waterloo Maple Inc. All rights reserved.`; INFO := INFO; plot := proc () if nargs = 0 then error...

>    u1:-plot3d(u(t,x),t=0..5,axes=normal,orientation=[-35,65]);

[Maple Plot]

>    u1:-animate(u(t,x),t=0..5,frames=30,
       labels=["x","u(x,t)"], labelfont=[TIMES,ROMAN,14],
       scaling=constrained);

[Maple Plot]

It has been the intent to show the basic methods for numerical methods for solving wave equations. But also, for the user, it is good to see the increased utility of using the numerical methods built into Maple.

>   

>   

>   

>   

EMAIL: herod@math.gatech.edu   or   jherod@tds.net

URL: http://www.math.gatech.edu/~herod

Copyright ©  2003  by James V. Herod

All rights reserved