Contrasting Numerical Schemes and Analytical Solutions of Fickian Diffusion with Dirichlet and Neumann Boundaries
Evaluating the performance of finite difference schemes.
--
1. Introduction
Analytical and numerical solutions differ drastically in derivation, efficiency, and implementation. This piece aims to provide a full-scale comparison between such models, from derivation to implementation, for a basic standard diffusion model.
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 to 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.
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
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.
