5
$\begingroup$

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<x<L$ $$ u_{tt}+bu_{t}=c^{2}u_{xx}\qquad0<x<L,t>0 $$ 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<x<\frac{L}{3}\\ h & & \frac{L}{3}<x<L \end{array} \right. \\ b & =\frac{4\pi}{3}\\ c & =4\\ L & =3\\ h & =\frac{1}{10} \end{align*} 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*)

Mathematica graphics

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

Mathematica graphics

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

Mathematica graphics

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}
 ]

enter image description here

Here is the animation from the analytical solution to compare with

enter image description here

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)

$\endgroup$
6
$\begingroup$

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. )

$\endgroup$
  • $\begingroup$ 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! $\endgroup$ – Ulrich Neumann Jul 11 at 8:20
  • 2
    $\begingroup$ @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. ) $\endgroup$ – xzczd Jul 11 at 8:28
  • 1
    $\begingroup$ @UlrichNeumann, just delete that Derivative[1, 0][u][L, t] == 0 boundary condition. $\endgroup$ – user21 Jul 11 at 8:29
  • 1
    $\begingroup$ @Nasser The default value of Piecewise is 0 i.e. $f(x)=0$ when $x\ge L$ under your definition. $\endgroup$ – xzczd Jul 11 at 8:33
  • 1
    $\begingroup$ @xzczd Thanks, FiniteElement seems to be much more robust (concerning the definition of f[x]) $\endgroup$ – Ulrich Neumann Jul 11 at 9:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.