Numerical Solution of linear Second Order Partial Differential Equation

Numerical Solution of linear Second Order Partial Differential Equation

Haruna Usman Idriss¹, *Ali Bulama Mamman² and Bashir Mai Umar³

^{1, 2, 3}Department of Mathematics, Federal University Gashua, Yobe State, Nigeria.

*Correspondence Author

DOI: https://doi.org/10.51584/IJRIAS.2025.10040081

Received: 10 April 2025; Accepted: 16 April 2025; Published: 20 May 2025

ABSTRACT

This study explores the analytical and numerical solutions of partial differential equations (PDEs), focusing on parabolic (heat). The first part presents their analytical solutions using initial and boundary conditions and delves into the finite difference method (FDM), discussing forward, backward, and central difference schemes. These methods are applied to numerically solve one- and two-dimensional heat. The Crank-Nicolson method, recognized for its unconditional stability, is employed to improve the accuracy of heat equation solutions, overcoming limitations of explicit and implicit schemes. We then analyze the performance, strengths, and weaknesses of FDM through numerical simulations of one-dimensional heat. Due to computational constraints, Crank-Nicolson for 1D simulation, was not executed. Results indicate that the implicit backward difference method demonstrates superior stability by allowing unrestricted step sizes compared to the explicit forward difference method. These findings contribute to a deeper understanding of numerical PDE solutions and stability considerations in computational mathematics.

Keywords: PDEs, parabolic (heat) Equation, Crank-Nicolson Method

Mathematics Subject Classification: xx, xx

INTRODUCTION

Many problems in physical phenomena, such as physics, applied science, and engineering, can be modelled mathematically with the help of techniques of partial differential equations (John H and Fink, 1992). When a function comprises two or more independent variables, the differential equation is aptly called partial differential equation. Since the functions of multitudinous variables are inherently more problematic compared those of one variable, partial differential equations may lead to some challenging tasks in numerical problems. In fact, finding numerical solutions to those problems requires a type of scientific calculation which needs the help of a computing system (Cheney and Kincaid, 2004).

Numerous physical phenomena, which include electrostatic problems, heat conduction, fluid dynamics, electrodynamics, gravitational potential, can be modelled mathematically using  partial differential equations (or PDEs) with a set of initial conditions or initial boundary conditions (Everstine, 2010). PDEs form the solid bedrock of many mathematical models related to chemical, physical, and biological phenomena, and more recently their application has diffused into the fields of financial forecasting, economics, image processing and others. (Morton and Meyers, 2005). The partial differential equation depends on two or more independent variables. These variables can be time and one or more coordinates in space or plane (Erwin, 1976).

In the study of PDE’s such as heat and wave equations, there are particular kinds of boundary condition commonly associated with above equations. For the heat equation in term of boundary condition, the initial value of the solutions are defined but in a bounded domain. The Dirichlet condition on the boundary of the domain takes positive time (t). However, the classical boundary problem in the case of the wave is equation is the Cauchy problem because it defines both the initial position and initial velocity at.

The stand for saying whether a boundary condition is appropriate for a particular PDE’s is physically difficult to understand. However, it can be explained by fundamental mathematical insight (Hadamard, 1923 as cited in Brezis and Browder, 1997).

The onset of finite difference techniques in numerical application began in in the early 1950s and the emergence of computers that presented a convenient framework for dealing with complex problems of science and technology gave impetus to their development. The theoretical result has been found during the last five decades in term of the accuracy, consistency, stability and convergence of the finite difference method for partial differential equations (Fadugba and Adegboyegun, 2013). Furthermore, finite difference methods can be used to solve partial differential equations. This is done by approximating the differential equations over the area of integration by a system of algebraic equations. They also add that the finite difference approximations happened to be one of the simplest and oldest to solve partial differential equations. It was known by L. Euler since (1707-1783) ca. 1768, in one-dimensional space and was probably extended to two-dimensional space by C. Runge in the years (1858-1927) ca. 1908.

MATERIALS AND METHODS

Formulation of Finite Difference

Figure 1.1

The Finite Difference Method (FDM) relies on Taylor’s theorem. For a function \( f \) with continuous derivatives and equal spacing \( h \) in the x-direction:

