Crank-Nicolson solver for heat equation

Question

My ultimate goal is to solve the 1D radial diffusion equation

$$\frac{\partial u(t,x)}{\partial x}=x^2\frac{\partial}{\partial x} \left(\frac{D(t)}{x^2} \frac{\partial u(t,x)}{\partial x } \right)$$ with a time-dependent diffusion coefficient $D(t)$, but as a first step I want to solve the simple 1D diffusion equation $$\frac{\partial u(t,x)}{\partial x}=D\frac{\partial^2 u(t,x)}{\partial x^2}$$

with a constant diffusion coefficient $D$ (K in code), which I'll set equal to 1.

The initial condition is $$u(0,x)=\sqrt{\frac{2}{\pi L_t}}e^{\frac{-x^2}{L_t^2}}$$ and here I also set the scale parameter $L_t=1$ (Lt = 1 in code) for now.

The boundary conditions are $$\left\{ \begin{aligned} {u(t,-10)=0\\ u(t,10)=0} \end{aligned} \right.$$

I want to use a Crank-Nicolson solver and I've used the code given here.

I've done some small adjustments, for example added an option for the MaxStepSize and my complete code reads as follows. For solving the equation I've adapted the example given (behind the same link) for solving the wave equation.

Clear["Global`*"]

(*Define options and method initialization for the Crank–Nicolson*)

Options[CrankNicolson] = {MaxIterations -> 10, Tolerance -> Automatic,
    MaxStepSize -> Automatic};
CrankNicolson /: 
 NDSolve`InitializeMethod[CrankNicolson, stepmode_, sd_, rhs_, state_,
   OptionsPattern[CrankNicolson]] := 
 Module[{prec, rtol, maxit}, maxit = OptionValue[MaxIterations];
  prec = state@"WorkingPrecision";
  rtol = OptionValue[Tolerance];
  If[rtol === Automatic, rtol = 10^(-prec*3/4)];
  CrankNicolson[maxit, rtol]]

(*Define the step function:*)

CrankNicolson[maxit_, rtol_]["Step"[f_, h_, t0_, x0_, f0_]] := 
  Module[{J, LU, t1 = t0 + h, x1, f1, residual, err, done = False, 
    tol = rtol, count = 0},
   x1 = x0 + h f0;
   f1 = f[t1, x1];
   x1 = x0 + (h/2) (f0 + f1);
   J = f["JacobianMatrix"[t1, x1]];
   LU = IdentityMatrix[Length[x1], SparseArray] - (h/2) J;
   LU = LinearSolve[LU];
   While[(count <= maxit) && ! done, f1 = f[t1, x1];
    residual = x1 - x0 - (h/2)*(f0 + f1);
    err = Norm[residual, Infinity];
    If[err < tol, done = True
     (*else*), x1 = x1 - LU[residual];
     count++;
     ]
    ];
   If[count > maxit, Message[CrankNicolson::cvmit, maxit];
    x1 = $Failed];
   {x1, f1}
   ];

(*This specifies for NDSolve the inputs, outputs, difference order, and 
  step mode for the Crank–Nicolson method:*)

CrankNicolson[___]["StepInput"] = {"F"["T", "X"], "H", "T", "X", "XP"};
CrankNicolson[___]["StepOutput"] = {"X", "XP"};
CrankNicolson[___]["DifferenceOrder"] := 2;
CrankNicolson[___]["StepMode"] := "Fixed";

(*Model paramaters*)
K = 1;
Lt = 1;
tmax = 10^2;
msf = 0.001;
k = 0.01;
xlim = 10;

(*Run the solver*)
uF = First[
   u /. NDSolve[{D[u[t, x], t] == K*D[u[t, x], x, x], 
      u[0, x] == Sqrt[2/(Pi*Lt)]*Exp[-x^2/Lt^2], u[t, -xlim] == 0, 
      u[t, xlim] == 0}, u, {t, 0, 1*tmax}, {x, -xlim, xlim}, 
     Method -> {"DoubleStep", 
       Method -> {CrankNicolson, 
         MaxStepSize -> {Automatic, 10^-7}}}]];

(*Plot results in 2D*)
Plot3D[
 Evaluate[uF[t, x]], {t, 0, tmax}, {x, -xlim, xlim}, 
 AxesLabel -> {"t", "x", "u"}, PlotRange -> All, Mesh -> All]

I now have two problems. First, the distribution should approach a uniform one in the spatial dimension as $t \rightarrow \infty$. However, when I look at the plot given by the code, the distribution seems to approach zero. I don't understand this, since in the diffusion equation I have no loss term, so the integral over $x$ should be constant, if I'm, not mistaken. However, this isn't the case right now. I've tried to decrease the MaxStepSize down to 10^-8 in both dimensions, but the result doesn't seem to change or the computation aborts. Below is the resulting graph I get, but I think that the "plateu" should be higher, approaching about 0.07 from above to conserve the spatial integral.

A second problem is that when I set tmax = 600 or larger, the solution starts to show a distribution that looks really strange, and setting tmax = 10^4 the behavior is clearly visible.

Answer 1

Let me extend my comment to an answer. 3 issues here.

1. The boundary conditions are incorrect

… the distribution seems to approach zero. I don't understand this, since in the diffusion equation I have no loss term, so the integral over x should be constant, if I'm, not mistaken. "

Sadly, you're mistaken. Indeed, there's no loss term inside the PDE, but you've set Dirichlet b.c. u[t, -xlim] == 0, u[t, xlim] == 0. This is amount to cooling down the whole domain with 2 walls whose temperatures are 0. You need to set adiabatic boundary condition i.e. Neumann b.c. whose RHS is 0:

K = 1;
Lt = 1;
tmax = 10^2;
msf = 0.001;
k = 0.01;
xlim = 10;

eq = D[u[t, x], t] == K*D[u[t, x], x, x];
ic = u[0, x] == Sqrt[2/(Pi*Lt)]*Exp[-x^2/Lt^2];
bc = D[u[t, x], x] == 0 /. {{x -> -xlim}, {x -> xlim}}

uFref = NDSolveValue[{eq, ic, bc}, 
    u, {t, 0, 1*tmax}, {x, -xlim, xlim}]; // AbsoluteTiming
(* {0.0242869, Null} *)

Animate[Plot[uFref[t, x], {x, -xlim, xlim}, AxesLabel -> {"x", "u"}, 
  PlotRange -> {0, 0.8}], {t, 0, 50}]

As we can see, now the "plateu" approaches about 0.07 as expected.

You can use CrankNicolson of course. The result looks the same, and the corresponding timing is about 0.117249 second. I'm not sure why you need Crank-Nicolson, but do remember the default ODE solver of NDSolve is quite robust and efficient, if it fails to solve a problem, classical method like RK4, Crank-Nicolson, etc. probably won't help either. ODE solver of NDSolve should always be the last thing to adjust.

2. The MaxStepSize doesn't have any effect

I've done some small adjustments, for example added an option for the MaxStepSize…

But this MaxStepSize in Options[CrankNicolson] is meaningless. Just try:

uF1 = NDSolveValue[{eq, ic, bc}, u, {t, 0, tmax}, {x, -xlim, xlim}, 
   Method -> {DoubleStep, 
     Method -> {CrankNicolson, MaxStepSize -> {Automatic, 10^-7}}}];

uF2 = NDSolveValue[{eq, ic, bc}, u, {t, 0, tmax}, {x, -xlim, xlim}, 
   Method -> {DoubleStep, Method -> {CrankNicolson, MaxStepSize -> ahhhhhh}}];

uF1 === uF2
(* True *)

Why? Because you've misunderstood how the Option of a function works. I know beginners may guess that, once Options[f]={…} is set, the function f will be clever enough to use those options in certain miraculous way, but it's not true. If you read the source code of CrankNicolson carefully, you'll see the option values of MaxIterations and Tolerance are explicitly refered with OptionValue!

"OK, then what's the correct way to set MaxStepSize for CrankNicolson?" There doesn't seem to be a simple way, because step size control isn't coded in the refered CrankNicolson, but MaxStepSize is already an option of NDSolve, why not use it? For example:

uF = 
  NDSolveValue[{eq, ic, bc}, u, {t, 0, tmax}, {x, -xlim, xlim}, 
   MaxStepSize -> {Automatic, 2 xlim/200}, 
   Method -> {DoubleStep, Method -> {CrankNicolson}}];

Dimensions /@ uF["Coordinates"]
(* {{30}, {201}} *)

I've used the undocumented "Coordinates" argument to extract grid coordinates from InterpolatingFunction[…]. To learn more about extracting data from InterpolatingFunction[…], check this and this post.

3. StiffnessTest should be turned off

Now, if you try CrankNicolson with tmax = 10^4, you'll see

NDSolveValue::ndstf At t == 111.35824722214431`, system appears to be stiff. Methods Automatic, BDF, or StiffnessSwitching may be more appropriate.

Observing the solution around t == 111, it doesn't seem to be stiff, so the stiffness detector of DoubleStep hasn't done a good job here. Then how to fix? I know little about stiffness detection, but according to my limited knowledge, Crank-Nicolson, which is an implicit method, doesn't need it, so let's turn it off:

tmax = 10^4; 
uFlong = 
 NDSolveValue[{eq, ic, bc}, u, {t, 0, tmax}, {x, -xlim, xlim}, 
  Method -> {DoubleStep, Method -> {CrankNicolson, MaxIterations -> 40}, 
    StiffnessTest -> False}];

Alternatively we can build this into CrankNicolson. Just execute

CrankNicolson[___]["StiffMethodQ"] := True;

The usage of "StiffMethodQ" is mentioned in tutorial "DoubleStep" Method for NDSolve. Now we don't need StiffnessTest option:

tmax = 10^4; 
uFlong2 = 
 NDSolveValue[{eq, ic, bc}, u, {t, 0, tmax}, {x, -xlim, xlim}, 
  Method -> {DoubleStep, 
             Method -> {CrankNicolson, MaxIterations -> 40}}]; // AbsoluteTiming
(* {5.92179, Null} *)

For comparison, timing of default method is 0.0297851 second.