How to solve time dependent pde?

Question

I have set out to try and solve a time dependent wave equation with specific boundary conditions, but I am running into a bit of a snag. As you can see, I am trying to solve for the wave function u[t, x] traveling at velocity 1, with my boundary conditions being a fixed point at u[t, 30], and a Gaussian wave packet at t = 0.

    difeqns = {D[u[t, x, {t, 2}] - D[u[t, x], {x, 2}] ==  0, 
    u[t, 30] == 0, 
    u[0, x] == E^(-x^2), 
    Derivative[1, 0][u][0, x]== 0}

I am plotting my solution as a variable of time holding n constant at 20 using the code:

    solv = NDSolve[difeqns, u[t, x], {t, 0, 500}, 
    {x, -30, 30},  
    WorkingPrecision -> 5, InterpolationOrder -> All];  
    x= 20
    Plot[u[t, x]/. solv, {t, 0, 500}, PlotRange -> All]

Since there is a fixed point set at u[t, 30], I expect the initial Gaussian wave packet to reflect off of the fixed point on the right side and then continue onto negative infinity, so as time goes on, no oscillation at all should occur. However, my oscillation is not going to zero at x=20 as time increases, and when I animated my solution using ListAnimate, it seemed like the wave packet was reflecting off of a boundary on the left side even though no such boundary should exist. What have I done wrong in my numerical approach to this problem?

Answer 1

Summary

If formulated consistently the problem can be solved analytically. It is completely conservative. Hence there are no decaying oscillations but just a pair of profiles moving, with velocity +- 1, in opposite directions and being reflected at the symmetric boundary.

The numerical problem of the OP is mainly due to due to the low WorkingPrecision; for small "a" the exact solution is susceptible to the inconsitency between the boundary condition and the intial condition, as the latter does not obey the boundary condition. Numerically this effect is less pronounced.

In the following we present first the exact and then the numerical solution.

Exact solution

The solution procedure is completely standard but it will be reprocuced here for clarity. Making a separation ansatz we find particular solutions (eigenfunction) of the form

uc[k_,t_,x_,a_] := Cos[ x/a Pi (1/2 + k) ] Cos[ t/a Pi (1/2 + k) ]

These satisfy the boundary condition

(1) u[t,x=a] = 0

and the intial condition

(2) D[u[t,x],t]/.t->0 = 0

The general solution is a superposition of the eigenfunctions

u[t_, x_] := 
 Sum[c[k] Cos[x/a \[Pi] (1/2 + k)] Cos[t/a \[Pi] (k + 1/2)], {k, 
   0, \[Infinity]}]

the coefficients c[k] must be determined from the initial condition

(3) u[ t=0, x ] = f[x]

where f[x] is assumed to be symmetric (f[-x] = f[x]).

Notice that the other possible eigenfunction

u1[k_,t_,x_,a_] := Sin[ x/a Pi k] Cos[ t/a Pi k ]

is ruled out by the symmetry of f.

Summarizig up to this point: the solution is symmetric in x which means that it vanishes not only at x = a but also at x = -a.

Now an intial condition consistent with the boundary conditions at x = +- a requires that

f[x = -+a] = 0

This is not the case for the Gaussian profile of the OP. Hence we modify the original initial condition

f[x_]= Exp[-x^2/2];

and set

f1[x_] = Exp[-x^2/2] - Exp[-a^2/2];

The coefficients are found using the orthogonality property

Integrate[
 1/a Cos[x/a \[Pi] (k + 1/2)] Cos[x/a \[Pi] (m + 1/2)], {x, -a, a}, 
 Assumptions -> k == m]
Simplify[%, m \[Element] Integers]

(* Out[598]= 1 - Sin[2 m \[Pi]]/(\[Pi] + 2 m \[Pi])

Out[599]= 1 *)

Integrate[
 1/a Cos[x/a \[Pi] (k + 1/2)] Cos[x/a \[Pi] (m + 1/2)], {x, -a, a}, 
 Assumptions -> k != m]
Simplify[%, {k, m} \[Element] Integers]

(* Out[603]= ((1 + 2 m) Cos[m \[Pi]] Sin[k \[Pi]] - (1 + 2 k) Cos[
   k \[Pi]] Sin[m \[Pi]])/((k - m) (1 + k + m) \[Pi])

Out[604]= 0 *)

to be

c1[k_] = Integrate[f1[x] 1/a Cos[x/a \[Pi] (k + 1/2)], {x, -a, a}]

(* Out[611]= -((4 E^(-(a^2/2)) Cos[k \[Pi]])/(\[Pi] + 2 k \[Pi])) + (
 I E^(-((\[Pi] + 2 k \[Pi])^2/(8 a^2))) Sqrt[\[Pi]/
  2] (Erfi[(-2 I a^2 + \[Pi] + 2 k \[Pi])/(2 Sqrt[2] a)] - 
    Erfi[(2 I a^2 + \[Pi] + 2 k \[Pi])/(2 Sqrt[2] a)]))/a *)

c[k_] = Integrate[f[x] 1/a Cos[x/a \[Pi] (k + 1/2)], {x, -a, a}]

(* Out[614]= (I E^(-((\[Pi] + 2 k \[Pi])^2/(
  8 a^2))) Sqrt[\[Pi]/2] (Erfi[(-2 I a^2 + \[Pi] + 2 k \[Pi])/(
    2 Sqrt[2] a)] - 
   Erfi[(2 I a^2 + \[Pi] + 2 k \[Pi])/(2 Sqrt[2] a)]))/a *)

Defining the sums with kk terms as

u1[t_, x_, kk_] := 
 Sum[c1[k] Cos[x/a \[Pi] (1/2 + k)] Cos[t/a \[Pi] (k + 1/2)], {k, 0, 
   kk}]

u[t_, x_, kk_] := 
 Sum[c[k] Cos[x/a \[Pi] (1/2 + k)] Cos[t/a \[Pi] (k + 1/2)], {k, 0, 
   kk}]

we can compare the two expressions in a plot. The difference is pronounced only for very small a:

a = 1; t0 = 0; Plot[{u1[t0, x, 5], u[t0, x, 20], f[x], 
  f1[x]}, {x, -1.1 a, 1.1 a}, PlotStyle -> {Red, Blue, Yellow, Green},
  PlotRange -> All, AxesLabel -> {"x", "u(0,x)"}, 
 PlotLabel -> 
  "Comparison of exact solutions with a = 1 at t = 0\nfor consistent \
(red)\nand inconsistent (blue) initial conditions\nAlso the green \
curve shows f[x], the red curve f1[x]\n"]

In what follows we shall discard the inconsistent case and plot the solution in the range -a<x<a for various instants of time

a = 30; Plot[{u1[60, x, 30], u1[55, x, 30], u1[35, x, 30], 
  u1[28, x, 30], u1[3, x, 30], u1[0, x, 30]}, {x, -a, a}, 
 PlotRange -> {-1, 1}, 
 PlotLabel -> 
  "The spatial dependence of the solution for various instants of \
time.\nThe intial yellow peak (t = 0) splits into two blue peaks (t = \
3)\nThese peaks have almost reached the boundary (red, t = 28),\nturn \
into green peaks after reflexion (t = 35),\nand move back to x = 0 as \
yellow peaks (t = 55)\njoining to exactly the negative of the initial \
peak (blue, t = 60)\n"]

Numerical solution

The main part of the code is a modification of the code of the OP in this respect

• symmetric boundary conditions
• inital condition consistent with boundary conditions
• take WorkingPrecision as a parameter

Usage instructions

First enter the parameters and run (solution) Then you can choose to plot either the time dependence at a specific location or the spatial dependence at a specific instant of time using the appropiate parts of the code.

(* solution *)
parm = {a = 5, tmax = 15, wp = 5, x0 = 0, cons = True};
difeqns1 = {D[un[t, x], {t, 2}] - D[un[t, x], {x, 2}] == 0, 
   un[t, a] == 0, un[t, -a] == 0, 
   un[0, x] == Exp[-x^2/2] - Boole[cons]*Exp[-a^2/2], 
   Derivative[1, 0][un][0, x] == 0};
solv = NDSolve[difeqns1, un[t, x], {t, 0, tmax}, {x, -a, a}, 
  WorkingPrecision -> wp, InterpolationOrder -> All]; 


(* plot time dependence *)
Plot[
 un[t, x] /. solv /. x -> x0, {t, 0, tmax}, PlotRange -> All, 
 PlotLabel -> 
  "Numerical solution of a time dependent PDE (1D wave equation)\n\
Time depencence at a specific point\n (x0 = " <> ToString[x0] <> 
   ", a = " <> ToString[a] <> ", tmax = " <> ToString[tmax] <> 
   ",\nWorkingPrecision = " <> ToString[wp] <> ", consistency = " <> 
   ToString[cons] <> ")\n"]

(* plot spatial dependence *)
t0 = 3;
Plot[un[t, x] /. solv /. t -> t0, {x, -a, a}, PlotRange -> All, 
 PlotLabel -> 
  "Numerical solution of a time dependent PDE (1D wave equation)\n\
Spatial depencence at a specific instant\n (t0 = " <> ToString[t0] <> 
   ", a = " <> ToString[a] <> ", tmax = " <> ToString[tmax] <> 
   ",\nWorkingPrecision = " <> ToString[wp] <> ", consistency = " <> 
   ToString[cons] <> ")\n"]

Some typical results are

1) take the region a = 5 (in order to speed up the calculations) and compare the spatial dependence for different values of the WorkingPrecision (wp)

