# Error when extending 1-dimensional PDE to 2 dimensions

I want to calculate how magnetic flux is trapped in a superconductor near the interface superconductor/vacuum. This problem already was solved analytically by J. Pearl for cylindrical symmetry (if interested, see: “Structure of Superconductive Vortices near a Metal‐Air Interface” by J. Pearl, 1966). So my idea is, to begin with the numerical solution of this special problem and to compare it with the given analytical solution.

The general PDEs I want to solve are: (first eq. for z<0 (vacuum) and second for z>0 (superconductor)) For z--> +infinity, we are far away from the interface at z=0, and thus the z-dependency of f2 should disappear. Setting Phi0=Lambda=1, f=f2, the second equation reduces to:

pde1 = D[1/r*D[r*f[r], r], r] - f[r] + 1/(2*Pi*r) == 0;


The boundary conditions with the boundaries “start” and “stop” are known from physics:

bc1 = {f[start] == 0, f[stop] == 1/(2*Pi*stop)};


So the first script for calculating the one-dimensional special case z-->+infinity is:

Clear["Global*"];
start = 1/1000;
stop = 20;
wp = 50;
pde1 = D[1/r*D[r*f[r], r], r] - f[r] + 1/(2*Pi*r) == 0;
bc1 = {f[start] == 0, f[stop] == 1/(2*Pi*stop)};
solf = NDSolveValue[{pde1, bc1}, f, {r, start, stop}, WorkingPrecision -> wp];
Plot[solf[r], {r, start, stop}, PlotRange -> All]


The result fits perfectly to the analytical solution!

Now, I want to extend the problem to two dimensions. For this, I replace f(r) by g(r,z), leading to the new pde2 (which has the same structure as pde1):

pde2 = D[1/r*D[r*g[r, z], r], r] - g[r, z] + 1/(2*Pi*r) == 0;


The boundary conditions are generalized to two dimensions in a manner that should produce the same result as in the one-dimensional case, since g[r,z] is forced to reproduce the one-dimensional solution at each border. (Important: This type of boundary conditions can’t be solved by the method of lines and thus requires the FEM-method which is implemented since Mathematica 10, so it doesn’t run in earlier versions!)

bc2 = {g[start, z] == solf[start], g[stop, z] == solf[stop], g[r, zborder] == g[r, -zborder] == solf[r]};


For each value of z, the result found by NDSolve should therefore be equivalent to the one-dimensional case (--> invariant in z-direction). The script for the two-dimensional extension is:

Clear["Global*"];
start = 1/1000;
stop = 20;
zborder = 10;
wp = 50;
pde1 = D[1/r*D[r*f[r], r], r] - f[r] + 1/(2*Pi*r) == 0;
pde2 = D[1/r*D[r*g[r, z], r], r] - g[r, z] + 1/(2*Pi*r) == 0;
bc1 = {f[start] == 0, f[stop] == 1/(2*Pi*stop)};
bc2 = {g[start, z] == solf[start], g[stop, z] == solf[stop], g[r, zborder] == g[r,-zborder] == solf[r]};
solf = NDSolveValue[{pde1, bc1}, f, {r, start, stop}, WorkingPrecision -> wp];
solg = NDSolveValue[{pde2, bc2}, g, {r, start, stop}, {z, -zborder, zborder},    WorkingPrecision -> wp];
Plot3D[solg[r, z], {r, start, stop}, {z, -zborder, zborder}, PlotPoints -> {100, 30}]
Plot[{solg[r, 0], solf[r]}, {r, start, stop}, PlotRange -> All]  For whatever reason, solg[r,0] differs significantly from solf[r]. Since the difference between the two solutions increases with zborder, I assume it is caused by some numerical inaccuracy/instability. But increasing/decreasing WorkingPrecision shows no significant influence.

Are there some ideas where this deviation comes from and how to avoid it?

• I do not have the time to answer this in detail now, but use a finer mesh like Method -> {"FiniteElement", "MeshOptions" -> {"MaxCellMeasure" -> 0.1}} or even finer. See if that helps. – user21 Sep 3 '14 at 10:15
• Thanks for the answer. Your hint solved the problem! – Thomas Sep 3 '14 at 11:17
• Can someone mark this question as solved? – Thomas Sep 3 '14 at 11:21

Clear["Global*"];
start = 1/1000;
stop = 20;
zborder = 10;
wp = 50;
pde1 = D[1/r*D[r*f[r], r], r] - f[r] + 1/(2*Pi*r) == 0;
pde2 = D[1/r*D[r*g[r, z], r], r] - g[r, z] + 1/(2*Pi*r) == 0;
bc1 = {f[start] == 0, f[stop] == 1/(2*Pi*stop)};
bc2 = {g[start, z] == solf[start], g[stop, z] == solf[stop],
g[r, zborder] == g[r, -zborder] == solf[r]};
solf = NDSolveValue[{pde1, bc1}, f, {r, start, stop},
WorkingPrecision -> wp];
solg = NDSolveValue[{pde2, bc2},
g, {r, start, stop}, {z, -zborder, zborder}];

NDSolveValue::femcscd: The PDE is convection dominated and the result may not be stable. Adding artificial diffusion may help.


Since the issus is at the boundary what you can do is to specify a mesh with a fine boundary:

Needs["NDSolveFEM"]
mesh = ToElementMesh[Rectangle[{start, -zborder}, {stop, zborder}],
"MaxBoundaryCellMeasure" -> 1/2/2/2/2,
MeshElementType -> TriangleElement]
mesh["Wireframe"] This will improve the solution:

solg = NDSolveValue[{pde2, bc2}, g, {r, z} \[Element] mesh]; // Quiet
Plot[{solg[r, 0] - solf[r]}, {r, start, stop}, PlotRange -> All]
` If you want to further improve the solution you could just refine the left hand side of the domain by constructing a boundary mesh manually. Also, there is a similar example Stabilization of Convection-Dominated Equations