Plotting heat equation in a Manipulate expression [duplicate]

Question

I have a 2D heat equation $u_t = \alpha (u_{xx} + u_{yy})$ with conditions:

• $u(x, y, 0) = 300$,
• $u_y(x, 0, t) = \mu_1(x)$,
• $u_y(x, 1, t) = \mu_2(x)$,
• $u(0, y, t) = \mu_3(y)$,
• $u(1, y, t) = \mu_4(y)$,

where:

• $\mu_1(x) = c_1$ if $|x - 0.5| < 0.25,$ else $0$,
• $\mu_2(x) = c_2$ if $|x - 0.5| < 0.25,$ else $0$.
• $\mu_3(y) = \mu_4(y) = T_0$

I am trying to solve it with the code, given below:

c1 = 100.0; c2 = -500.0; T0 = 500;

a = 0.002;

phi[x_, y_] = 300;

mu1[x_] := If[Abs[x - 0.5] <= 0.25, c1, 0];
mu2[x_] := If[Abs[x - 0.5] <= 0.25, c2, 0];

mu3[y_] := T0;
mu4[y_] := T0;

Off[NDSolve::ibcinc];

Manipulate[
  Module[{plt, sol},
    sol = NDSolve[
      { 
        D[u[x, y, t], t] == a (D[u[x, y, t], x, x] + D[u[x, y, t], y, y]),
        u[x, y, 0] == phi[x, y],
        (D[u[x, y, t], y] /. y -> 0) == mu1[x],
        (D[u[x, y, t], y] /. y -> 1) == mu2[x],
        u[0, y, t] == mu3[y],
        u[1, y, t] == mu4[y]
      },
      u[x, y, t], {x, 0, 1}, {y, 0, 1}, {t, 0, 100}];
   plt = ContourPlot[
      u[x, y, t] /. sol /. t -> time, {x, 0, 1}, {y, 0, 1}, 
      Contours -> 15, 
      ColorFunction -> "TemperatureMap", 
      PlotLabel -> "Temperature profile", 
      ContourLabels -> True, 
      ImageSize -> Large]
  ],
  {{time, 5, "Time (t)"}, 5, 100, 0.5}
]

The problem is that I set $c_1, c_2$ to different values (even functions of $x$) but I can't see any difference on the plot. I think I should see some effect. If my code is incorrect, how do I fix it?

Answer 1

Update

Recently I revisited the obscure tutorial for "MethodOfLines" and figured out how to fix your problem i.e. make NDSolve use (D[u[x, y, t], y] /. y -> 0) == mu1[x] and (D[u[x, y, t], y] /. y -> 1) == mu2[x] for solving the equation rather than ignore them. What you need to do is:

(1) Add Method -> {"MethodOfLines", "DifferentiateBoundaryConditions" ->{True, "ScaleFactor" -> 1}} to NDSolve.

This manual option adjustment is necessary because NDSolve will automatically set "ScaleFactor" -> 0 for non-Dirichlet boundary condition, when i.c. and b.c. are inconsistent, this will make the Neumann inconsistent b.c. that's not a function of t be completely ignored.

(2) Use UnitStep instead of Piecewise or If to define mu1 and mu2.

It's probably a bug, luckily we have Simplify`PWToUnitStep and PiecewiseExpand so don't need to reconstruct the function by hand.

The following is the fixed code:

c1 = -200.0; c2 = -200.0;
T0 = 500;
a = 0.002;
phi[x_, y_] = 300;

mu1[x_] = Simplify`PWToUnitStep@PiecewiseExpand@If[Abs[x - 0.5] <= 0.25, c1, 0]
mu2[x_] = Simplify`PWToUnitStep@PiecewiseExpand@If[Abs[x - 0.5] <= 0.25, c2, 0]

mu3[y_] := T0;
mu4[y_] := T0;

sol = u /. First@
   NDSolve[{D[u[x, y, t], t] == a (D[u[x, y, t], x, x] + D[u[x, y, t], y, y]), 
     u[x, y, 0] == phi[x, y], (D[u[x, y, t], y] /. y -> 0) == 
      mu1[x], (D[u[x, y, t], y] /. y -> 1) == mu2[x], u[0, y, t] == mu3[y], 
     u[1, y, t] == mu4[y]}, u, {x, 0, 1}, {y, 0, 1}, {t, 0, 100}, 
    Method -> {"MethodOfLines", 
      "DifferentiateBoundaryConditions" ->{True, "ScaleFactor" -> 1}}]

Manipulate[ContourPlot[sol[x, y, time], {x, 0, 1}, {y, 0, 1}, Contours -> 15, 
  ColorFunction -> "TemperatureMap", PlotLabel -> "Temperature profile", 
  ContourLabels -> True, ImageSize -> Large]
 , {{time, 5, "Time (t)"}, 5, 100, 0.5}
 ]