We can see that for wp = 5 there are spurious oscillations which vanish at wp = 10

2) study the time dependence at a specific position for different values of the WorkingPrecision (wp)

For wp = 10 gradually appear spurious additional oscillations which are suppressed for wp = 15.

Finally, here is the case of the OP with a = 20, x0 = 20 for three increasing values of wp (5, 10, and 15).

Only at wp = 15 the accuracy is high enough for the solution to show the expected completely periodic behaviour generated by the two profiles moving with velocity 1 to and fro between the boundaries.

This graph is the correction of the graph in the solution of Young as of Jul 26 at 20:13 up to t = 250.

It is also instructive to study energy conservation. The extent of non-conservation can serve as a measure for numerical accuracy.

The energy is defined as the integral over the energy density as follows

ww[t_] = Integrate[
  1/2 D[uun[t, x], t]^2 + 1/2 D[uun[t, xx], xx]^2 /. xx -> x, {x, -a, 
   a}]

It can easily be shown to be conserved in time with the given boundary conditions.

The initial value is given by

ww0 = 
 Integrate[1/2 D[Exp[-xx^2/2], xx]^2 /. xx -> x, {x, -a, a}] 

(* Out[634]= 1/4 (-2 a E^-a^2 + Sqrt[\[Pi]] Erf[a])  *)

