Search code examples
wolfram-mathematicadifferential-equations

How to specify initial conditions at parts of domain in Mathematica?

I try to compare solution of my Finite Volume Method to solution of Mathematica for simple Heat Equation U_t = U_xx. For this I have to specify initial and boundary condition for NDSolve function in Mathematica. I would like to have on boundaries U = 90. As initial condition I would like to have 100 in all square domain except boundaries. How can I do this? Obviously my code doesn't work. 

FUNC = NDSolve[{D[T[x, y, t], t] == (D[T[x, y, t], x, x] + D[T[x, y, t], y, y]), 
    T[x, y, 0] == 100, T[0, y, t] == 90, T[9, y, t] == 90, 
    T[x, 0, t] == 90, T[x, 9, t] == 90}, 
   T, {x, 0, 9}, {y, 0, 9}, {t, 0, 6}];

It try to set initial condition at all domain equals to 100.

Solution

  • Actually, with a small modification your code does work.

    FUNC = T /. 
  NDSolve[{
     D[T[x, y, t], t] == (D[T[x, y, t], x, x] + D[T[x, y, t], y, y]), 
     T[x, y, 0] == 100, T[0, y, t] == 90, T[9, y, t] == 90, 
     T[x, 0, t] == 90, T[x, 9, t] == 90}, 
    T, {x, 0, 9}, {y, 0, 9}, {t, 0, 10}][[1]]

    Modifications are the T/. (ReplaceAll) part in front and the [[1]] (Part) at the end; you may want to look these op in the documentation. Thery are necessary to massage the output into the correct shape .

    You get an error message about inconsistency of boundary and initial conditions (which is correct). The result is nevertheless usable.

    Alternatively, you can change the initial condition so as to read

    FUNC = T /. 
  NDSolve[{
     D[T[x, y, t], t] == (D[T[x, y, t], x, x] + D[T[x, y, t], y, y]), 
     T[x, y, 0] == If[x == 0 || x == 9 || y == 0 || y == 9, 90, 100], 
     T[0, y, t] == 90, T[9, y, t] == 90, T[x, 0, t] == 90, 
     T[x, 9, t] == 90}, T, {x, 0, 9}, {y, 0, 9}, {t, 0, 10}][[1]]

    You can now do lots of things with this function like e.g. a Manipulate:

    Manipulate[ContourPlot[FUNC[x, y, t], {x, 0, 9}, {y, 0, 9}], {t, 0, 10}]

    Manipulate ContourPlot

    or an animated GIF:

    anim = Table[DensityPlot[FUNC[x, y, t], {x, 0, 9}, {y, 0, 9}, 
          ColorFunctionScaling -> False, PlotPoints -> 50,
          ColorFunction -> (Hue[Rescale[#, {50, 100}], 1, 1] &)], 
        {t, 0, 10, .2}];
Export[ToFileName[$UserDocumentsDirectory, "anim.gif"], anim, "GIF"]

    animated gif