Now mu1 and mu2 apparently plays roles in the solution.

---

Original Answer

Well, this answer is quite incomplete because I can't fix the problem. I post this answer just to point out that the true reason for the changeless sol is the BCs (D[u[x, y, t], y] /. y -> 0) == mu1[x] and (D[u[x, y, t], y] /. y -> 1) == mu2[x] are largely ignored by NDSolve. If you try:

Clear["`*"]
T0 = 500; a = 0.002; phi[x_, y_] = 300;

(*   I changed your definitions for mu1[x] and mu2[x] a little 
   since it doesn't work at least for version 8. *)
mu1[x_] := Piecewise[{{c1, Abs[x - 0.5] <= 0.25}}]; 
mu2[x_] := Piecewise[{{c2, Abs[x - 0.5] <= 0.25}}]; 
mu3[y_] := T0; 
mu4[y_] := T0;

sol[c1_, c2_] := NDSolve[{D[u[x, y, t], t] == a (D[u[x, y, t], x, x] + D[u[x, y, t], y, y]), 
                          u[x, y, 0] == phi[x, y], 
                          (D[u[x, y, t], y] /. y -> 0) == mu1[x], 
                          (D[u[x, y, t], y] /. y -> 1) == mu2[x], 
                          u[0, y, t] == mu3[y], 
                          u[1, y, t] == mu4[y]}, 
                          u[x, y, t], {x, 0, 1}, {y, 0, 1}, {t, 0, 100}];

Plot3D[Evaluate[D[u[x, y, t] /. sol[100, -500], y] /. y -> 1], 
       {x, 0, 1}, {t, 0, 100}, PlotRange -> All] 
Plot3D[Evaluate[D[u[x, y, t] /. sol[100, -500], y] /. y -> 0], 
       {x, 0, 1}, {t, 0, 100}, PlotRange -> All]

you'll see:

Apparently mu1[x] and mu2[x] are almostly (maybe essentially…) not used. You can try several other c1 and c2 or even take off the BCs for $y=0$ and $y=1$ to confirm.

I guess some changes for the Method will help but haven't succeeded until now… Look forward to someone who's more experienced in solving PDEs！

Answer 2

I think that one of the main issues related to the perceived responsiveness of Manipulate is still missing and should be addressed in this thread. Same as SynchronousUpdating, TrackedSymbols is an option which must be taken into account when using Manipulate.
In general, if you have a toy application for Manipulate you won't need any of the two options, however, the more complicated your application the more likely you will need to use these two options.
TrackedSymbols is to tell Manipulate when to trigger an update. If you don't specify which symbols to track Manipulate may update much more often than what you would expect. In general you may think of Manipulate as a Mathematica expression that is wrapped with Dynamic and that will thus reevaluate its contents whenever any symbols in the expression changes. To avoid continuous and unexpected reevaluations use TrackedSymbols to specify which symbols shuld trigger an update and your Manipulate will most probably start behaving as you expect.

One more comment: the more complicated your application is the more probably Manipulate will start becoming a bad tool to address your problem. In this case you should start considering using DynamicModule which is basically the underlying backbone of Manipulate. I have worked extensively with Dynamic tools and widgets (this screencasts shows a good example) and I can tell that starting to use DynamicModule requires a fairly large amount of practice and study time.

Answer 3

There are two things about your Manipulate that stand out as possible sources of trouble.

1. The use of Module to localize variables within a Manipulate isn't a good idea. It is better to introduce an invisible control to create a local variable. Also note that the variable plt isn't needed at all.

2. For complex computations such as you are carrying out, it is a good idea to use the SynchronousUpdating -> False option. I suspect that not having this option set to False is the main source of your problem.

Making the required changes to your code to correct these two problems produces a Manipulate in which the plot the dynamically changes when the time slider is moved.

Manipulate[
  sol = NDSolve[{
    D[u[x, y, t], t] == a (D[u[x, y, t], x, x] + D[u[x, y, t], y, y]), 
    u[x, y, 0] == phi[x, y], 
    (D[u[x, y, t], y] /. y -> 0) == mu1[x], 
    (D[u[x, y, t], y] /. y -> 1) == mu2[x], 
    u[0, y, t] == mu3[y], 
    u[1, y, t] == mu4[y]}, 
    u[x, y, t], 
    {x, 0, 1}, {y, 0, 1}, {t, 0, 100}];
  ContourPlot[u[x, y, t] /. sol /. t -> time, {x, 0, 1}, {y, 0, 1},
    Contours -> 15,
    ColorFunction -> "TemperatureMap", 
    PlotLabel -> "Temperature profile",
    ContourLabels -> True],
  {{time, 5, "Time (t)"}, 5, 100, 0.5},
  {sol, ControlType -> None},
  SynchronousUpdating -> False]

Here is the plot a t =5.

Here is the plot a t =30.