I’m using a standard Crank-Nicholson algorithm to solve this system of two coupled diffusion equations:

$$\dot{u} - \dot{v} = \frac{\partial}{\partial x} \left( \alpha(x) \frac{\partial u}{\partial x} \right) \label{1} \tag{1} $$

$$\dot{v} - \dot{u} = \frac{\partial}{\partial x} \left( \beta(x) \frac{\partial v}{\partial x} \right) \label{2} \tag{2} $$

where $\alpha$ and $\beta$ are fully determined functions of $x$. BCs are Dirichlet conditions only. I ordered the solution vector as $u_1, v_1, u_2, v_2, \ldots$ to obtain a band matrix.

The linear system is obtained by writing discretized $\eqref{1}$ then $\eqref{2}$ for the first node, then these two equations for the second node, and so on.

The results I get are not completely nonsensical, but oscillations are so huge that the data cannot be exploited. Basically, solutions will be rather OK for half of the time steps, and the other ones will be extremely low negative values with little physical sense. There is no improvement when changing $\Delta t$ or $\Delta x$ or trying different BCs.

I failed at understanding where the phenomenon comes from. Could the problem be ill-posed? Or does it have something to do with the way equations are written/ordered?


2 Answers 2


I think it might be ill-posed, since the time-dependent parts are linearly dependent. If you add your two time-dependent equations together, you get a time-independent equation:

$(\alpha(x)u_x)_x + (\beta(x) v_x)_x = 0.$

Do your initial conditions satisfy this equation?

  • $\begingroup$ I did not pay much attention to the initial conditions in the debug phase, so this condition might very well be unsatisfied! I will check. I was also surprised by the linearly dependent time-dependent parts... I guess this is the result of some kind of conservation law in the physical process (I did not derive the PDE themselves, they were established and discussed in some previous works. I am only dealing with the simulation part.) Thank you! $\endgroup$
    – user20536
    Commented May 31, 2016 at 10:12
  • $\begingroup$ Isn’t this equation automatically satisfied with constant IC ? $\endgroup$
    – user20536
    Commented May 31, 2016 at 10:40

I was curious about your problem and it was not difficult for me to test CN method for your two equations. Added - parts beginning with "added" were written afterwards to address the comment.

The point given by David is important. Here are three plots for initial condition (added - that do not fulfill the suggested condition $\alpha u_{xx} + \beta v_{xx}=0$), the first time step and the 5th one, computed by implicit Euler method, the blue curve is u, the other one is v:

Initial condition - blue is u, the other one is v

u and v after small first time step

u and v after 5 time steps

The diffusion coefficients $\alpha=\beta=1$ are constant. Added - The solution changes rapidly in the first time step for arbitrary small time step that is unusual. Afterwards standard numerical diffusion can be observed (the 5th time step).

The same computations with Crank-Nicolson are strange as you observed (the 1st and 5th time step):

The first time step with Crank-Nicolson

The 5th time step with Crank-Nicolson

I suppose someone can give you later comment on your PDEs directly, but from this simplified numerical analysis I would guess that you have no chance to obtain reasonable solution with Crank-Nicolson method (added - if your initial functions do not fulfill $(\alpha u_x)_x + (\beta v_x)_x=0$). I doubt if the numerical solution with implicit Euler method is reasonable, but if yes, it suggests that for initial condition that do not fulfill this PDE described by David, you obtain for arbitrary small time a jump in the solution, so you have some "infinite" speed of evolution. Added - The situation is even worse for explicit Euler methods as it is noted in the comment.

This is only heuristic answer, but at least it confirms your doubts about Crank-Nicolson method for your problem. Below is the code in Mathematica that produce the pictures.

Added - if your initial functions fulfill the required stationary PDE (but in the approximative form using e.g. finite difference method as in the provided code) then you may use implicit ($\theta=1$ in the code) or explicit ($\theta=0$) Euler method or Crank-Nicolson method ($\theta=0.5$) for time discretization, they will work because these methods preserve this property also for evolved numerical functions if it is fulfilled for initial functions.

Model data
t0 = 0.; T = 0.1; L = 10.; alpha = 1.; beta = 1.; 
u0[x_] := If[Abs[x - 3] <= 1., 1., 0.]; 
v0[x_] := If[Abs[x - 7] <= 1., -1., 0.]; 

Method's input
Nd = 5; M = 20; theta = 0.5; \[Tau] = T/Nd; h = L/(M + 1); 
Mat = SparseArray[{}, {2*M, 2*M}]; 
ps = Array[{}, {2*M}]; 

For[i = 1, i <= M, i++, If[i != 1, Mat[[i,i - 1]] = (-alpha)*\[Tau]*theta]; Mat[[i,i]] = h^2 + 2*alpha*\[Tau]*theta; 
   If[i != M, Mat[[i,i + 1]] = (-alpha)*\[Tau]*theta]; Mat[[i,i + M]] = -h^2; ]
For[i = M + 1, i <= 2*M, i++, If[i != M + 1, Mat[[i,i - 1]] = (-beta)*\[Tau]*theta]; Mat[[i,i]] = h^2 + 2*beta*\[Tau]*theta; 
   If[i != 2*M, Mat[[i,i + 1]] = (-beta)*\[Tau]*theta]; Mat[[i,i - M]] = -h^2; ]

solver = LinearSolve[Mat]; 

For[i = 0, i <= M + 1, i++, u[i, 0] = u0[i*h]; v[i, 0] = v0[i*h]]

For[n = 1, n <= Nd, n++, 
  For[i = 1, i <= M, i++, ps[[i]] = h^2*(u[i, n - 1] - v[i, n - 1]) + (1 - theta)*\[Tau]*alpha*(u[i - 1, n - 1] - 2*u[i, n - 1] + 
         u[i + 1, n - 1]); ps[[M + i]] = h^2*(v[i, n - 1] - u[i, n - 1]) + (1 - theta)*\[Tau]*beta*(v[i - 1, n - 1] - 2*v[i, n - 1] + 
         v[i + 1, n - 1])]; r = solver[ps]; u[0, n] = u[M + 1, n] = v[0, n] = v[M + 1, n] = 0.; 
   For[i = 1, i <= M, i++, u[i, n] = r[[i]]; v[i, n] = r[[M + i]]]; ]

n = 0; ListLinePlot[{Table[{i*h, u[i, n]}, {i, 0, M + 1}], Table[{i*h, v[i, n]}, {i, 0, M + 1}]}, PlotRange -> All]
n = 1; ListLinePlot[{Table[{i*h, u[i, n]}, {i, 0, M + 1}], Table[{i*h, v[i, n]}, {i, 0, M + 1}]}, PlotRange -> All]
n = 5; ListLinePlot[{Table[{i*h, u[i, n]}, {i, 0, M + 1}], Table[{i*h, v[i, n]}, {i, 0, M + 1}]}, PlotRange -> All]
  • $\begingroup$ Thank you very much for this very detailed answer. If this is really hopeless I may test the explicit method, just to try. $\endgroup$
    – user20536
    Commented May 31, 2016 at 10:06
  • $\begingroup$ If you initial functions do not fulfill the stationary PDE as suggested by David, it makes no sense to use explicit method. To address it I added several explanation to my answer. I suggest you "reformulate" your origin problem as one stationary PDE and one (any of two) nonstationary PDE. Good luck! $\endgroup$ Commented May 31, 2016 at 14:20
  • $\begingroup$ The reformulation you suggest is a great idea that I will try as soon as possible $\endgroup$
    – user20536
    Commented May 31, 2016 at 15:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.