3

How can I obtain the 3D numerical solution for the potential (or field) due to a point charge inside a cubic domain with periodic boundary conditions in all directions? I guess I can use NDSolveValue with first argument:

-Laplacian[u[x, y, z], {x, y, z}] == diracD

to represent Gauss' Law $\nabla^2 \phi(r)=-q\delta(r)/\epsilon_0$. The function diracD approximates a Dirac delta function, e.g.: $diracD(x,y,z,\epsilon)=(\epsilon/ \pi)/(x^2 + y^2+z^2 + \epsilon^2)$ with some small value for $\epsilon$. I probably need to also make use of PeriodicBoundaryCondition.

kotozna
  • 319
  • 1
  • 6
  • You can't really use DiracDelta in numerical computation. It works only with DSolve. This is an analytical entity and has no meaning for numerical work. Math people call it a distribution. You need to approximate it with a real function that NDSolve can understand. Like a very large spike at some location with a very area it sits on or something like this. – Nasser Oct 15 '23 at 20:31
  • @Nasser Thanks, I replaced it with a numerical approximation to a Dirac delta. – kotozna Oct 15 '23 at 20:55
  • For this simple case, you can solve it symbolically with DSolve[Laplacian[u[r], {r, th, phi}] == k DiracDelta[r], u, {r}], then plug numbers in to the solution: {{u -> Function[{r}, C[1] + r C[2] + k r HeavisideTheta[r]]}}. Is there a reason you want to do this numerically? –  Oct 16 '23 at 01:16
  • @jdp Thanks but I want to solve it numerically because I will use this as a starting point for more complicated cases. The solution also needs to be periodic. – kotozna Oct 16 '23 at 02:12
  • 1
    Okay. I found an example in the documentation which is close to what you want, except for the delta function. ref/PeriodicBoundaryCondition under Scope/3D Problems shows how to set up a Laplacian in a cube, using Cartesian coordinates. For the delta function, try n*HeavisidePi[x/n], where n is some suitable constant. The x/n determines the width, and multiplying by n fixes the normalization,which should be unity. –  Oct 16 '23 at 02:53
  • 1D analogue of this problem was discussed here https://mathematica.stackexchange.com/questions/209169/solve-a-one-dimensional-heat-transfer-problem-with-ndsolve/210124#210124 May be it can help you – Oleksii Semenov Oct 16 '23 at 07:47
  • I'm missing the definition of the periodic domain – Ulrich Neumann Oct 16 '23 at 08:14

1 Answers1

4

Based on @jdp's helpful hint, try:

eps = .01
dirac = Function[x, Exp[-(x^2/(2 eps))]/Sqrt[2 Pi eps]]

cube=Cuboid[{0, 0, 0}, {1, 1, 1}]

Examplary solution for DiracDelta at x=y=z=1/2

U  = NDSolveValue[{-Laplacian[u[x, y, z], {x, y, z}] == 
dirac[x - 1/2] dirac[y - 1/2] dirac[z - 1/2], 
PeriodicBoundaryCondition[u[x, y, z], x == 1, 
TranslationTransform[{-1, 0, 0}]],
PeriodicBoundaryCondition[u[x, y, z], x == 0, 
TranslationTransform[{ 1, 0, 0}]],
PeriodicBoundaryCondition[u[x, y, z], y == 1, 
TranslationTransform[{0, -1 , 0}]],
PeriodicBoundaryCondition[u[x, y, z], y == 0, 
TranslationTransform[{ 0, 1, 0}]],
PeriodicBoundaryCondition[u[x, y, z], z == 1, 
TranslationTransform[{0, 0, -1}]],
PeriodicBoundaryCondition[u[x, y, z], z == 0, 
TranslationTransform[{ 0, 0, 1}]]}
, u, {x, y, z} \[Element] cube]

SliceContourPlot3D[U[x, y, z], 
{{"XStackedPlanes", {0.25, 0.5, .75}}, {"YStackedPlanes",1}, {"ZStackedPlanes", 1}}, {x, y, z} [Element] cube, ColorFunction -> "TemperatureMap", Boxed -> False, Axes -> False]

enter image description here

Hope it helps!

Ulrich Neumann
  • 53,729
  • 2
  • 23
  • 55