How to get the answer I expect with NDSolve

Question

I have been struggling with FEM method in NDSolve. I get the test problem "working" however it does not match well with the exact solution. My original problem has angle dependence as well. I was wondering how I can improve this result. There are couple of things that I tried:

1-Decreasing the Max cell measure from .01 to .001. It literally only changed things on the 4th digit.

2-I tried to change the interpolation order however, I got an error saying that it can not be less than the mesh order.

3- Why do I always have to install the FEM package?

4- I know you do not need to impose vanishing Neumann boundary conditions. However when we have more than one unknown functions why does not that automatically assume that we have the vanishing Neumann for that unknown function even if we did not add a zero Neumann value for the corresponding equation to that unknown. (Here I am talking about coupled elliptic PDE's.)

Would this complaint cause poor results?

    Needs["NDSolve`FEM`"]
    eq111 = D[D[\[Psi][r, \[Theta]], r], r] + 
    Sin[\[Theta]]/r^2*
    D[1/Sin[\[Theta]]*D[\[Psi][r, \[Theta]], \[Theta]], \[Theta]];
    a = 1;
    b = 2;
    \[CapitalOmega] = Rectangle[{a, 0}, {b, Pi}];
    RegionPlot[\[CapitalOmega], AspectRatio -> .5]
    mesh = ToElementMesh[\[CapitalOmega], 
    MaxCellMeasure -> {"Length" -> .01}];
    ui = NDSolveValue[
   {eq111 == 
   0, {DirichletCondition[\[Psi][r, \[Theta]] == -1/(4*Pi*a), 
   r == a && 0 < \[Theta] <= Pi], 
   DirichletCondition[\[Psi][r, \[Theta]] == -1/(4*Pi*b), 
   r == b && 0 < \[Theta] <= Pi]}},
   \[Psi], Element[{r, \[Theta]}, mesh], 
   Method -> {"FiniteElement", "InterpolationOrder" -> {\[Psi] -> 2}}
   ]
   ContourPlot[ui[r, \[Theta]], {r, a, b}, {\[Theta], 0, Pi}, 
   AspectRatio -> 0.5, ColorFunction -> "TemperatureMap", 
   PlotLegends -> Automatic] 

This looks like it is working in the sense that the solution should not depend on the angle and it does not. However we know the exact solution is -1/4*Pi*r.

    ui[1.5, Pi/8]

    -1/(4*Pi*1.5) // N

They are quite different from each other.

Another question that I have is would transforming the spherical mesh result into cartesian result mess up the final answer? I appreciate the time and the suggestions.

Answer 1

NDSolveValue produces the correct solution but you provide wrong boundary conditions for what you want to do. Here is the corrected and simplified solution:

eq111 = D[D[p[r, t], r], r] + 
   Sin[t]/r^2*D[1/Sin[t]*D[p[r, t], t], t];
a = 1; b = 2;
reg = Rectangle[{a, 0}, {b, Pi}];
asol[r_, t_] := -1/4*Pi*r
ui = NDSolveValue[{eq111 == 0, 
    DirichletCondition[p[r, t] == asol[r, t], r == a || r == b]
    }, p, Element[{r, t}, reg]];
Plot3D[asol[r, t] - ui[r, t], Element[{r, t}, reg]]

This is a result with 10^-13 difference to your analytical function.

Now to the questions:

1) Because NDSolveValue did compute the solution to the problem you set up. It could not do much better by adding more elements.

2) Please read the documentation. The mesh order is changed via the MeshOrder option of ToElementMesh. The change of InterpoaltionOrder is a very specialized functionality that is documented in the section Ordering of Dependent Variable Names - all of this is not needed in your case.

3) You do not need that - as shown in the answer

4) Please ask a new question with an example of your claim.

If I used

asol[r_, t_] := -1/(4*Pi*r)

I get:

eq111 /. p -> asol // InputForm
-1/(2*Pi*r^3)