Crank–Nicolson method
Example: 1D diffusion with advection for steady flow, with multiple channel connections

This is a solution usually employed for many purposes when there is a contamination problem in streams or rivers under steady flow conditions, but information is given in one dimension only. Often the problem can be simplified into a 1-dimensional problem and still yield useful information.

Here we model the concentration of a solute contaminant in water. This problem is composed of three parts: the known diffusion equation (D x {\displaystyle D_{x}} chosen as constant), an advective component (which means that the system is evolving in space due to a velocity field), which we choose to be a constant U x {\displaystyle U_{x}}, and a lateral interaction between longitudinal channels (k {\displaystyle k}):

where C {\displaystyle C} is the concentration of the contaminant, and subscripts N {\displaystyle N} and M {\displaystyle M} correspond to previous and next channel.

The Crank–Nicolson method (where i {\displaystyle i} represents position, and j {\displaystyle j} time) transforms each component of the PDE into the following:

Now we create the following constants to simplify the algebra:

and substitute (2), (3), (4), (5), (6), (7), α {\displaystyle \alpha }, β {\displaystyle \beta } and λ {\displaystyle \lambda } into (1). We then put the new time terms on the left (j + 1 {\displaystyle j+1}) and the present time terms on the right (j {\displaystyle j}) to get

To model the first channel, we realize that it can only be in contact with the following channel (M {\displaystyle M}), so the expression is simplified to

In the same way, to model the last channel, we realize that it can only be in contact with the previous channel (N {\displaystyle N}), so the expression is simplified to

To solve this linear system of equations, we must now see that boundary conditions must be given first to the beginning of the channels:

For the last cell of the channels (z {\displaystyle z}), the most convenient condition becomes an adiabatic one, so

This condition is satisfied if and only if (regardless of a null value)

Let us solve this problem (in a matrix form) for the case of 3 channels and 5 nodes (including the initial boundary condition). We express this as a linear system problem:

where

Now we must realize that AA and BB should be arrays made of four different subarrays (remember that only three channels are considered for this example, but it covers the main part discussed above):

where the elements mentioned above correspond to the next arrays, and an additional 4×4 full of zeros. Please note that the sizes of AA and BB are 12×12:

The d vector here is used to hold the boundary conditions. In this example it is a 12×1 vector:

To find the concentration at any time, one must iterate the following equation:
startend:end_checker