1
$\begingroup$

I want to solve the Laplace problem $-\Delta u=f$ ("analyst's Laplacian") for a given $f$ and unknown $u$ on the torus $[0,2\pi]/\sim$, i.e. a rectangle with opposite sides identified. The solution is unique if one imposes $u(0,0)=0$. As a test I chose $f(x,y)=\sin(x)$. I modified code from the documentation and now have the following code:

ufun=NDSolveValue[{-Laplacian[u[x,y],{x,y}]==Sin[ x],DirichletCondition[u[x,y]==0,x==0&&y==0],PeriodicBoundaryCondition[u[x,y],0<=y<=2 Pi&&x==2 Pi,FindGeometricTransform[{{0,0},{0,2 Pi}},{{2 Pi,0},{2 Pi,2 Pi}}][[2]]],PeriodicBoundaryCondition[u[x,y],0<=x<=2 Pi&&y==2 Pi,FindGeometricTransform[{{0,0},{2 Pi,0}},{{0,2 Pi},{2 Pi,2 Pi}}][[2]]]},u,{x,y}\[Element]Rectangle[{0,0},{2Pi,2 Pi}]];
Plot3D[ufun[x,y],{x,y}\[Element]Rectangle[{0,0},{2 Pi,2 Pi}]]

The result is the following plot:

Resulting plot

This is unexpected for me, because the solution to the equation $-\Delta u=\sin(x)$ with $u(0,0)=0$ is $u(x,y)=\sin(x)$. That is different from the plot. What went wrong?

$\endgroup$
4
  • 2
    $\begingroup$ The general solution is your sin(x) plus anything that solves Laplacian[u[x,y],{x,y}]==0 (any Harmonic function) $\endgroup$ Dec 28, 2023 at 20:37
  • $\begingroup$ @CraigCarter indeed, but how does mathematica pick this specific solution? $\endgroup$
    – chris
    Dec 28, 2023 at 22:33
  • 1
    $\begingroup$ @chris Evaluate ufun["ElementMesh"], we see that problem solved with FEM on the mesh of 400 quad element. The explanation for why the erroneous result was obtained apparently lies in the fact that the implementation of periodically boundary conditions is incorrect in this case. $\endgroup$ Dec 29, 2023 at 3:58
  • 1
    $\begingroup$ @CraigCarter The only harmonic function on the torus with $u(0,0)=0$ is the constant zero function, so this should shouldn't add any ambiguity. $\endgroup$
    – user505117
    Dec 29, 2023 at 16:03

2 Answers 2

3
$\begingroup$

New answer (06 Jan 2024)

My old answer was not clean about the way the unicity of the solution is imposed with a Dirichlet condition.

Here is a clean approach :

It imposes the solution to be 0 at a arbitrary point called dirichletConditionPoint which is randomly taken to be at position {4,4}. The point is the small black ball in the graphic. One can impose a value different from 0, but in this case the value would be different from the specified value (this is normal)

Get["NDSolve`FEM`"]

dirichletConditionPoint = {4, 4};

offset = 1;

myBoundaryMesh = 
  ToBoundaryMesh[Rectangle[{-offset, -offset}, {2 Pi, 2 Pi}], 
   "IncludePoints" -> {dirichletConditionPoint }];
myMesh = ToElementMesh[myBoundaryMesh];

ufun = NDSolveValue[{-Laplacian[u[x, y], {x, y}] == Sin[x]
    , DirichletCondition[u[x, y] == 0, 
     Norm[{x, y} - dirichletConditionPoint] < 10^-3]
    , PeriodicBoundaryCondition[u[x, y], x == 2 Pi, 
     TranslationTransform[{-2 Pi, 0}]]
    , PeriodicBoundaryCondition[u[x, y], y == 2 Pi, 
     TranslationTransform[{0, -2 Pi}]]
    , PeriodicBoundaryCondition[u[x, y], x == -offset, 
     TranslationTransform[{2 Pi, 0}]]
    , PeriodicBoundaryCondition[u[x, y], y == -offset, 
     TranslationTransform[{0, 2 Pi}]]}
   , u
   , {x, y} \[Element] myMesh];

Show[
 Plot3D[ufun[x, y], {x, y} \[Element] myMesh]
 , Graphics3D[{Black, 
   Ball[Append[dirichletConditionPoint, 
     Apply[ufun, dirichletConditionPoint]], .2]}]]  

enter image description here

Old answer

Here is the very triky approach I explain in this answer :

