2
$\begingroup$

I'v been trying to solve for the steady state of the heat diffusion equation numerically but I cant seem to get it to work. My code is

pde = D[u[x, y], {x, 2}] + D[u[x, y], {y, 2}] == 0
sol = NDSolve[{pde, u[x, 0] == 100, u[x, 10] == 400, u[0, y] == 0, 
u[10, y] == 0}, u[x, y], {x, 0, 10}, {y, 0, 10}] 

and I get this error

NDSolve::ivone: Boundary values may only be specified for one independent
variable. Initial values may only be specified at one value of the other 
independent variable.

I think It has something to do with there not being an initial value condition, but I don't have one since I'm solving for the steady state. Could anyone please tell me where the above error comes from? Any help would be appreciated. Thanks.

$\endgroup$
  • $\begingroup$ what is u[10,10]? 0 or 400? $\endgroup$ – chris Jun 24 '15 at 17:47
  • $\begingroup$ You can cook up a boundary condition without the corner discontinuity issue and you still get that same error message. $\endgroup$ – george2079 Jun 24 '15 at 17:52
  • $\begingroup$ It works by me ! Mathematica 10.0.2 on Windows 7 $\endgroup$ – andre314 Jun 24 '15 at 17:56
  • $\begingroup$ works here as well: Mathematica 10.1 on Mac OS X 10.10.4 $\endgroup$ – sunt05 Jun 24 '15 at 17:59
  • $\begingroup$ On Mathematica 8.0.4 (Windows 7), there is effectively the error $\endgroup$ – andre314 Jun 24 '15 at 18:48
1
$\begingroup$

If you have v.10 you can explicitly use the finite element method:

Needs["NDSolve`FEM`"]
mesh = ToElementMesh[Rectangle[{0, 0}, {10, 10}]]
sol = First@NDSolveValue[{Laplacian[w[x, y], {x, y}] == 0,
   DirichletCondition[w[x, y] == 100, y == 0],
   DirichletCondition[w[x, y] == 400, y == 10],
   DirichletCondition[w[x, y] == 0, x == 0],
   DirichletCondition[w[x, y] == 0, x == 10]}, {w}, {x, y} \[Element] mesh];

ContourPlot[sol[x, y], {x, 0, 10}, {y, 0, 10}]

enter image description here

$\endgroup$
  • 1
    $\begingroup$ You don't have to explicitly tell Mathematica 10 to use the finite element method. With the syntax the OP uses, it switches automatically to the finite element method. $\endgroup$ – andre314 Jun 24 '15 at 19:47
4
$\begingroup$

