# Numerical solution of PDE does not satisfy initial conditions

I want to solve the PDE for the variable $u(r,z,t)$

$$\frac{\partial^2 u}{\partial t^2} = a(r) \left(\frac{\partial^2 u}{\partial r^2} + \frac{1}{r}\frac{\partial u}{\partial r} - \frac{1}{r^2}u + \frac{\partial^2 u}{\partial z^2} \right)$$

subject to the initial conditions $u(r,z,0) = re^{-(r^2+z^2)}$, $\frac{\partial u}{\partial t}(r,z,0) = 0$ and the boundary condition $u(0,z,t) = 0$.

The function $a(r)$ has the form of a step function.

This PDE describes an MHD wave propagating along an infinite cylinder embedded in an infinite medium, hence the physical domain is $0\leq r\leq\infty$, $-\infty\leq z\leq\infty$. Of course, I want to use a smaller numerical domain, e.g., $0\leq r\leq 4$, $-4\leq z\leq 4$. Although in this example I consider $0\leq t\leq 2$, where $t$ is time, I would like to solve the PDE for much longer times.

The code below is one of a few attempts I have tried to solve this problem with NDSolveValue and NDSolve. When I plot the numerically obtained $u$ for fixed $r$ and $z$ as a function of $t$, I obtain solutions that I know to be wrong, so my guess is that I am not correctly using NDSolveValue and NDSolve.

In the code below I simply plot the numerical solution and its error at $t=0$ along the line $r=r0$ and the line $z=z0$. There is a non-negligible difference between the numerical solution and the initial condition. In addition, I am quite surprised to find that this error changes if I modify the numerical domain in which the PDE is solved.

I will thank any ideas on how to improve the accuracy of the solution.

PDE = Derivative[0, 0, 2][u][r, z, t] ==
a[r] (Derivative[2, 0, 0][u][r, z, t] +
1/r Derivative[1, 0, 0][u][r, z, t] - 1/r^2 u[r, z, t] +
Derivative[0, 2, 0][u][r, z, t]);
f[r_, z_] = r Exp[-(r^2 + z^2)];
ic = {u[r, z, 0] == f[r, z], Derivative[0, 0, 1][u][r, z, 0] == 0};
bc = u[0, z, t] == 0;

a[r_] = UnitStep[r - 1] + 1;

rmin = 0;
rmax = 4;
zmin = -4;
zmax = -zmin;
tmax = 2;

usol = NDSolveValue[{PDE, ic, bc},
u, {r, rmin, rmax}, {z, zmin, zmax}, {t, 0, tmax},
Method -> {"MethodOfLines", "TemporalVariable" -> t,
"SpatialDiscretization" -> "FiniteElement"},
InterpolationOrder -> All];

r0 = 1;
z0 = 0;
{Plot[{usol[r0, z, 0], f[r0, z]}, {z, -5, 5}, AxesLabel -> {"z", "u"}],
Plot[f[r0, z] - usol[r0, z, 0], {z, -5, 5},
AxesLabel -> {"z", "Error"}, PlotRange -> All]}
{Plot[{usol[r, z0, 0], f[r, z0]}, {r, 0, 4},
AxesLabel -> {"r", "u"}],
Plot[f[r, z0] - usol[r, z0, 0], {r, 0, 4},
AxesLabel -> {"r", "Error"}, PlotRange -> All]}


Error at t=0 for r=0 using original code Error at t=0 for r=0 using "MaxCellMeasure" -> 0.01 • Can you show the plot of the error? What version are you using? – user21 Oct 18 '17 at 23:12
• Have tried refining the mesh: "SpatialDiscretization" -> {"FiniteElement", "MeshOptions" -> {"MaxCellMeasure" -> 0.01}} – user21 Oct 19 '17 at 2:17
• I am using version 11.1.1.0 – Ramon Oliver Oct 19 '17 at 8:57
• @user21 I am trying to figure out how to insert an image, or a link to it, in my reply. This is my first post in stackexchange. The plots of the errors I get with the original code and with the MaxCellMeasure option (I have used two values: 0.01, as you suggested, and 0.001) show clearly that the MaxCellMeasure option helps reduce the error of the solution at t=0 by 1 or 2 orders of magnitude. – Ramon Oliver Oct 19 '17 at 10:37
• You can add images by editing your post and select the small image icon - right, next to the {} braces for code. Are you satisfied with the MeshCellMeasure as a solution to your issue? – user21 Oct 19 '17 at 10:53

One way to improve this is to use a finer mesh via an option like:

"SpatialDiscretization" -> {"FiniteElement", "MeshOptions" -> {"MaxCellMeasure" -> 0.01}}


Since you have a discontinuity at r==1 you could also add

{"MaxCellMeasure" -> 0.01, "IncludePoints"->{{1,0}}}


To force a line in the mesh at r==1 Play a bit with that and see if it helps.

• MaxCellMeasure works perfect for my purposes. Using IncludePoints, on the other hand, does not make any difference on the solution. Perhaps Mathematica "learns" that some of the PDE coefficients are discontinuous there. Thanks!! – Ramon Oliver Oct 20 '17 at 8:53