# Singularities forming on boundary while solving system of pde's

Backslide introduced after 9.0.1, fixed in 12.0.1 or earlier.

This is a follow up of a previous question I asked regarding solving a system of coupled, non-linear partial differential equations, 2D spatially + time. The equations (shown below) model a magnetic field diffusing into a plasma column.

Note: the equation is actually in polar coordinates. I made the changes r -> x and theta -> y for simplicity in the code I post below.

Before my previous question, I was having problems solving the system for any measure of time. Thanks to a recommendation, I rescaled the time derivative which helped me solve a solution well past, in time, what I was capable of previously.

Now, however, I'm running into apparent singularities appearing along the boundary x = 0. Frankly, this was probably my problem the whole time and I'm only just now realizing it. I'll illustrate what I mean by discussing two cases.

The equations in code form are:

azExt[x_, y_, t_] = -bo*a*(1 - E^(-f*t*2*Pi/(w*T)))*x*Sin[2*Pi*f*t - y];
eqns = {c*D[az[x, y, t], t] == (Laplacian[az[x, y, t], {x, y}, "Polar"] + 1/x*d*p*(D[az[x, y, t], x]*D[bz[x, y, t], y] - D[az[x, y, t], y]*D[bz[x, y, t], x])),
c*D[bz[x, y, t], t] == (Laplacian[bz[x, y, t], {x, y}, "Polar"] + 1/x*1/d*p* (D[az[x, y, t], y]* D[Laplacian[az[x, y, t], {x, y}, "Polar"], x] - D[az[x, y, t], x]* D[Laplacian[az[x, y, t], {x, y}, "Polar"], y])),
(*Initial Conditions*)
az[x, y, 0] == 0,
bz[x, y, 0] == 1,
(*Boundary Conditions*)
bz[x, Pi, t] == bz[x, -Pi, t],
az[x, Pi, t] == az[x, -Pi, t],
bz[1, y, t] == 1,
az[1, y, t] == 0,
(D[bz[x, y, t], x] /. x -> x0) == 0,
(D[az[x, y, t], x] /. x -> 1) == 2*azExt[1, y, t]/(bo*a) - az[1, y, t],
(D[D[az[x, y, t], x], x] /. x -> 1) == (D[az[x, y, t], x] /. x -> 1)};


Where the constants are:

bo = 0.005;
w = 5*10^6;
T = 0.4*10^-6;
a = 0.02;
f = 4.*10^-2;
t0 = 3.8595/(f)
c = 40/f;
d = 10;
p = 31;
x0 = 10^-6;


Here, f is a scaling factor for the time derivatives, so it can change from solution to solution without affecting the physics of the system.

The equations are solved using NDSolve and plotted using the following code:

sol = NDSolve[eqns, {az, bz}, {x, x0, 1}, {y, -Pi, Pi}, {t, 0, t0}, "ExtrapolationHandler" -> {Indeterminate &, "WarningMessage" -> False}, Method -> {"MethodOfLines", "SpatialDiscretization" -> {"TensorProductGrid", "MinPoints" -> 20, "MaxPoints" -> 20, "DifferenceOrder" -> 2}}]

{{ContourPlot[Evaluate[az[x, 0, t] /. sol], {x, x0, 1}, {t, 0, t0}, Contours -> 20, ColorFunction -> Hue, PlotLegends -> Automatic, PlotRange -> All, FrameLabel -> {"x", "t"}, PlotLabel -> "Az"],
ContourPlot[Evaluate[bz[x, 0, t] /. sol], {x, x0, 1}, {t, 0, t0}, Contours -> 20, ColorFunction -> Hue, PlotLegends -> Automatic, PlotRange -> All, FrameLabel -> {"x", "t"}, PlotLabel -> "Bz"]},
{ContourPlot[Evaluate[az[x, y, t0] /. sol], {x, x0, 1}, {y, -Pi, Pi}, Contours -> 20, ColorFunction -> Hue, PlotLegends -> Automatic, PlotRange -> All, FrameLabel -> {"x", "y"}, PlotLabel -> "Az"],
ContourPlot[Evaluate[bz[x, y, t0] /. sol], {x, x0, 1}, {y, -Pi, Pi}, Contours -> 20, ColorFunction -> Hue, PlotLegends -> Automatic, PlotRange -> All, FrameLabel -> {"x", "y"}, PlotLabel -> "Bz"]}}


