5
$\begingroup$

I would like to solve an axisymmetric Poisson equation on a disk with FEM.

The code for Cartesian coordinates works fine:

NDSolve[{Laplacian[u[x, y], {x, y}, "Cartesian"] == x^2 + y^2, DirichletCondition[u[x, y] == 0, True]}, u, {x, y} \[Element] Disk[], Method -> "FiniteElement"]

But going axisymmetric (which should make equation easier) leads to error:

NDSolve[{Laplacian[u[r], {r, phi}, "Polar"] == r^2, DirichletCondition[u[r] == 0, r == 1]}, u, {r, 0, 1}, Method -> "FiniteElement"]
(* NDSolve::femcnmd: "The PDE coefficient {{1/r}} does not evaluate to a numeric matrix of dimensions {1,1}. " *)

It seems that Mathematica doesn't like 1/r term in polar laplacian. Multiplying by r the both sides of the equation solves the problem, but not always. Moreover, I would like to keep code free of such "patches" to be flexible to change the coordinate system without rewriting code.

What the reason of the error? Do you have any suggestion about how to solve it?

$\endgroup$
7
$\begingroup$

Update:

In version 11.0 this works out of the box:

if = NDSolveValue[{Laplacian[u[r], {r, phi}, "Polar"] == r^2, 
    DirichletCondition[u[r] == 0, r == 1]}, u, {r, 0, 1}, 
   Method -> "FiniteElement"];
Plot[if[r], {r, 0, 1}]

enter image description here

Old Answer:

If you evaluate:

Laplacian[u[r],{r,phi},"Polar"]/.r->0

Power::infy: "Infinite expression 1/0 encountered."

You can use:

if = NDSolveValue[{Laplacian[u[r], {r, phi}, "Polar"] == r^2, 
    DirichletCondition[u[r] == 0, r == 1]}, u, {r, 10^-6, 1}, 
   Method -> "FiniteElement"];
Plot[if[r], {r, 10^-6, 1}]

To get:

enter image description here

Here is a different way to do it:

Needs["NDSolve`FEM`"]
mesh = ToElementMesh[ImplicitRegion[True, {{x, 0, 1}}]];
if = NDSolveValue[{Laplacian[u[r], {r, phi}, "Polar"] == r^2, 
    DirichletCondition[u[r] == 0, r == 1]}, u, {r} \[Element] mesh, 
   Method -> {"FiniteElement"}];

Why does this work? The coefficients given in the PDE need to be test evaluated and this happens for 1D problems at the left boundary. This maybe not have been such a good choice on my side but it is correct as 0 is part of the domain. (I have an idea for an improvement in mind). In the ElementMesh case the center of mass of thee first element is used which in this case is not at zero.

I understand this is unfortunate but testing the coefficients for correctness is absolutely essential. Note, also, that you one could overwrite the test coordinate used, but for that one has to use the low level functions. Maybe I expose that at the NDSolve level as well. The jury is still out on that one.

Update

Another way to do it is to change the bound to go from -1 to 0:

if = NDSolveValue[{Laplacian[u[r], {r, phi}, "Polar"] == r^2, 
    DirichletCondition[u[r] == 0, r == -1]}, u, {r, -1, 0}, Method -> "FiniteElement"];
Plot[if[r], {r, -1, 0}]
$\endgroup$
  • $\begingroup$ As I know, FEM doesn't need function values at each point, but integral of the function over the subregions. Value at a point can be used for integration. Does anybody know which integration rule Mathematica's FEM uses and how it can be changed? $\endgroup$ – Ivan Mar 30 '15 at 17:02
  • $\begingroup$ @chris, hm, not sure I understand your photo comment. $\endgroup$ – user21 Mar 30 '15 at 17:51
  • $\begingroup$ @chris, thanks for the help. We must have miss understood each other. I was not looking for the options to NDSolve but thinking of a way to change the test coordinate for the test of the coefficient should be exposed at the level of NDSolve. I.e. in a sub option to "PDEDiscretization" ->.{"FiniteElement", "TestCoordinate"->Automatic}... $\endgroup$ – user21 Mar 30 '15 at 18:11
  • $\begingroup$ @chris, thank you! The second approach is very interesting. Do I understand correct, that NDSolve needs to chech that the coefficients are numeric and, by default, it evaluates them at the left point? By the way, for me this code also works: if = NDSolveValue[{Laplacian[u[r], {r, phi}, "Polar"] == r^2, DirichletCondition[u[r] == 0, r == 1]}, u, {r} \[Element] ImplicitRegion[0 <= r <= 1, {r}], Method -> {"FiniteElement"}]. I find this a little bit more elegant, since we change only the domain definition. $\endgroup$ – Ivan Mar 31 '15 at 17:09
  • 1
    $\begingroup$ @Ivan, that's correct. In the 1D case it uses the left coordinate of the region bound. Good that you found that the` ImplicitRegion` works as well $\endgroup$ – user21 Mar 31 '15 at 18:44

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.