# NDSolve cannot evolve Gaussian wave packet with small width

I want to solve a 1d wave equation with initial conditions: $$u(x,0) = f(x) = e^{-(x - x0)^2/2\sigma^2} \ \ \text{and}\\ u_t(x,0) = - f_x(x)$$ The second condition is the advection equation and describes a wave packet moving to the right.

    xmin = 0; xmax = xmin + 50; tmax = 100; σ = 5; x0 = (xmin + xmax)/2;
f[x_] := Exp[-(x - x0)^2/(2*σ^2)];
g[x_] := (E^(-((x - x0)^2/(2 σ^2))) (x - x0))/σ^2;
weqn = D[u[t, x], {t, 2}] == D[u[t, x], {x, 2}];

uifWave =
NDSolveValue[{weqn, u[0, x] == f[x],
Derivative[1, 0][u][0, x] == g[x],
DirichletCondition[u[t, x] == f[xmin], x == xmin],
DirichletCondition[u[t, x] == f[xmax], x == xmax]},
u, {t, 0, tmax}, {x, xmin, xmax}, PrecisionGoal -> 100]


I am also imposing Dirichlet boundary conditions and setting the function to $$0$$ at both ends.

This first animation was made by choosing $$\sigma = 5$$, $$x_{min} = 0$$, and $$x_{max} = 50$$. It behaves as expected.

However, if the width of the pulse is shorter (for instance, $$\sigma = 1$$), the solution is simply numerical error. I have tried setting PrecisionGoal->100 to no avail.

$\sigma = 1$">

Is there anything I can do to make it work for shorter widths?

• Try to use a MaxCellMeasure to specify a finer mesh. – user21 May 29 at 18:56
• @user21 I put Method -> {"FiniteElement", "MeshOptions" -> MaxCellMeasure -> 0.001} inside NSDolve and I am getting the error InitializeBoundaryConditions::fembdnl How should I do it properly? Thank you. – Thiago May 30 at 0:45
• @Thiago This is not a mathematical definition It behaves as expected. For the wave equation, you can put three types of boundary conditions: periodic, with reflection, with absorption. You probably want to get a wave with a reflection. – Alex Trounev May 30 at 5:32

You can use this:

(*xmin=0;xmax=xmin+50;tmax=100;\[Sigma]=5;x0=(xmin+xmax)/2;*)
xmin = \
0; xmax = xmin + 50; tmax = 100; \[Sigma] = 1; x0 = (xmin + xmax)/2;

f[x_] := Exp[-(x - x0)^2/(2*\[Sigma]^2)];
g[x_] := (E^(-((x - x0)^2/(2 \[Sigma]^2))) (x - x0))/\[Sigma]^2;
weqn = D[u[t, x], {t, 2}] == D[u[t, x], {x, 2}];

uifWave =
NDSolveValue[{weqn, u[0, x] == f[x],
Derivative[1, 0][u][0, x] == g[x],
DirichletCondition[u[t, x] == f[xmin], x == xmin],
DirichletCondition[u[t, x] == f[xmax], x == xmax]},
u, {t, 0, tmax}, {x, xmin, xmax},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"FiniteElement",
"MeshOptions" -> {MaxCellMeasure -> 0.1}}}]
Manipulate[
Plot[uifWave[t, x], {x, 0, 50}, PlotRange -> {1, -1}], {t, 0, tmax}]


I also would like to direct you to the Acoustics in the Time Domain tutorial which has lots and lots of information on the wave equation, how to model with it, how to set up boundary condition etc. (The web page does not look good, but the in product page is fine. Paste this into the V12 help system: PDEModels/tutorial/AcousticsTimeDomain)

• Thank you! That was really helpful. – Thiago May 31 at 15:52