Now, here's the odd thing. If I specify t0 = 3.8594/f (or 96.485) , the solution converges with pretty much no problem, just a eerr warning about spatial errors. I get something that looks like this, which makes physical sense:

Great! Now, if I increase t0 to 3.8595/f, a few things happen. For one, the solution is not able to converge. In fact, it stops at t = 73.18 (as opposed to 96.4875). This is odd as in the previous example, the solution was able to solve well past this time stamp. Further more, a singularity forms at t = 73.18 as shown below:

As you can see, the singularity begins forming along the boundary near x = 0 and I get the following errors:

NDSolve::ndsz: At t == 73.18237747079812, step size is effectively zero; singularity or stiff system suspected.
NDSolve::eerr: Warning: scaled local spatial error estimate of 29.87465573060914 at t = 73.18237747079812 in the direction of independent variable x is much greater than the prescribed error tolerance. Grid spacing with 20 points may be too large to achieve the desired accuracy or precision. A singularity may have formed or a smaller grid spacing can be specified using the MaxStepSize or MinPoints method options.


My question is, how can I either avoid or rectify this singularity? I was thinking that perhaps I could specify a constant region for Az and Bz near the boundary. I've tried increasing the max points but I either get the singularity at a different time and/or I get large (>10^3) spatial errors, which also seems odd to me since I would think refining the grid would decrease the spatial errors while also helping to tackle the singularity problem. I'm at a loss here as to what to do. I'm still fairly new at Mathematica and solving PDE's numerically so any help would be greatly appreciated. Please let me know if I missed anything in my code or if something doesn't make sense to you. Thank you!

• You introduced the singularities yourself by using polar coordinates: While the actual PDE on the disk may have smooth solutions, the polar coordinates pulls this solution back to something with a singularity for $r= 0$ because the differential of the polar coordinate mapping is singular for $r=0$. So any discrete method for solving the PDE will exhibit large errors because the fundamental assumptions in all convergence proofs for discrete methods exploit that the solution of the continuous problem is smooth. I'd suggest to apply FEM (by using NDSolve) directly on the original disk. – Henrik Schumacher Jan 18 '19 at 15:12
• @HenrikSchumacher Sorry, I should have clarified. Instead of using r = 0, which I realize will introduce a singularity, I use a small value. In this case, 10^-6, to avoid any infinite expressions. Is that what you are refering to? Also, from my understanding FEM does not work for spatial derivatives with an order higher than 2. Mine has spatial derivatives up to 3. – Woody Jan 18 '19 at 15:25
• @xzczd Sorry, should say "increase t0 to 3.8595/f". I've edited the question. – Woody Jan 18 '19 at 15:38
• Just tested in v12.0.1, the ndsz warning no longer pops up, at least in this specific example. – xzczd Aug 29 '20 at 3:59

The problem is version related. I manage to reproduce the warning in v11.2, but not in v9.0.1. So the ndsz warning may be related to improper time step choosing in v11. Following this guess, I find MaxStepSize -> {Automatic, Automatic, 0.1} helps:

sol = NDSolve[eqns, {az, bz}, {x, x0, 1}, {y, -Pi, Pi}, {t, 0, t0},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid", "MinPoints" -> 20,
"MaxPoints" -> 20, "DifferenceOrder" -> 2}},
MaxStepSize -> {Automatic, Automatic, 0.1}]; // AbsoluteTiming


The first 2 Automatic is the grid size specification for x and y. This isn't the only solution. The fix function also helps:

(* Definition of fix isn't included in this post,
`