# Numerical Schemes and Analytical Solutions of Fickian Diffusion with Dirichlet and Neumann Boundaries

## Evaluating the performance of finite difference schemes.

## 1. Introduction

Suppose we have a system where concentration diffuses from a singular point across a one-dimensional physical domain, modeled by the parabolic partial differential equation

often referred to as Fick’s Second Law, where diffusivity is denoted by the constant *D*. Initial conditions can be expressed as

Consider the Dirichlet boundary conditions

as well as the Neumann boundary conditions

where the system is said to have reached completion when the concentration profile at a particular iteration first reaches a linear condition.

The differential equation and boundaries may be solved mathematically and achieve effective analytical models. However, numerical methods offer an alternative in approximating the differential equation’s solutions by adopting finite difference schemes for derivatives. Such schemes develop large algebraic matrices that are extremely efficient for computer algorithms to solve. Of the numerical solutions, we are primarily concerned with explicit, implicit and Crank-Nicholson schemes, derived from Euler Methods. MATLAB will be used to calculate the time needed for the concentration profile to reach completion as well as provide graphical representations. The truncation error and stability of difference schemes will be discussed.

## 2. Analytical Solutions

Non-homogeneous partial differential equations can be solved via Laplace transformations. Fick’s Second Law will first be solved for the Dirichlet conditions, followed by the Neumann conditions.

## 2.1 Dirichlet Boundary Solution

When the Laplace transformation is taken with respect to time, we arrive at

where *U *represents the transformed concentration variable and *s *represents the time substitute. Applying the initial condition, the equation simplifies to

and the transformed boundaries appear as

The general solution to homogeneous ordinary differential equations is

Substituting the transformed boundary conditions into the general solution allows the following calculations,

For

Taylor’s expansion of the denominator allows us to discretize the equation, thus viable for inverse Laplace transformation.

Given the following inverse Laplace transformation identity,

where *erfc(z) *represents the complementary error function, we arrive at

# 2.2 Neumann Boundary Solution

When solving for the Neumann boundary problem, the methodology is strikingly similar, with the obvious exception of the transformed boundary conditions. The initial transformed equation remains

while the transformed boundaries now appear as

Substituting the transformed boundary conditions into the general solution allows the following calculations,

The denominators are expressed as a Taylor series for discretization, hence

## 3. Finite Difference Schemes

Given our existing analytical models, it may appear unnecessary to derive numerical schemes that offer approximations rather than exact solutions. However, it must be noted that analytical methodology grows progressively futile with problems of greater complexity. Numerical schemes, namely finite differences for our purposes, are imperative to computing and analyzing higher order models of increasing intricacy. The following numerical analysis is applicable to a far greater range of problems than the previously derived analytical solutions. Discrete approximations simply limit the numerical solution to a finite number of points. As such, discrete approximations will result in sets of algebraic equations where certain discrete unknowns, or “nodes”, may be evaluated. A “mesh” is often used to visualize the locations of nodes, dictated by a distance-step, measuring the local distance between adjacent locations in space, and a time-step, measuring the local distance between adjacent iterations. Hence, derivatives with respect to time or space may be replaced with difference equations, written in terms of discrete values located on the mesh. Though a mesh is not used in this paper, its principles are applied to the derivation and analysis of the following finite difference schemes.

## 3.1 Explicit Scheme

Concerning the following numerical schemes, the discrete approximation of the given problem is based on varying finite difference schemes. The explicit scheme of Fick’s Second Law can be calculated through the Forward Euler method. Unless otherwise mentioned, assume *n *as the distance-step and *j *as the time-step.

Fick’s Second Law can, therefore, be rewritten as

In some cases, expressing the scheme in a matrix vector form is more convenient. Conversion to matrix form does not give the explicit scheme a significant advantage during computation. Regardless, the basic matrix form is given below but is elaborated further in the implicit and Crank-Nicolson scheme.

## 3.2 Implicit Scheme

The implicit scheme gives the advantage of unconditional stability, as proven in (4.3.2). It can be calculated through the Backward Euler method,

which must be converted into matrix form because the values of concentration are now coupled.

For the given Dirichlet conditions, an additional matrix must be added to the right-hand side.

The Neumann conditions require an adjustment to the left-hand side matrix concerning the coefficients.

The new coefficient of

must be included at *(N+1,N) *for a *(N+1, N+1) *matrix.

## 3.3 Crank-Nicolson Scheme

The Crank-Nicolson scheme proves to be useful because of its ability to combine the unconditional stability of an implicit method with the accuracy of a second-order scheme in space and time, as proven in (4.1.3).

Allow

with the Dirichlet boundaries accounted for in a similar fashion. Previously solved for in the simple implicit scheme, Neumann boundaries are included in the coefficient matrix on the left-hand side.

## 4. Local Truncation Error and Stability of Difference Schemes

## 4.1 Local Truncation Error

