18

I am new to Mathematica, a friend recommended this software and started using it, in fact download the trial version to know.

I recently did a program in C to calculate numerically the solution to the Laplace equation in two dimensions for a set of points as in the figure.

enter image description here

The result was very good, finding the image below.

enter image description here

This program took me about 100 lines in C, my friend told me that Mathematica could do it in a couple of lines, which seemed quite interesting. I studied a bit and found that Mathematica can solve the Laplace and Poisson equations using NDSolve command.

However, this command requires to be given to the specific boundary conditions. The boundary condition in which $\phi = 0$, it is quite easy to introduce. But on the inside border, where $\phi = 100$, I failed to get the condition. My idea is to make this problem as follows:

            uval = NDSolveValue[{D[u[x, y], x, x] + D[u[x, y], y, y] == 
            0, 
            u[x, 0] == u[x, 100] == u[0, y] == u[100, y] == 0, 
            u[40, y] == u[60, y] == u[x, 40] == u[x, 60] == 100}, 
            u, {x, 0, 100}, {y, 0, 100}]

But this does not work, and I could not give the boundary conditions on the inner square. Any help would be the most grateful.

Santi Carmesí
  • 183
  • 1
  • 5
  • 1
    Aren't there conflicting boundary conditions? For instance, u[x,40] equals u[100,y] at {x,y}=={100,40} making you state 0==100 at that point. – Sjoerd C. de Vries Sep 11 '14 at 19:03

2 Answers2

26

You need a DirichletCondition (new in V10) here. Using regions (also from V10):

Ω = RegionDifference[Rectangle[{0, 0}, {100, 100}], Rectangle[{40, 40}, {60, 60}]];
sol = NDSolveValue[{
   D[u[x, y], x, x] + D[u[x, y], y, y] == 0,
   DirichletCondition[u[x, y] == 100., 
     x == 40 && 40 <= y <= 60 || 
     x == 60 && 40 <= y <= 60 || 
     40 <= x <= 60 && y == 40 || 
     40 <= x <= 60 && y == 60],
   u[x, 0] == u[x, 100] == u[0, y] == u[100, y] == 0
   }, u, {x, y} ∈ Ω]

DensityPlot[sol[x, y], {x, y} ∈ Ω,  Mesh -> None, ColorFunction -> "Rainbow", 
            PlotRange -> All, PlotLegends->Automatic]

Mathematica graphics

Sjoerd C. de Vries
  • 65,815
  • 14
  • 188
  • 323
  • Thank you very much. I wondered, with a couple of lines !!! Just one more question friend, if not a lot to ask. How could plot the vector field phi? That is, graphic $-\nabla sol(x,y)$. Mi idea es hacer: 
          field = -Grad[sol[x, y], {x, y}]
      VectorPlot[campo[x, y], {x, 0, 100}, {y, 0, 100}]

    But I shows a blank graph.

     – Santi Carmesí Sep 11 '14 at 22:17
  • And I get the messages:

    InterpolatingFunction::dmval: Input value {42.8643,42.8643} lies outside the range of data in the interpolating function. Extrapolation will be used. >>

    General::stop: Further output of InterpolatingFunction::dmval will be suppressed during this calculation. >>

     – Santi Carmesí Sep 11 '14 at 22:19
  • What's going on? I'll be even more grateful. – Santi Carmesí Sep 11 '14 at 22:20
  • 1
    @SantiCarmesí For the gradient and field strength, you can use my answer here. The error message in the potential plot appears because the inner boundary isn't sampled quite accurately. One can get rid of this by first defining rm=RegionMember[\[CapitalOmega]] and then replacing the plot argument sol[x,y] by If[rm[{x,y}],sol[x,y],100]. – Jens Sep 12 '14 at 03:33
  • @SantiCarmesí The gradient field issue is discussed in great detail here. – Mark McClure Sep 12 '14 at 09:15
  • I was able to figure out my particular geometry, but how did you get this plot to look so smooth? – Young Jun 25 '16 at 19:48
  • but why when you define the boundary condition it is || rathet than && – dcydhb May 21 '18 at 14:02
  • @dcydhb since you’re describing a general point on the box it can be on each of the four sides. Using && would require a point that is on all four sides at once, which is impossible. – Sjoerd C. de Vries May 21 '18 at 14:08
  • @SjoerdC.deVries but when you use || it means or,so it means we define one side or two side or three side or four side is 100 not all of them is 100 – dcydhb May 21 '18 at 14:50
  • @SjoerdC.deVries i mean how you distinct you defined one side or you defined all of them is 100 – dcydhb May 21 '18 at 14:53
  • @dcydhb it’s u that should equal 100. That’s what is stated in the first line of the condition. The next 4 lines define where this condition should hold, that is, anywhere on the box. So, that’s on the top of the box OR on its left side OR on it’s right side OR on its bottom. – Sjoerd C. de Vries May 21 '18 at 14:56
  • @SjoerdC.deVries so you means it is only one side is 100 rather than four of them?if so,the result should show different of one side from other side,but it diesn't show and by the way,how to define four of thrm are ALL 100 – dcydhb May 21 '18 at 15:08
  • @SjoerdC.deVries also in the question,all four side are 100 and only one side is 100 should be different and how to distinguish them – dcydhb May 21 '18 at 15:29
  • @dcydhb it looks like you have some difficulties with the logical OR. It’s not an exclusive or. The condition states that a general boundary point may be on any of the square’s sides. Any point on any side is OK. But you cannot demand that a general point is on all sides at the same time, which is what an AND would do. – Sjoerd C. de Vries May 21 '18 at 16:08
  • @SjoerdC.deVries yes but in the question it doesn't mean all four side of all points should be 100? and if we define all of them are 100 and if we drfine only one side is 100 and other three side is 0,they are two kind of situation,how to distinguish,ut is what i want to know – dcydhb May 22 '18 at 01:36
  • @SjoerdC.deVries and also i have tried that define four side respectively and makes it are all 100,and the result is the same as your ||, why || not && – dcydhb May 22 '18 at 07:09
17

Building on @SjoerdC.deVries answer you can use:

sol = NDSolveValue[{D[u[x, y], x, x] + D[u[x, y], y, y] == 0, 
   DirichletCondition[u[x, y] == 100., 
    x == 40 && 40 <= y <= 60 || x == 60 && 40 <= y <= 60 || 
     40 <= x <= 60 && y == 40 || 40 <= x <= 60 && y == 60], 
   u[x, 0] == u[x, 100] == u[0, y] == u[100, y] == 0}, 
  u, {x, y} \[Element] \[CapitalOmega], 
  Method -> {"FiniteElement", 
    "MeshOptions" -> {"BoundaryMeshGenerator" -> "Continuation"}}]

Adding this option results in a better boundary approximation and you will not get the warning messages about extrapolation at the boundary.

user21
  • 39,710
  • 8
  • 110
  • 167