# Problem with using NDSolve on wave PDE on string

I am trying to verify my hand solution for wave PDE in 1D, on a string. The solution I get from NDSolve does not match the analytical solution.

This is a string, fixed on left end, and free to move on the right. There is damping involved. (from air motion). The damping is set to be critical for mode 1 only. The string is pulled up for initial position, but has zero initial velocity:

Solving the wave PDE on string $$t>0,0 $$u_{tt}+bu_{t}=c^{2}u_{xx}\qquad00$$ Boundary conditions \begin{align*} u\left( 0,t\right) & =0\\ \left. \frac{\partial u}{\partial x}\right\vert _{x=L} & =0 \end{align*} Initial conditions, $$t=0$$ \begin{align*} u_{t}\left( x,0\right) & =0\\ u\left( x,0\right) & =f\left( x\right) \end{align*} Using \begin{align*} f\left( x\right) & =\left\{ \begin{array} [c]{ccc} \frac{3h}{L}x & & 0 Hence the PDE becomes $$u_{tt}+\frac{4\pi}{3}u_{t}=16u_{xx}$$.

I get warnings from NDSolve which I need help to get rid of them. I think they are what cause the numerical solution not to look the same as the analytical.

Here is what I tried

ClearAll[x, t, f];
L = 3; c = 4; h = 1/10; b = Pi*c/L;
f[x_] := Piecewise[{{3*h/L*x, 0 < x < L/3}, {h, L/3 < x < L}}];
Plot[f[x], {x, 0, L}]  (*Initial position*)


pde = D[u[x, t], {t, 2}] + b D[u[x, t], t] == c^2* D[u[x, t], {x, 2}];
bc = {u[0, t] == 0, Derivative[1, 0][u][L, t] == 0};
ic = {Derivative[0, 1][u][x, 0] == 0, u[x, 0] == f[x]};


Now, if I just type

numericalSol = NDSolve[{pde, ic, bc}, u, {x, 0, L}, {t, 0, 2}]


It hangs and I get

Using

numericalSol =
NDSolve[{pde, ic, bc}, u, {x, 0, L}, {t, 0, 2},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> 200}}]


It does not hang, and gives solution with these warnings

But the solution when animated does not look correct:

Manipulate[
Quiet@Chop@Plot[Evaluate[u[x, t] /. numericalSol], {x, 0, L},
BaseStyle -> 15,
ImageMargins -> 3,
PerformanceGoal -> "Speed",
PlotRange -> {{0, L}, {-0.1, 0.11}},
ImageSize -> 500,
AxesLabel -> {"x", "U"}
],
{{t, 0, "time"}, 0, 2, 0.01}
]


Here is the animation from the analytical solution to compare with

I understand that NDSolve gives Warning: boundary and initial conditions are inconsistent. but I do not see why, as I think they are consistent. Free right end mean $$u_x(L)=0$$ and that is what I used.

Note that DSolve can't solve this PDE in V 12:

pde = D[u[x, t], {t, 2}] + b D[u[x, t], t] == c^2* D[u[x, t], {x, 2}];
bc = {u[0, t] == 0, Derivative[1, 0][u][L, t] == 0};
ic = {Derivative[0, 1][u][x, 0] == 0, u[x, 0] == f[x]};
DSolve[{pde, ic, bc}, u[x, t], {x, t}]


(no solution)

Question is: What option to use to make NDSolve remove these warning with the hope the result will look like my solution (I think the analytical solution is correct, but there is of course possibility I made error, that is why I was trying to verify it with NDSolve)

This problem is related to this one. The definition for f is improper, because the default value of Piecewise is 0, which leads to $$f(L)=0$$. Modifying it to
f[x_] = Piecewise[{{3 h/L x,  x < L/3}}, h]

will remove the ibcinc warning and make NDSolve produce the expected result. (It's still necessary to control the spatial grid size as you've already done, of course. )
• I tried your approach together with "SpatialDiscretization" -> {"FiniteElement",...} but I get a warning CoefficientArrays::poly: (u^(1,0))[3,t] is not a polynomial. What could be the reason? Thanks! – Ulrich Neumann Jul 11 '19 at 8:20
• @UlrichNeumann "FiniteElement" method cannot handle the b.c. Derivative[1, 0][u][L, t] == 0 (at least now), one needs to rewrite it with NeumannValue. (BTW, some related discussion can be found in the comments under this post. ) (Edit: Oh I forgot it's not necessary to explicitly set zero NeumannValue, one can just remove that b.c., as suggested by user21 below. ) – xzczd Jul 11 '19 at 8:28
• @UlrichNeumann, just delete that Derivative[1, 0][u][L, t] == 0 boundary condition. – user21 Jul 11 '19 at 8:29
• @Nasser The default value of Piecewise is 0 i.e. $f(x)=0$ when $x\ge L$ under your definition. – xzczd Jul 11 '19 at 8:33
• @xzczd Thanks, FiniteElement seems to be much more robust (concerning the definition of f[x]) – Ulrich Neumann Jul 11 '19 at 9:43