Local truncation error pertains to the error in approximating differential equations, calculated by the difference between said differential equation and its finite difference representation at a particular point in space and time. From such calculations, the local accuracies of various difference schemes can be compared.

## 4.1.1 Local Truncation Error of Explicit Scheme

Taylor’s expansion permits the following definitions,

Therefore,

Given

the main portion of the local truncation error is

Thus,

A simplified and perhaps more elegant methodology is shown below.

## 4.1.2 Local Truncation Error of Implicit Scheme

Taylor’s expansion provides the following,

Substituting (3), (5) and (6) into (4) yields

Once again,

and the main portion of the local truncation error is

Thus, we arrive at the same result as in the explicit scheme,

## 4.1.3 Local Truncation Error of Crank-Nicolson Scheme

Substituting (1–3), (5), and (6) into (4) yields

Again,

and the main portion of the local truncation error is,

Thus, the local truncation error for this scheme can be expressed by

## 4.2 Stability and Convergence of Difference Schemes

Finite difference schemes vary in their stability. For the case of difference schemes of differential equations, a scheme may be deemed stable if it does not allow the magnification of round-off error or minor fluctuations in initial data as the iterative process continues.

Convergence can be defined as when a discretized equation approaches the analytical solution of the differential equation as the time and distance-step approach zero. A mathematical representation of convergence is shown below.

where

are fixed values.

## 4.3 Von Neumann Stability Analysis

The von Neumann stability analysis is a method used to verify the stability of finite difference schemes, particularly when applied to linear partial differential equations. The method is possible through the decomposition of round-off errors into a Fourier series. Recall that Fick’s Second Law may be discretized into

Allow round-off error to be defined by

where

represents the ideal case of the discretized equation calculated without round-off error and

represents the numerical solution that suffers from round-off error. If

is assumed to satisfy the discretized equation, we have

Thus, the variation of round-off error can be expressed in the finite Fourier series

where

as usual. The amplitude of error can be assumed as an exponential function of time,

Since the behavior of each term in the difference equation is the same as the series, it is appropriate to express the growth of round-off error of a single term as

Amplification, often referred to as the growth factor, must be equal or less than 1 if the system is to remain stable. For example,

## 4.3.1 Explicit Scheme

The Forward Euler Method provides the iterative equation

Since

Because

The greatest value of alpha that satisfies the inequality is 0.5. Thus, the explicit scheme is conditionally stable so long

## 4.3.2 Implicit Scheme

The Backward Euler Method provides the iterative equation

Through very similar calculations as that of the explicit, we arrive at

Observe that because

will always hold true. Therefore, the implicit scheme is unconditionally stable.

## 4.3.3 Crank-Nicolson Scheme

Yet again, the von Neumann analysis of the Crank-Nicolson Scheme is very similar to that of the explicit.

Since

will always hold true. Therefore, the Crank-Nicolson scheme is unconditionally stable.

## 5. Numeric and Graphical Results

## 5.1 Analytical Results

The Newton-Raphson method allows for the calculation of increasingly accurate approximations of the solutions to the analytical expressions. For the Dirichlet and Neumann solutions, the initial approximation of

respectively, enables the Newton-Raphson method to determine successive approximations that converge to the following analytical solutions.

Where* D=1,*

and where *D=.5,*

## 5.2 Numerical Results

However, the difference schemes must be solved through iterative calculations. The following charts depicts the time in seconds for the concentration profile to reach a linear state. Alpha is fixed at 0.5 while various values of diffusivity are used.

Where *D=1,*

Where *D=0.5,*

## 5.3 Graphical Results

**Analytical — Neumann Boundary**

**Analytical — Dirichlet Boundary**

**Numerical — Neumann Boundary**

The legend indicates the time in seconds at each iteration.

**Numerical — Dirichlet Boundary**

## 6. Conclusion

The explicit, implicit, Crank-Nicolson and analytical solutions to Fick’s Second Law yielded reasonable results within our established boundaries. Among the finite difference schemes, the explicit scheme is perhaps the most simplistic but faces restrictions due to its conditionally stable nature. The explicit scheme’s local truncation error is of the same order as the implicit scheme. Though the implicit and Crank-Nicolson schemes were both proven to be unconditionally stable and therefore converging, the Crank-Nicolson is faster and more accurate given its higher order of local truncation error. The analytical solution, of course, yields the highest degree of accuracy. However, it’s computational speed is substantially slower than the algebraic sets that difference formulas require.

## Software Specifications & Miscellaneous Notes

Simulations were run and plotted on “MATLAB R2016b”, licensed for academic use.

LaTeX was used for producing scientific or mathematical symbols and expressions. Thus, the conversion to a Medium story resulted in odd spacing and slightly discolored images.

The rights of this report are strictly reserved to the author. Please reach out if you are interested in any referencing logistics or the MATLAB algorithms.