Since you are specifically asking about versions below 10, it may be useful to point out that this problem is equivalent to the electrostatics problem of finding the potential in a region bounded by conductors held at fixed voltages. This can be solved, e.g., with the simple relaxation method I implemented in this answer, where I actually allow for lots of other complications. The potential is solved on a rectangular grid, and the main point was to play around with a graphical method for specifying the boundary conditions. You can apply that method to the question here as follows (I'm first copying the relevant definitions from the lined answer):

defaultGridSize = 120;
step = 
  Compile[{{phi, _Real, 2}, {orig, _Real, 2}, {mask, _Real, 
     2}, {dchiX, _Real, 2}, {dchiY, _Real, 2}}, 
   Module[{f = (RotateRight[phi] + RotateLeft[phi] + 
          RotateRight[phi, {0, 1}] + RotateLeft[phi, {0, 1}])*mask + 
       orig, ex, ey}, 
    ex = (RotateRight[phi, {0, 1}] - RotateLeft[phi, {0, 1}]);
    ey = (RotateRight[phi] - RotateLeft[phi]);
    f - (dchiX*ex + dchiY*ey)*mask]];

iterate = 
  Compile[{{gridArray, _Real, 2}, {originalArray, _Real, 
     2}, {maskArray, _Real, 2}, {dchiX, _Real, 2}, {dchiY, _Real, 
     2}, {tol, _Real}}, 
   FixedPoint[step[#, gridArray, maskArray, dchiX, dchiY] &, 
    originalArray, 100000, 
    SameTest -> (Max@Abs@Flatten[#1 - #2] < tol &)]];

digitize[gr_, n_] := 
  N@ImageData@
    ColorConvert[
     Image[Show[gr, Background -> Black, 
       BaseStyle -> {Antialiasing -> False}], ImageSize -> n], 
     "GrayScale"];

createLandscape[conductors_, chargePlus_, chargeMinus_, 
   suceptibility_, nGrid_] := 
  Module[{gridConductors, gridRho, gridChi, maskList}, 
   gridConductors = digitize[conductors, nGrid];
   maskList = N[1. - Unitize[gridConductors]];
   gridRho = (digitize[chargePlus, nGrid] - 
       digitize[chargeMinus, nGrid])*maskList;
   gridChi = digitize[suceptibility, nGrid]*maskList;
   {gridConductors, gridRho, gridChi, maskList}];

poissonSolver[conductors_, chargePlus_, chargeMinus_, suceptibility_, 
   nGrid_: defaultGridSize, tolerance_: 10^(-6)] := 
  Block[{averagePotential, gridConductors, gridRho, gridChi, gridEps, 
    gridList, maskList, mask4List, dChiYList, dChiXList, 
    initialGrid}, {gridConductors, gridRho, gridChi, maskList} = 
    createLandscape[conductors, chargePlus, chargeMinus, 
     suceptibility, nGrid];
   averagePotential = Mean[Select[Flatten@gridConductors, Positive]];
   initialGrid = averagePotential*maskList + gridConductors;
   gridEps = 1. + gridChi;
   gridList = gridConductors + gridRho/(4.*gridEps);
   mask4List = maskList/4.;
   dChiYList = (RotateLeft[gridChi] - 
       RotateRight[gridChi])/(2. gridEps);
   dChiXList = (RotateLeft[gridChi, {0, 1}] - 
       RotateRight[gridChi, {0, 1}])/(2. gridEps);
   Reverse@
    iterate[gridList, initialGrid, mask4List, dChiXList, dChiYList, 
     tolerance]];

Now for the specific boundary conditions. The region extends from 0 to 10 in both directions, but I need to add some padding on the outside in order to define areas with the desired boundary potential. So the simulation area is the square extending from -1 to 11 in both directions (this is where the figure shows rectangles in various gray shades):

conductors = 
 Graphics[{{GrayLevel[.5], Rectangle[{0, 10}, {10, 11}]},
   {GrayLevel[.1], Rectangle[{10, 0}, {11, 11}]}, {GrayLevel[.2], 
    Rectangle[{-1, -1}, {11, 0}]},
   {GrayLevel[.1], Rectangle[{-1, 0}, {0, 11}]}},
  PlotRangePadding -> 0, ImagePadding -> None]

conductors

To explain how I chose the conductor potentials: in my graphical approach, the color black corresponds to the absence of any object. All other potentials are encoded in the grayscale value of the region. This means the value 0 can't be used as a potential, and therefore I shift all potentials by a value 0.1 to make them non-zero. Also, I rescaled all potentials so that they fit into the interval from .1 to 1 that can be directly translated to GrayLevel. This means I have to use the following correspondence between your boundary conditions and mine:

  • 0 $\mapsto$ .1
  • 100 $\mapsto$ .2
  • 400 $\mapsto$ .5.

There is no need to specify any of the other electrostatic quantities, except to define their corresponding graphics objects to be empty. Then I invoke the solver:

chargePlus = 
  Graphics[{}, PlotRange -> (PlotRange /. FullOptions[conductors])];

chargeMinus = 
  Graphics[{}, PlotRange -> (PlotRange /. FullOptions[conductors])];

susceptibility = 
  Graphics[{}, PlotRange -> (PlotRange /. FullOptions[conductors])];

Timing[potential = 
   poissonSolver[conductors, chargePlus, chargeMinus, susceptibility];]

(* ==> {5.65001, Null} *)

ListPlot3D[potential, PlotRange -> All, 
 PlotStyle -> {Orange, Specularity[White, 10]},
  DataRange->{{-1,11},{-1,11}}]

potential

For more details, you may want to consult the link I gave above.

$\endgroup$

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.