Example a = 5 wp = 5 and 10

After running the solution code we proceed as follows:

The numerical value of the initial energy is

ww0n = ww0[a] // N

(* Out[659]= 0.443113 *)

The energy values at equidistant instances are (replacing the very slow spatial integrals)

wwt = Table[{t, ww[t]}, {t, 0, 15, 0.5}] // N;

for comparison the inital value

wwt0 = Table[{t, ww0n}, {t, 0, 15, 0.5}];

Plotting the result with

ListLinePlot[{wwt0, wwt}, PlotRange -> {0, 1.3}, 
 AxesLabel -> {"t", "energy"}, 
 PlotLabel -> 
  "Numerical solution of a time dependent PDE (1D wave equation)\n\
Time depencence of the energy\n (a = " <> ToString[a] <> ", tmax = " <>
    ToString[tmax] <> ",\nWorkingPrecision = " <> ToString[wp] <> 
   ", consistency = " <> ToString[cons] <> ")\n"]

gives "bad" energy conservation for wp = 5

and almost perfect energy conservation for wp = 10

Answer 2

Clear[x]

solv = NDSolve[{
    D[u[t, x], {t, 2}] - D[u[t, x], {x, 2}] == 0,
    u[t, 30] == 0, u[0, x] == E^(-x^2),
    Derivative[1, 0][u][0, x] == 0},
   u[t, x], {t, 0, 2000}, {x, -30, 30},
   Method -> {"MethodOfLines", 
     "SpatialDiscretization" -> {"TensorProductGrid", 
       "MinPoints" -> 50, "MaxPoints" -> 50, "DifferenceOrder" -> 2},
     Method -> {LSODA, "MaxDifferenceOrder" -> 1}}];

x = 20;
Plot[u[t, x] /. solv, {t, 0, 2000}, PlotRange -> All]