(1) \( f(x+h) = f(x) + h f'(x) + \frac{h^2}{2!}f”(x) + \cdots \)

(2) \( f(x-h) = f(x) – h f'(x) + \frac{h^2}{2!}f”(x) – \cdots \)

• \( f'(x) \) is the derivative evaluated at \( x \)
• There exists a truncation error term \( \mathcal{O}(h^2) \)

Combining equations (1) and (2):

(3) \( f'(x) = \frac{f(x+h) – f(x-h)}{2h} + \mathcal{O}(h^2) \)

Taylor Series and Forward Difference

Using the Taylor series expansion of \( f(x+h) \):

(4) \( f(x+h) = f(x) + h f'(x) + \frac{h^2}{2!}f”(x) + \cdots \)

Truncating the series gives the forward difference:

(5) \( f'(x) \approx \frac{f(x+h) – f(x)}{h} \)

With truncation error:

(6) \( \mathcal{O}(h) \)

Similarly, backward difference:

(7) \( f'(x) \approx \frac{f(x) – f(x-h)}{h} \)

With truncation error:

(8) \( \mathcal{O}(h) \)

Combining forward and backward difference:

(9) \( f'(x) \approx \frac{f(x+h) – f(x-h)}{2h} \)

Truncation error is \( \mathcal{O}(h^2) \), meaning quadratic convergence when spacing is constant.

(10) \( f”(x) \approx \frac{f(x+h) – 2f(x) + f(x-h)}{h^2} + \mathcal{O}(h^2) \)

1D Heat Equation

Consider the 1D parabolic heat equation:

(12) \( \frac{\partial u}{\partial t} = \alpha \frac{\partial^2 u}{\partial x^2} \)

Explicit Method (Forward Difference)

Using central difference in space:

(13a) \( \frac{\partial^2 u}{\partial x^2} \approx \frac{u_{i+1}^n – 2u_i^n + u_{i-1}^n}{\Delta x^2} \)

And forward difference in time:

(14a) \( \frac{\partial u}{\partial t} \approx \frac{u_i^{n+1} – u_i^n}{\Delta t} \)

Substituting gives:

(16) \( \frac{u_i^{n+1} – u_i^n}{\Delta t} = \alpha \frac{u_{i+1}^n – 2u_i^n + u_{i-1}^n}{\Delta x^2} \)

Rewriting:

(17c) \( u_i^{n+1} = u_i^n + \lambda (u_{i+1}^n – 2u_i^n + u_{i-1}^n) \)

Where \( \lambda = \alpha \frac{\Delta t}{\Delta x^2} \)

Matrix Form

(18) \( \mathbf{U}^{n+1} = A \mathbf{U}^n \)

Implicit Method (Backward Difference)

Backward in time:

(19) \( \frac{u_i^{n+1} – u_i^n}{\Delta t} = \alpha \frac{u_{i+1}^{n+1} – 2u_i^{n+1} + u_{i-1}^{n+1}}{\Delta x^2} \)

Rewriting and rearranging:

(21) \( -\lambda u_{i-1}^{n+1} + (1 + 2\lambda)u_i^{n+1} – \lambda u_{i+1}^{n+1} = u_i^n \)

Matrix Form

(22) \( A \mathbf{U}^{n+1} = \mathbf{U}^n \)

Crank-Nicolson Method

This method averages explicit and implicit methods:

(23) \( \frac{u_i^{n+1} – u_i^n}{\Delta t} = \alpha \frac{1}{2} \left[ \frac{u_{i+1}^{n} – 2u_i^n + u_{i-1}^{n}}{\Delta x^2} + \frac{u_{i+1}^{n+1} – 2u_i^{n+1} + u_{i-1}^{n+1}}{\Delta x^2} \right] \)

Simplified to tridiagonal matrix form:

(25b) \( A \mathbf{U}^{n+1} = B \mathbf{U}^n + \mathbf{b} \)

With coefficients:

(26) \( a = -\lambda/2, \quad b = 1 + \lambda, \quad c = -\lambda/2 \)

And

See the below grid:

Figure 1.2

Figure 1.2 mesh point for Crank-Nicolson method. The open circles are the unknowns at each time step and the filled circles are known representing the initial and boundary conditions.

RESULTS

This section has solved example of heat (explicit forward difference and implicit backward difference methods with an exclusion of Crank Nicolson) equations numerically only one dimension by using Finite Difference Method (FDM).

Heat Diffusion equation in one Dimension

Suppose the initial and boundary conditions for the heat equation (12) are given below:

Table 1. Using the Forward-difference with 0.5

Figure 1a

Figure 1b Using the Forward difference method with r=0.5

Table 2. Using the Forward-difference with 0.55

Figure 2a

Figure 2b

Table 3. Using Backward difference with y=0.69

Using Back ward difference with y=0.69

Figure 3a

Figure 3b

The formula used is in equation (3.7d). This FDM method is stable for y=0.69 which satisfies its criteria of being not restricted to half (1/2) and can be used to form accurate approximation to the solution  for and as given in Table 3. The consecutive rows and columns are given in Figure (4.3). A three-dimensional representation of data in table (3) given in figure (3a) and (3b).

DISCUSSION

Strengths and Weaknesses Of The Finite Difference Methods

This section will give a brief account on how strength and weak finite difference methods are based on analysis of the stability of forward difference (explicit method) and its criteria of being stable and give a proof for stability criteria. Also, the stability of the backward difference (implicit method) their truncation error and Crank Nicolson method as a remedy for short comings of the first two methods.

Stability Analysis for Forward Difference (Explicit Method)

The strange behaviour occurred in the previous heat simulation gave rise to a problem. The solution to partial differential equations by the implicit forward difference method, can take a good care of error amplification or magnification for practical step size. This happened to be a crucial and pivotal aspect of stable and efficient solution. Here, the maximum value of the ratio of our step sizes (y) is half (1/2) as in figure (1a) and (1b) but for proof see equation .. When the ratio value (y) exceeds half (1/2), the explicit forward difference is said to be conditionally stable. The reason is its stability depends on the choice of step sizes.  This method is first order accurate in time and second in space.

The discretized form of equation (17c) is one of the contributing factors for truncation error because of the approximations of the derivative of (12) and also error magnification due to the method itself. Von Newman stability analysis measures the error amplification or magnification. Looking more closely at what the finite difference method is doing help us to investigate this magnification. To have a stable method, step size must be selected in such a way that the amplification factor should not be larger than 1(Sauer T, 2014).

In an explicit scheme, the temperature at time n+1 depends explicitly on the temperature at time n.

Stability Analysis for Backward Difference (Implicit Method)

In order to improve the stability of the explicit forward difference method, we use implicit backward difference method (Euler method) in equation (20d). Even though the magnitude of the truncation error of the explicit in (17c) is similar to that of implicit of (21) with different matrix arrangement, the nature analysis of stability of the implicit backward difference is similar to explicit forward difference case. Although, the value of the ratio of our step sizes (y) is not restricted to half (1/2) for stability which is rather bigger than that of explicit method. Hence the implicit method is stable for all value of the ratio of step sizes (y). The backward difference method is unconditionally stable and is first order accurate in time and second space. In term of two-dimensional heat, the scheme for implicit method provides second order convergence because only a very few iterations per time were needed.

Crank Nicolson Method

This method is unconditionally stable and also second order accurate in respect of both time and space. Previously, in the heat equation we noticed that, the explicit method is some time stable and the implicit method is all the times stable. When stable they both have error of order. The step size  requires being fairly small to have good accuracy (Sauer T, 2014).

Crank Nicolson is suitable finite difference method for the heat equation because it is unconditionally stable and second order convergence. Deriving this method cannot be achieved directly due to the first partial derivative  appeared in the equation (Sauer T, 2014). Although the explicit method is unconditionally stable with a serious draw back, the time step is very small.

In improving accuracy of Crank Nicolson method, one may employ higher order multistep methods or implicit Runge-Kutta methods to improve the Crank Nicolson method in time. To improve accuracy of Crank Nicolson, one can use analogy with Nomerov’s method which is out of the confines of this method.

Observations

• The error of numerical solutions increases with number of steps.
• These errors are called Accumulative errors.
• Step size has strong effect on the accuracy on the finite difference method.
• Trade-off between computational effort and step size is an issue in any numerical technique such as FDM.

Stability of the Finite Difference method for the heat equation Consider the following approximation to the 1-D heat equation (Gordon D. Smith, 2004):

\[
\frac{u_i^{n+1} – u_i^n}{\Delta t} = \alpha \frac{u_{i+1}^n – 2u_i^n + u_{i-1}^n}{(\Delta x)^2}
\]

According to the notation used throughout this work:

\[
r = \frac{\alpha \Delta t}{(\Delta x)^2}
\]

Let \( u_i^n = \xi^n e^{i k x} \), then:
\[
\xi = 1 – 4r \sin^2\left(\frac{k \Delta x}{2}\right) \tag{4.1}
\]

The amplification factor:
\[
G = \xi = 1 – 4r \sin^2\left(\frac{k \Delta x}{2}\right) \tag{28}
\]

Taking modulus:
\[
|G| = |1 – 4r \sin^2\left(\frac{k \Delta x}{2}\right)| \tag{29}
\]

According to the trigonometric identity:
\[
1 – \cos(k \Delta x) = 2 \sin^2\left(\frac{k \Delta x}{2}\right)
\]

Therefore,
\[
G = 1 – 2r(1 – \cos(k \Delta x)) \tag{30}
\]

\[
G = 1 – 4r \sin^2\left(\frac{k \Delta x}{2}\right) \tag{31}
\]

For the sake of stability, we need \( |G| \leq 1 \), such that:

\[
-1 \leq 1 – 4r \sin^2\left(\frac{k \Delta x}{2}\right) \leq 1 \tag{32}
\]

The right inequality gives:
\[
1 – 4r \sin^2\left(\frac{k \Delta x}{2}\right) \leq 1 \Rightarrow 0 \leq 4r \sin^2\left(\frac{k \Delta x}{2}\right) \tag{33}
\]

This is always true. The left inequality:
\[
-1 \leq 1 – 4r \sin^2\left(\frac{k \Delta x}{2}\right) \Rightarrow r \leq \frac{1}{2 \sin^2\left(\frac{k \Delta x}{2}\right)} \tag{34}
\]

Since this condition must hold for all \( k \), the most stringent case is when:
\[
\sin^2\left(\frac{k \Delta x}{2}\right) = 1 \Rightarrow r \leq \frac{1}{2} \tag{35}
\]

By rearranging the inequality:
\[
\frac{\alpha \Delta t}{(\Delta x)^2} \leq \frac{1}{2} \tag{36}
\]

This gives the stability condition for the explicit finite difference method.

CONCLUSIONS

This research has introduced method for solving 1-D linear second order partial differential equation for Parabolic Equation using finite difference method (FDM). The FDM with three types as the forward, central and backward differences are used to solve heat equation in 1 Dimension. The unconditionally Crank Nicolson method is applied solve heat problems in both the dimensions. Strengths and weaknesses (in term of stability and accuracy) of these methods were checked by numerical solution.

The implicit method is more stable than explicit and the Crank Nicolson is the in term of stability and accuracy. The Crank Nicolson is unconditionally stable. In FDM, the step taking for solving is convergent and accurate. Successive over relaxation method is applied in elliptic equation to speed up the rate of convergence. When this method is applied, the number of iterations reduces drastically.

In this research, we consider numerical methods for PDEs and obtain the approximations for finite difference method.

ACKNOWLEDGMENTS

This research was funded by the Tertiary Trust Fund (TETFund) Research Scholarship.

Conflict of interest: The authors declare no conflicts of interest.

REFERENCES

1. Abdullah, Abdul Rahman Bin. “The study of some numerical methods for solving parabolic partial differential equations.” PhD diss., Loughborough University, 1983.
2. Data Engineering and Applications: Proceedings of the International Conference, IDEA 2K22, Volume 1.Singapore: Springer Nature Singapore, 2024.
3. De Weck, O., and I. Y. Kim. “Finite Element Method. Engineering Design and Rapid Prototyping.” Massachussetts Institute of Technology (2004).
4. Dhatt, Gouri, Emmanuel Lefrançois, and Gilbert Touzot. Finite element method. John Wiley & Sons, 2012.
5. Hyman, James M., and Bernard Larrouturou. “The numerical differentiation of discrete functions using polynomial interpolation methods.” Applied Mathematics and Computation 10 (1982): 487-506.
6. Johnson, Claes. Numerical solution of partial differential equations by the finite element method. Courier Corporation, 2009.
7. Liu, Gui-Rong, and Siu Sin Quek. The finite element method: a practical course. Butterworth-Heinemann, 2013.
8. Mathews, John H., and Kurtis D. Fink. Numerical methods using MATLAB. 4. Upper Saddle River, NJ: Pearson prentice hall, 2004.
9. Pozrikidis, C. “Finite and spectral element methods using Matlab.” University of California, San Diego (2005).
10. Rao, Singiresu S. The finite element method in engineering. Elsevier, 2010.
11. Smith, Gordon D. Numerical solution of partial differential equations: finite difference methods. Oxford university press, 1985.