How to define the boundary condition in 1D Heat transfer

Question

I am trying to calculate the heat transfer among a 1-D rod, with one end insulated while the right end is immersed in constant temperature surface T=0. Assume that the initial temperature of the rod is T=1. The rod length is 5. I set up the equation like this:

$$ \frac {\partial^2 u(x,t)}{\partial x^2}-\frac {\partial u(x,t)}{\partial t} = 0\\ u(x,0)=1\\ \frac {\partial u(5,t)}{\partial x} = 0\\ u(0,t)=0 $$

sol = NDSolve[{
       eqn = D[u[t, x], t] - D[u[t, x], {x, 2}] == 0,
       u[0, x] == 1,
       u[t, 0] == 0,
       (D[u[t, x], x] /. x -> 5) == 0
       }, u, {t, 0, 50}, {x, 0, 5}]

    Plot3D[Evaluate[u[t, x] /. %], {t, 0, 50}, {x, 0, 5}, 
     PlotRange -> All]

Unfortunately, I got something like this:

NDSolve::ibcinc: Warning: boundary and initial conditions are inconsistent.

Can anyone help me with the boundary value problem?

Answer 1

To get rid of the inconsistency between your BC and IC, you can quickly ramp down your BC from 1 to 0, like so:

sol = NDSolve[{eqn = D[u[t, x], t] - D[u[t, x], {x, 2}] == 0, 
    u[0, x] == 1, 
    u[t, 0] == Exp[-1000 t], (D[u[t, x], x] /. x -> 5) == 0}, 
   u, {t, 0, 50}, {x, 0, 5}];
Plot3D[Evaluate[u[t, x] /. sol], {t, 0, 50}, {x, 0, 5}, 
 PlotRange -> All, PlotPoints -> 100, MaxRecursion -> 6]

Answer 2

From experience, I know that despite the Warning, the solution is OK.

But in any case it is more comfortable not to have a irrelevant Warning : Here, it suffices to replace u[t, 0] == 0 by u[t, 0] == If[t > 0, 0, 1]

sol = NDSolve[{
   eqn = D[u[t, x], t] - D[u[t, x], {x, 2}] == 0
   , u[0, x] == 1
   , u[t, 0] == If[t > 0, 0, 1]
   , (D[u[t, x], x] /. x -> 5) == 0}
  , u, {t, 0, 50}, {x, 0, 5}]

Plot3D[Evaluate[u[t, x] /. %], {t, 0, 50}, {x, 0, 5}, 
 PlotRange -> All]  

Answer 3

Another options is to use DSolve

ClearAll[u, x, t];
pde = D[u[x, t], t] == D[u[x, t], {x, 2}];
ic = u[x, 0] == 1;
bc = {u[0, t] == 0, Derivative[1, 0][u][5, t] == 0};
sol[0] = DSolve[{pde, ic, bc}, u[x, t], {x, t}];
sol[1] = sol[0] /. K[1] -> n;

$$ u(x,t)\to \frac{2}{5} \underset{n=1}{\overset{\infty }{\sum }}-\frac{10 e^{-\frac{1}{100} (2 n-1)^2 \pi ^2 t} \sin \left(\frac{1}{10} (2 n-1) \pi x\right)}{\pi -2 n \pi } $$

 sol[2] = Activate[sol[1] /. Infinity -> 300];
 Plot3D[Evaluate[u[x, t] /. sol[2]], {t, 0, 50}, {x, 0, 5}, PlotRange -> All]

Manipulate[
 Quiet@Plot[Evaluate[u[x, t] /. sol[2] /. t -> t0], {x, 0, 5}, 
   PlotRange -> {Automatic, {0, 1.1}}, GridLines -> Automatic, 
   GridLinesStyle -> LightGray, PlotStyle -> Red, 
   AxesLabel -> {"x", "u(x,t"}],
 {{t0, 0.01, "time"}, 0, 20, 0.001, Appearance -> "Labeled"},
 TrackedSymbols :> {t0}
 ]