https://mathematica.stackexchange.com/a/174514/5467

There is one more difficulty : Without any Dirichlet condition the solution is unique up to a constant.
With a Dirichlet condition on a small part of the boundary, it works :

offset = 1;

ufun = NDSolveValue[{-Laplacian[u[x, y], {x, y}] == Sin[x], 
    DirichletCondition[
     u[x, y] == -offset, -offset < x < -offset + .2 && y == -offset], 
    PeriodicBoundaryCondition[u[x, y], x == 2 Pi, 
     TranslationTransform[{-2 Pi, 0}]
    (* TranslationTransform[{-2 Pi, 0}] is equivalent to OP's code but is cleaner *)],
    PeriodicBoundaryCondition[u[x, y], y == 2 Pi, 
     TranslationTransform[{0, -2 Pi}]],
    PeriodicBoundaryCondition[u[x, y], x == -offset, 
     TranslationTransform[{2 Pi, 0}]], 
    PeriodicBoundaryCondition[
     u[x, y], (x > -offset + .3) && y == -offset, 
     TranslationTransform[{0, 2 Pi}]]
    }, u, {x, y} \[Element] 
    Rectangle[{-offset, -offset}, {2 Pi, 2 Pi}]];

Plot3D[ufun[x, y], {x, y} \[Element] Rectangle[{0, 0}, {2 Pi, 2 Pi}]]

enter image description here

Notice that the restriction of my extended domain [-1,2 Pi]X[-1,2 Pi] to your domain [0,2 Pi]X[0,2 Pi] is exactly periodic in your domain.

Disclaimer : I'm not a specialist. There may be problems due to the topology of the domain. Nevertheless I think that the code is worth to be shared.

$\endgroup$
4
  • $\begingroup$ The answer works, but I don't understand why. Why did you choose .2 in the line "-offset < x < -offset + .2"? My code also had a Dirichlet condition that made the solution to the equation unique. $\endgroup$
    – user505117
    Jan 3 at 19:14
  • $\begingroup$ A Dirichlet condition can't be at the same location than a periodic condition (if you try, you get a explicit error message), but it can be anywhere else. Here I have arbitrarily put it on a small segment of the boundary where I have removed the periodic constraint. It would be cleaner to put it somewhere inside the domain. It's probably feasable but would involve more modifications of your code (I have tried quickly but didn't have enough time to do it) $\endgroup$
    – andre314
    Jan 3 at 19:35
  • $\begingroup$ Internal Dirichlet condition are possible.Here is an example $\endgroup$
    – andre314
    Jan 3 at 19:49
  • $\begingroup$ Answer updated : The Dirichlet condition is now on one point somewhere inside the domain. $\endgroup$
    – andre314
    Jan 6 at 17:37
-2
$\begingroup$

The inhomogenity $\sin x$ is a local sink and source term.

The solution yields a current on the torus between the positive and negative areas of $\sin$.

To clamp down the potential at a point in the corners, is somehow an artefact, that generates infinite current densities.

Do VectorPlot of the Gradient without and with Dirichlet.

From physical reasons, a contact with constant potential has to have a finite contact lenght in order to generate a regular solution. Experiment with a line segment and a circle of constant potential in the center.

$\endgroup$
3
  • 2
    $\begingroup$ I don't understand this answer because I don't know some of the words (local sink, local source, infinite current density) and I don't the facts you state (e.g. solution yields a current). In your instruction "Do VectorPlot of the Gradient" I don't know which gradient you mean and I don't know what you mean by Dirichlet. In your instruction "Experiment with a line segment and a circle" I don't know what kind of experiment you're suggesting. It seems I'm lacking a lot of knowledge to understand this answer and there may be no way to make it understandable for me. $\endgroup$
    – user505117
    Dec 29, 2023 at 16:26
  • $\begingroup$ A note about "... that generates infinite current densities." The current is discretized, not anormaly infinite. I would be infinite in a continuous domain, but the principle of F.E.M. is to discretize the domain, and after discretization, it's perfectly right to impose a Dirichlet condition at only one point. To go back to the continuous domain, consider that the discretized current is the integral of the distributed current over some continuous small domain (roughly a few elements of the mesh around the point). $\endgroup$
    – andre314
    Dec 29, 2023 at 20:22
  • $\begingroup$ @Roland F I think your text originates from some considerations about this Question. As such it is not understable in the present question. I have created yesterday this chatroom where one can freely discuss without restrictions. $\endgroup$
    – andre314
    Dec 29, 2023 at 20:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.