Is this system of diffusion equations well-posed?

Question

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?

Answer 1

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?

Answer 2

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:

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):

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}]; 

Method
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]