Abstract

Translate

The paper expands the finite difference method to the complex plane, and thus obtains an improvement in the resolution of differential equations with an increase in numerical precision and a generalization in the mathematical modeling of problems. The article begins with a selection of the best techniques for obtaining finite difference coefficients for approximating derivatives in the real domain. Then, the calculation is expanded to the complex domain. The research expands forward, backward, and central difference approximations of the real case by a quadrant approximation in the complex plane, which facilitates the use in boundary conditions of differential equations. The article shows many real and complex finite difference equations with their respective order of error, intended to serve as a basis and reference, which have been tested in practical examples of solving differential equations used in engineering. Finally, a comparison is made between the real and complex techniques of finite difference methods applied in the Theory of Elasticity. As a surprising result, the article shows that the finite difference method has great advantages in numerical precision, diversity of formulas, and modeling generalities in the complex domain when compared to the real domain.

Full text

Turn on search term navigation

Translate

1. Introduction

The finite difference method can be employed to solve boundary value or initial value problems involving ordinary or partial differential equations. Thus, this method can be applied to solve the equations of models using concentrated or distributed parameters [1,2,3]. The technique consists of replacing each derivative or differential of the differential equations by approximation of finite differences or finite addition of the variables, as shown in Equation (1) below:

(1) $\{\begin{matrix} d x ≅ ∆ x, d y ≅ ∆ y, d u ≅ ∆ u \\ \frac{d y}{d x} ≅ \frac{∆ y}{∆ x} ≅ \frac{\nabla y}{\nabla x} ≅ \frac{δ y}{δ x} ≅ \frac{μ δ y}{μ δ x}, \frac{d^{2} y}{d x^{2}} ≅ \frac{∆^{2} y}{∆ x^{2}} ≅ \frac{\nabla^{2} y}{\nabla x^{2}} ≅ \frac{δ^{2} y}{δ x^{2}}, \frac{d^{3} y}{d x^{3}} ≅ \frac{∆^{3} y}{∆ x^{3}} ≅ \frac{\nabla^{3} y}{\nabla x^{3}} ≅ \frac{{μ δ}^{3} y}{δ x^{3}} ≅ \frac{{μ δ}^{3} y}{μ δ x^{3}} \\ \frac{\partial u}{\partial x} ≅ \frac{∆ u}{∆ x} ≅ \frac{\nabla u}{\nabla x} ≅ \frac{δ u}{δ x} ≅ \frac{μ δ u}{μ δ x}, \frac{\partial^{2} u}{\partial x^{2}} ≅ \frac{∆^{2} u}{∆ x^{2}}, \frac{\partial^{3} u}{\partial x^{3}} ≅ \frac{∆^{3} u}{∆ x^{3}}, \frac{\partial^{2} u}{\partial x \partial y} ≅ \frac{∆^{2} u}{\partial x \partial y}, \frac{\partial^{3} u}{\partial x \partial y} ≅ \frac{∆^{3} u}{\partial x^{2} \partial y} \end{matrix}$

The error is the difference in the absolute value of the exact value and of the approximate value (Error = |Value_exact—Value_approx|). In principle, shortening the calculation step reduces the number of errors. The ratio of step to error decrease is the order of error (n), indicated by O(hⁿ), where h = Δx = x_i₊₁₋x_i is the constant step of the independent variables [1,4]. Thus, if a numerical method is of error order O(h⁴), it means that if you shorten the step by half, the error decreases (1/2)⁴ = 1/16; that is, the error decreases approximately 16 times. If you reduce the step by 10 times, the error drops 10⁴ = 10,000 times. Therefore, the higher the error order of a method is, the more accurate the method. Most numerical methods contain errors; that is, they are approximate methods. The interesting question is to be able to reduce these errors to a level that meets the demand of the problem under investigation. Thus, this approach makes the method as accurate as one wishes [5].

The idea of replacing differentials with finite additions is derived from the origin of differential and integral calculus, where at the time of its creation, the idea of limit was not very developed. Thus, mathematicians approximated the differentials by very small positive values near zero. Therefore, the finite difference method is very simple and intuitive mathematics.

In recent years, the finite difference method has been losing ground to other numerical methods, such as the finite element method and the boundary element method. It is believed that one of the reasons for this trend is the difficulty of obtaining approximations of the derivatives by finite difference equations, especially for high-order derivatives with several variables, and finite difference formulas that have a high order of error. The present work attempts to correct this gap and starts with a review of numerical finite difference operators.

2. Review of Numerical Finite Difference Operators

The operator is a symbol that represents, in short, a set of operations to be performed on a variable [1]. The domain or field of definition of an operator is the set, to which the operator can be applied [6]. Numeric operators have as an operand or domain the images of a discrete (or tabulated) function, because numerical methods work mainly with tables of points rather than algebraic or analytic functions [4]. A discrete function (tabulated function) will be represented by the set of points {(x₁, y₁), (x₂, y₂), …, (x_n, y_n)} or {[x₁, u(x₁)], [x₂, u(x₂)], …, [x_n_, u(x_n)]}, where y_n= u(x_n), and h = Δx = x_i₊₁-x_i is a constant step of the independent variables (h = Δx = x₂ − x₁ = x₃ − x₂ = x₄ − x₃ = x_i₊₁ − x_i = … = x_n − x_n−₁). Table 1 defines the main numeric operators.

In addition, the differential operator is usually denoted by D (the differential operator (D) is not a numeric operator), where $D f (x) = \frac{d}{d x} f (x) = f ’ (x) o r D^{2} f (x) = \frac{d^{2}}{{d x}^{2}} f (x) = f ´ ´ (x) o r D^{3} f (x) = \frac{d^{3}}{{d x}^{3}} f (x) = f ´ ´ ´ (x) o r D^{4} f (x) = \frac{d^{4}}{{d x}^{4}} f (x) = f^{(4)} (x) o r \dots o r D^{n} f (x) = \frac{d^{n}}{{d x}^{n}} f (x) = f^{(n)} (x)$ . Additionally, as the spacing of the variable x is constant, h = Δx = ∇x = δx = μδ x.

The numerical operators can be successively applied, which is referred to as the potence of operators, analogous to the operation of algebraic potentiation. Table 2 shows the successive application of the numeric operators, where n is any number.

3. Review of Finite Difference Equations for Derivative Approximations of a Single Variable

Finite difference equations are used to approximate derivatives of any order at any point, as long as there is a sufficient number of points [2,7,8,9,10]. To calculate the numerical value of a derivative, finite differences use nearby points and their locations. The locations of these sampled points receive the name of finite difference stencil.

To determine the coefficients for a difference equation:

(2) $\frac{d^{4} u}{d x^{4}} \approx A u (x - 2 h) + B u (x - h) + C u (x) + D u (x + h) + E u (x + 2 h)$

The formula for Taylor expansion is:

(3) $u (b) = u (a) + \frac{u ´ (a)}{1!} (b - a) + \frac{u ´ ´ (a)}{2!} {(b - a)}^{2} + \dots + \frac{u^{(n)} (a)}{n!} {(b - a)}^{n} + \dots$

By expanding the first term, A of u(x − 2h), with b = x − 2h and a = x, we obtain:

(4) $\begin{matrix} A u (x) + A u ´ (x) (- 2 h) + A \frac{1}{2} u ´ ´ (x) {(- 2 h)}^{2} + A \frac{1}{6} u ´ ´ ´ (x) {(- 2 h)}^{3} + \\ A \frac{1}{24} u ´ ´ ´ ´ (x) {(- 2 h)}^{4} + A \frac{1}{120} u ´ ´ ´ ´ ´ (x) {(- 2 h)}^{5} + \dots \end{matrix}$

We then did the same for all the other terms in the numerator:

(5) $\begin{matrix} + B u (x) + B u ´ (x) (- h) + B \frac{1}{2} u ´ ´ (x) {(- h)}^{2} + B \frac{1}{6} u ´ ´ ´ (x) {(- h)}^{3} \\ + B \frac{1}{24} u ´ ´ ´ ´ (x) {(- h)}^{4} + B \frac{1}{120} u ´ ´ ´ ´ ´ (x) {(- h)}^{5} + \dots \\ + C u (x) \\ + D u (x) + D u ´ (x) (h) + D \frac{1}{2} u ´ ´ (x) {(h)}^{2} + D \frac{1}{6} u ´ ´ ´ (x) {(h)}^{3} \\ + D \frac{1}{24} u ´ ´ ´ ´ (x) {(h)}^{4} + D \frac{1}{120} u ´ ´ ´ ´ ´ (x) {(h)}^{5} + \\ + E u (x) + E u ´ (x) (2 h) + E \frac{1}{2} u ´ ´ (x) {(2 h)}^{2} + E \frac{1}{6} u ´ ´ ´ (x) {(2 h)}^{3} \\ + E \frac{1}{24} u ´ ´ ´ ´ (x) {(2 h)}^{4} + E \frac{1}{120} u ´ ´ ´ ´ ´ (x) {(2 h)}^{5} + \dots \end{matrix}$

Now, to choose A, B, C, D, and E such that the summation of all of these terms results in the cancellation of all u(x), u′(x), u″(x), and u‴(x) terms and in the u’’’’(x) coefficients summing to one, it is necessary to solve the linear system below:

(6) $\{\begin{matrix} A + B + C + D + E = 0 \\ - 2 A - B + D + 2 E = 0 \\ \begin{matrix} 4 A + B + D + 4 E = 0 \\ \begin{matrix} - 8 A - B + D + 8 E = 0 \\ 16 A + B + D + 16 E = \frac{24}{h^{4}} \end{matrix} \end{matrix} \end{matrix}$

These can be rewritten in general form as $[w i t h φ (0) = 1 a n d \forall n | n \neq 0, φ (n) = 0]$ :

(7) ${(- 2)}^{n} A + {(- 1)}^{n} B + C φ (n) + {(1)}^{n} D + {(2)}^{n} E = \frac{4!}{h^{4}} φ (n - 4) f o r 0 \leq n \leq 4$

Or in matrix form as:

(8) $[\begin{matrix} 1 & 1 & 1 & 1 & 1 \\ - 2 & - 1 & 0 & 1 & 2 \\ 4 & 1 & 0 & 1 & 4 \\ - 8 & - 1 & 0 & 1 & 8 \\ 16 & 1 & 0 & 1 & 16 \end{matrix}] [\begin{matrix} A \\ B \\ C \\ D \\ E \end{matrix}] = \frac{1}{h^{4}} [\begin{matrix} 0 \\ 0 \\ 0 \\ 0 \\ 24 \end{matrix}]$

Now, this matrix is invertible, and the solution is:

(9) $[\begin{matrix} A \\ B \\ C \\ D \\ E \end{matrix}] = \frac{1}{h^{4}} {[\begin{matrix} 1 & 1 & 1 & 1 & 1 \\ - 2 & - 1 & 0 & 1 & 2 \\ 4 & 1 & 0 & 1 & 4 \\ - 8 & - 1 & 0 & 1 & 8 \\ 16 & 1 & 0 & 1 & 16 \end{matrix}]}^{- 1} [\begin{matrix} 0 \\ 0 \\ 0 \\ 0 \\ 24 \end{matrix}] = \frac{1}{h^{4}} [\begin{matrix} 1 \\ - 4 \\ 6 \\ - 4 \\ 1 \end{matrix}]$

Substituting these values into the original equation, we arrive at the finite difference equation:

(10) $\frac{d^{4} u}{d x^{4}} \approx \frac{u (x - 2 h) - 4 u (x - h) + 6 u (x) - 4 u (x + h) + u (x + 2 h)}{h^{4}}$

It is noticed that this equation can be generalized to obtain the finite difference equation from any finite difference stencil considering the desired derivative. With a stencil, s, of length, N, and derivative order d < N, the coefficients, c’s, are obtained by the finite difference coefficients equation $[w i t h φ (0) = 1 a n d φ (a n y o t h e r v a l u e) = 0]$ :

(11) $s_{1}^{n} c_{1} + s_{2}^{n} c_{2} + \dots + s_{N}^{n} c_{N} = \frac{d!}{h^{d}} φ (n - d) f o r 0 \leq n \leq N - 1$

The solution to this equation can be written in matrix form as:

(12) $[\begin{matrix} c_{1} \\ ⋮ \\ c_{N} \end{matrix}] = \frac{1}{h^{d}} [\begin{matrix} s_{1}^{0} & \dots & s_{N}^{0} \\ ⋮ & ⋱ & ⋮ \\ s_{1}^{N - 1} & \dots & s_{N}^{N - 1} \end{matrix}] [\begin{matrix} 0 \\ ⋮ \\ \begin{matrix} d! \\ \begin{matrix} ⋮ \\ 0 \end{matrix} \end{matrix} \end{matrix}]$

The error order of finite difference equations is provided by the order of the Taylor expansion.

Using Equation (12) to find the coefficients of the approximations of the derivatives by finite differences is very practical and quick. Furthermore, solving the linear system in (12) generates very general finite difference formulas, with any discrete input points. The Taylor series-based approximations are in closed forms. It was shown that a finite difference approximation of a derivative of a function u(x) at a reference mesh point x = x₁ can be represented as:

(13) $\frac{d}{d x_{1}} u (x_{1}) \approx \frac{1}{∆ x_{1}} \sum_{k} c_{k} u_{i 1 + k}$

where the coefficients c_k and the iterator k are defined based on the order and the type of the approximations. In the following discussion, we will denote c_k as c_k^F, c_k^B, and c_k^C for forward, backward, and central difference approximations, respectively [11,12,13,14]. Working algebraically with the determinants of the linear system (12), explicit equations can be deduced for calculating the coefficients of the finite difference formulas. For a forward approximation, 0 ≤ k ≤ N, where N is order of the approximation, and:

(14) $c_{k}^{F} = \{\begin{matrix} - \sum_{j = 1}^{N} 1 / j, k = 0, \\ \frac{{(- 1)}^{k + 1} N!}{k (N - k)! k!}, 1 \leq k \leq N \end{matrix}$

For backward approximations, −N ≤ k ≤ 0, and c_k^B = −c_k^F, i.e., the coefficients are additive inverse of those for forward approximations. For central approximations, −N ≤ k ≤ N, and:

(15) $c_{k}^{C} = \{\begin{matrix} 0, k = 0, \\ \frac{{(- 1)}^{k + 1} {(N!)}^{2}}{k (N - k)! (N + k)!}, - N \leq k \leq N, k \neq 0 \end{matrix}$

These formulas in Equations (14) and (15) can be used to find the finite difference approximations of any type and order. Explicit formulas for the coefficients of finite difference approximations of first-degree derivatives have been derived mathematically from Taylor series [15]. Although it is very simple to use Equations (14) and (15) for finite difference coefficients, the linear system in (12) is also easy to use and is even more general, as it is possible to choose which discrete points will be used in the finite difference formulas. The use of the linear system in (12) will also allow the expansion of finite difference coefficients for the complex field.

Table A1, Table A2 and Table A3 in Appendix A show some finite difference formulas obtained for the first derivatives in forward, backward, and centered forms for various orders of error. To use a more general notation and for several variables, we chose to exchange y for u as the independent variable and x for x₁ or x₂ or x₃ or… for the dependent variable, as well as the index i by i₁ or i₂ or i₃ or… for the indices of each dependent variable. Forward finite difference formulas have terms with an index greater than or equal to zero. Backward finite difference formulas have indices less than or equal to zero. Centered finite difference formulas have terms with indices equally on the positive and negative sides.

The proposed method for finding finite difference coefficients consists of using Taylor expansion for derivatives of a single variable with any sequence of points, and numerical operators to generalize finite difference formulas for partial derivatives and higher-order derivatives. The objective is to expand the approximations of derivatives by finite differences involving several different sequences of points and increasing the options of formulas with a high order of error that can be used. Thus, generalizing the finite difference equations that can be applied to solve differential equations.

It is observed that in Table A3 in Appendix A, there is an approximation with order 30, which is an order much higher than any order of error obtained, for example, for the Runge–Kutta method, which currently only has knowledge up to order 14 and is still very difficult to obtain (Terry Feagin’s order 14 Runge–Kutta scheme) [16].

4. Finite Difference Equations for Derivative Approximations with Numeric Operators

First-derivative finite difference formulas can now be applied to develop other finite difference approximations for higher-order derivatives and partial derivatives of more than one independent variable with the same error order as the original formulas. Here, a structure of the operational calculation is employed.

For example, for the second-order centered derivative (16), one has:

(16) $\begin{matrix} \frac{d u (x_{1})}{d x_{1}} = \frac{- u_{i 1 - 1} + u_{i 1 + 1}}{2 ∆ x_{1}}, O (∆ x_{1}^{2}) \to \frac{d u (x_{1})}{d x_{1}} = (\frac{- E^{- 1} + E}{2 ∆ x_{1}}) u_{i 1} \\ \frac{d^{2} u (x 1)}{d x_{1}^{2}} = {(\frac{- E^{- 1} + E}{2 ∆ x_{1}})}^{2} u_{i 1} = (\frac{E^{- 2} - 2 + E^{2}}{4 ∆ x_{1}^{2}}) u_{i 1} = (\frac{E^{- 2} u_{i 1} - 2 u_{i 1} + E^{2} u_{i 1}}{4 ∆ x_{1}^{2}}) = \\ (\frac{u_{i 1 - 2} - 2 u_{i 1} + u_{i 1 + 2}}{4 ∆ x_{1}^{2}}), O (∆ x_{1}^{2}) \end{matrix}$

Alternatively, for a second-order centered third derivative (17), one has:

(17) $\begin{matrix} \frac{d u (x_{1})}{d x_{1}} = \frac{- u_{i 1 - 1} + u_{i 1 + 1}}{2 ∆ x_{1}}, O (∆ x_{1}^{2}) \to \frac{d u (x_{1})}{d x_{1}} = (\frac{- E^{- 1} + E}{2 ∆ x_{1}}) u_{i 1} \\ \frac{d^{3} u (x_{1})}{{d x 1}^{3}} = {(\frac{- E^{- 1} + E}{2 ∆ x_{1}})}^{3} u_{i 1} = (\frac{{- E}^{- 3} + 3 E^{- 1} - 3 E + E^{3}}{8 ∆ x_{1}^{3}}) u_{i 1} = \\ (\frac{{- E}^{- 3} u_{i 1} + 3 E^{- 1} u_{i 1} - 3 E u_{i 1} + E^{3} u_{i 1}}{8 ∆ x_{1}^{3}}) = \\ (\frac{- u_{i 1 - 3} + 3 u_{i 1 - 1} - 3 u_{i 1 + 1} + u_{i 1 + 3}}{8 ∆ x_{1}^{3}}), O (∆ x_{1}^{2}) \end{matrix}$

In another example, for a second-order centered partial derivative with second-order error (18), one has:

(18) $\begin{matrix} \frac{d u (x_{1})}{d x_{1}} = \frac{- u_{i 1 - 1} + u_{i 1 + 1}}{2 ∆ x_{1}}, O (∆ x_{1}^{2}) \to \frac{d u (x_{1})}{d x_{1}} = (\frac{- {E_{x 1}}^{- 1} + E_{x 1}}{2 ∆ x_{1}}) u_{i 1} \\ \frac{\partial^{2} u (x_{1}, x_{2})}{\partial x_{1} \partial x 2} = (\frac{- {E_{x 1}}^{- 1} + E_{x 1}}{2 ∆ x_{1}}) (\frac{- {E_{x 2}}^{- 1} + E_{x 2}}{2 ∆ x_{2}}) u_{i 1, i 2} = \\ (\frac{{E_{x 1}}^{- 1} {E_{x 2}}^{- 1} - {E_{x 1}}^{- 1} E_{x 2} - E_{x 1} {E_{x 2}}^{- 1} + {E_{x 1} E}_{x 2}}{4 ∆ x_{1} ∆ x_{2}}) u_{i 1, i 2} \\ = \frac{{E_{x 1}}^{- 1} {E_{x 2}}^{- 1} u_{i 1, i 2} - {E_{x 1}}^{- 1} E_{x 2} u_{i 1, i 2} - E_{x 1} {E_{x 2}}^{- 1} u_{i 1, i 2} + {E_{x 1} E}_{x 2} u_{i 1, i 2}}{4 ∆ x_{1} ∆ x_{2}} \\ = \frac{u_{i 1 - 1, i 2 - 1} - u_{i 1 - 1, i 2 + 1} - u_{i 1 + 1, i 2 - 1} + u_{i 1 + 1, i 2 + 1}}{4 ∆ x_{1} ∆ x_{2}}, O (∆ x_{1}^{2} ∆ x_{2}^{2}) \end{matrix}$

This same procedure could be performed without using the Shift E operator. The successive application of various operators is equivalent to plugging one formula into another, and for partial derivatives the operators act independently on different variables. In other words, a successive application of the operational calculation as the derivative operator, as follows:

$\begin{matrix} \frac{d u (x_{1})}{d x_{1}} = \frac{- u_{i 1 - 1} + u_{i 1 + 1}}{2 ∆ x_{1}}, O (∆ x_{1}^{2}) \\ \frac{\partial^{2} u (x_{1}, x_{2})}{\partial x_{1} \partial x_{2}} = \frac{\partial (\frac{\partial u (x_{1}, x 2)}{\partial x_{1}})}{\partial x_{2}} = \frac{\partial (\frac{- u_{i 1 - 1, i 2} + u_{i 1 + 1, i 2}}{2 ∆ x_{1}})}{\partial x_{2}} = \\ \frac{- \frac{(- u_{i 1 - 1, i 2 - 1} + u_{i 1 + 1, i 2 - 1})}{2 ∆ x_{1}} + \frac{(- u_{i 1 - 1, i 2 + 1} + u_{i 1 + 1, i 2 + 1})}{2 ∆ x_{1}}}{2 ∆ x_{2}} = \\ \frac{u_{i 1 - 1, i 2 - 1} - u_{i 1 - 1, i 2 + 1} - u_{i 1 + 1, i 2 - 1} + u_{i 1 + 1, i 2 + 1}}{4 ∆ x_{1} ∆ x_{2}}, O (∆ x_{1}^{2} ∆ x_{2}^{2}) \end{matrix}$

These are simple and straightforward manipulations, and the operator approach merely reproduces them using indices. There is nothing wrong with this procedure, but it is believed that for larger expressions and to facilitate understanding, as with automation, the use of numerical operators is much simpler and easier because it transforms the determination of derivative approximations with finite differences into quick algebraic manipulations of polynomials with the variable E (numeric operator). In other words, the development of finite difference formulas has become simpler than solving linear systems or evaluating formulas with sums and producers. The objective of this research is to provide alternatives of how to obtain the coefficients of approximations of derivatives by finite differences in a simple, quick, easy, and accurate way.

For another example, of a third-order centered derivative with two independent variables (19), one has:

(19) $\begin{matrix} \frac{d u (x_{1})}{d x_{1}} = \frac{- u_{i 1 - 1} + u_{i 1 + 1}}{2 ∆ x_{1}}, O (∆ x_{1}^{2}) \to \frac{d u (x_{1})}{d x_{1}} = (\frac{- {E_{x 1}}^{- 1} + E_{x 1}}{2 ∆ x_{1}}) u_{i 1} \\ \frac{\partial^{3} u (x 1, x_{2})}{\partial x_{1}^{2} \partial x_{2}} = {(\frac{- {E_{x 1}}^{- 1} + E_{x 1}}{2 ∆ x 1})}^{2} (\frac{- {E_{x 2}}^{- 1} + E_{x 2}}{2 ∆ x_{2}}) u_{i 1, i 2} \\ = (\frac{{- E}_{x 1}^{- 2} {E_{x 2}}^{- 1} + {E_{x 1}}^{- 2} E_{x 2} + 2 {E_{x 2}}^{- 1} {- 2 E}_{x 2} {- E}_{x 1}^{2} {E_{x 2}}^{- 1} + {E_{x 1}}^{2} E_{x 2}}{8 ∆ x_{1}^{2} ∆ x_{2}}) u_{i 1, i 2} \\ = \frac{{- E}_{x 1}^{- 2} {E_{x 2}}^{- 1} u_{i 1, i 2} + {E_{x 1}}^{- 2} E_{x 2} u_{i 1, i 2} + 2 {E_{x 2}}^{- 1} u_{i 1, i 2} {- 2 E}_{x 2} {u_{i 1, i 2} - E}_{x 1}^{2} {E_{x 2}}^{- 1} u_{i 1, i 2} + {E_{x 1}}^{2} E_{x 2} u_{i 1, i 2}}{8 ∆ x_{1}^{2} ∆ x_{2}} \\ = \frac{- u_{i 1 - 2, i 2 - 1} + u_{i 1 - 2, i 2 + 1} + 2 u_{i 1, i 2 - 1} - 2 u_{i 1, i 2 + 1} - u_{i 1 + 2, i 2 - 1} + u_{i 1 + 2, i 2 + 1}}{8 ∆ x_{1}^{2} ∆ x_{2}}, O (∆ x_{1}^{2} ∆ x_{2}^{2}) \end{matrix}$

Now, for a third-order centered partial derivative with second-order error and three independent variables (20), we have:

(20) $\begin{matrix} \frac{d u (x_{1})}{d x_{1}} = \frac{- u_{i 1 - 1} + u_{i 1 + 1}}{2 ∆ x_{1}}, O (∆ x_{1}^{2}) \to \frac{d u (x_{1})}{d x_{1}} = (\frac{- {E_{x 1}}^{- 1} + E_{x 1}}{2 ∆ x_{1}}) u_{i 1} \\ \frac{\partial^{3} u (x_{1}, x_{2}, x_{3})}{\partial x_{1} \partial x_{2} \partial x_{3}} = (\frac{- {E_{x 1}}^{- 1} + E_{x 1}}{2 ∆ x_{1}}) (\frac{- {E_{x 2}}^{- 1} + E_{x 2}}{2 ∆ x_{2}}) (\frac{- {E_{x 3}}^{- 1} + E_{x 3}}{2 ∆ x_{3}}) u_{i 1, i 2, i 3} = \\ \frac{(\begin{matrix} {E_{x 1}}^{- 1} {E_{x 2}}^{- 1} {E_{x 3}}^{- 1} + {E_{x 1}}^{- 1} {{E_{x 2}}^{- 1} E}_{x 3} + {E_{x 1}}^{- 1} E_{x 2} {E_{x 3}}^{- 1} - {{{E_{x 1}}^{- 1} E}_{x 2} E}_{x 3} {+ E_{x 1} E}_{x 2}^{- 1} {E_{x 3}}^{- 1} \\ - E_{x 1} {E_{x 2}}^{- 1} E_{x 3} - E_{x 1} E_{x 2} {E_{x 3}}^{- 1} + E_{x 1} {E_{x 2} E}_{x 3} \end{matrix})}{8 ∆ x_{1} ∆ x_{2} ∆ x_{3}} u_{i 1, i 2, i 3} = \\ \frac{(\begin{matrix} {- u}_{i 1 - 1, i 2 - 1, i 3 - 1} + u_{i 1 - 1, i 2 - 1, i 3 + 1} + u_{i 1 - 1, i 2 + 1, i 3 - 1} - u_{i 1 - 1, i 2 + 1, i 3 + 1} + u_{i 1 + 1, i 2 - 1, i 3 - 1} \\ - u_{i 1 + 1, i 2 - 1, i 3 + 1} - u_{i 1 + 1, i 2 + 1, i 3 - 1} + u_{i 1 + 1, i 2 + 1, i 3 + 1} \end{matrix})}{8 ∆ x_{1} ∆ x_{2} ∆ x_{3}}, \\ O (∆ x_{1}^{2} ∆ x_{2}^{2} ∆ x_{3}^{2}) \end{matrix}$

Thus, having finite difference formulas for the first derivatives of a given error order, one can obtain approximations for successive derivatives and for partial derivatives of several variables with that same error order. For this, it is enough to use $u_{i p + n} = E_{i p}^{n}$ and make algebraic manipulations with a remarkable product of multiplication and potentiation of the operator E; at the end of the analytic expansions, the inverse transformation of $E_{x 1}^{n 1} E_{x 2}^{n 2} E_{x 3}^{n 3} \dots E_{x p}^{n p} u_{i} = u_{i 1 + n 1, i 2 + n 2, i 3 + n 3, \dots, i p + n p}$ . Using computer algebra software (Maple 2023 Academic Edition^®) makes it even easier to obtain finite difference approximations of successive derivatives of any order and of any number of independent variables.

With the use of finite difference formulas for first derivative approximations, obtaining finite difference formulas for approximation of higher-order derivatives of a variable or of several variables becomes a work of the relation between two operators and purely algebraic manipulation of polynomials, which can be automatically performed by mathematical algebra software for any error order.

There is also the possibility of mixing approximations with different error orders for each independent variable; thus, we used a formula with order error two for x₁ and with order error three for x₂. One can also mix the types of formulas by merging forward approximations of one variable with centered and backward approximations of another variable. The most common approach is to use finite difference approximations with the same type (forward, backward, and central) and with the same order of error.

5. Finite Difference Coefficients for Complex Variables

Analytic functions form a very important special case of functions defined over a 2D complex plane. A function u(x₁) with u and x₁ are complex variables and said to be analytic if $\frac{d u}{d x_{1}} = \lim_{∆ x_{1} \to 0} \frac{u (x_{1} + ∆ x_{1}) - u (x_{1})}{∆ x_{1}}$ is uniquely defined, no matter from which direction in the complex plane Δx₁ approaches zero [17,18]. Any function u(x) that possesses a Taylor expansion at some x-location can be extended to an analytic function, with a vast range of further consequences [17,18,19,20].

Now, instead of having forward, backward, and central difference approximations on the real line, we have first quadrant (Q1), second quadrant (Q2), third quadrant (Q3), fourth quadrant (Q4), and central (C) on the complex plane. Furthermore, the spacings form square matrices in the complex plane, with the number of points being squared powers (4, 9, 16, 25, 36, 49, 64, 81, 100, 121, 144, 169, etc.). Thus, there are spacings, as shown in Figure 1, Figure 2 and Figure 3, where h = Δx and I $= \sqrt{- 1}$ .

Following the same scheme for the real case, the linear system generated by the Taylor series for the complex domain is solved. In the complex case, the spacing points of the complex plane are used to generate the linear system, as shown in item 3. Table A4, Table A5, Table A6, Table A7 and Table A8 in Appendix A show the formulas obtained for the first derivative and a single variable. It should be noted that the Cauchy–Riemann equations were not used in the limit of defining the complex derivative, no matter from which direction the complex plane approaches zero. The only different thing used was that the operations of addition, subtraction, multiplication, and division were now in the complex domain.

This simple transition from the real domain to the complex domain in obtaining the approximations of the derivatives through finite differences is very important because it shows that the finite difference method can possibly also be easily expanded to the domain of variables.

Furthermore, with the resolution of complex linear systems, it is possible to obtain finite difference formulas for higher-order derivatives since these derivatives appear in the Taylor expansion. The problem is that there is a reduction in the order of the error as the order of the derivative increases, but with the advantage of having few points in the formulas, as shown in Table A9 and Table A10 in Appendix A.

A generalization for calculating finite difference coefficients via only linear systems, as described in item 3 of this article, would be to exchange Equation (3) for the Taylor series in several variables. It is also noted that the variables are now in the complex domain [2,3,4,10]. The linear systems generated are larger, but it is now possible to calculate finite difference approximations for partial derivatives and high-order derivatives:

(21) $\begin{matrix} u (b_{1}, b_{2}, \dots, b_{d}) = \sum_{n_{1} = 0}^{\infty} \dots \sum_{n_{d = 0}}^{\infty} \frac{{(b_{1} - a_{1})}^{n_{1}} \dots {(b_{d} - a_{d})}^{n_{d}}}{n_{1}! n_{2}! \dots n_{d}!} (\frac{\partial^{n_{1} + n_{2} \dots + n_{d}} u}{\partial x_{1}^{n_{1}} \partial x_{2}^{n_{2}} \dots \partial x_{d}^{n_{d}}}) (a_{1}, a_{2}, \dots a_{d}) \\ = u (a_{1}, a_{2}, \dots, a_{d}) + \sum_{j = 1}^{d} \frac{\partial u (a_{1}, a_{2}, \dots, a_{d})}{\partial x_{j}} (b_{j} - a_{j}) \\ + \frac{1}{2!} \sum_{j = 1}^{d} \sum_{k = 1}^{d} \frac{\partial^{2} u (a_{1}, a_{2}, \dots, a_{d})}{\partial x_{j} \partial x_{k}} (b_{j} - a_{j}) (b_{k} - a_{k}) + \frac{1}{3!} \sum_{j = 1}^{d} \sum_{k = 1}^{d} \sum_{l = 1}^{d} \frac{\partial^{3} u (a_{1}, a_{2}, \dots, a_{d})}{\partial x_{j} \partial x_{k} \partial x_{l}} (b_{j} - a_{j}) (b_{k} - a_{k}) (b_{l} - a_{l}) + \dots \end{matrix}$

In contrast, this research advocates using the Taylor series and the resolution of linear systems only for calculating finite difference coefficients only for first derivatives and with a single variable. To calculate the coefficients of finite difference approximations of the partial derivatives of several variables and derivatives of higher order than the first, it is always proposed to use the numerical operator Shift E and multiplication of polynomials. The reason for this is to reduce the number of mathematical operations necessary to calculate the finite difference coefficients, and to guarantee the order of the error in the approximations of derivatives by finite differences.

6. Second- and Higher-Order Derivatives and Derivatives with Several Variables

Again, with the complex finite difference approximations to the first-order derivatives, one can easily calculate the higher-order derivatives and the derivatives of several variables using the numerical operator Shift E and multiplication of polynomials. The use of operational calculation is perhaps the safest way to calculate other approximations of derivatives by finite differences, guaranteeing the order of the error. In practical engineering applications, the option is usually to maintain the quadrant and order of the error in finite difference approximations.

A difference in the complex case is that complex potentiation will be used in the operators. For example, $E^{n} u (x_{i}) = u (x_{i + n})$ can be generalized to $u [x_{i} + (n + m I) ∆ x] = u (x_{i + n + m I}) = E^{n + m I} u (x_{i})$ ; then, for example, $E^{1 + I} u (x_{i}) = u [x_{i} + (1 + I) ∆ x]$ and $E^{2 - 3 I} u (x_{i}) = u [x_{i} + (2 - 3 I) ∆ x]$ . Here, a structure of the operational calculation is employed.

For example, for the second-order quadrant 1 derivative (22), one has:

(22) $\begin{matrix} \frac{d u (x_{1})}{d x_{1}} = \frac{2 u_{i 1 + I} + (1 - I) u_{i 1 + 1 + I} + (- 3 + 3 I) u_{i 1} - 2 u_{i 1 + 1}}{2 ∆ x_{1}}, O (∆ x_{1}^{3}) \to \\ \frac{d u (x_{1})}{d x_{1}} = (\frac{2 E^{I} + (1 - I) E^{1 + I} + (- 3 + 3 I) - 2 E}{2 ∆ x_{1}}) u_{i 1} \\ \frac{d^{2} u (x_{1})}{d x_{1}^{2}} = {(\frac{2 E^{I} + (1 - I) E^{1 + I} + (- 3 + 3 I) - 2 E}{2 ∆ x_{1}})}^{2} u_{i 1} \\ \frac{d^{2}}{d x_{1}^{2}} u (x_{1}) = \frac{1}{2 ∆ x_{1}^{2}} (2 u_{i 1 + 2 I} + (2 - 2 I) u_{i 1 + 1 + 2 I} - I u_{i 1 + 2 + 2 I} + (- 6 + 6 I) u_{i 1 + I} + 2 I u_{i 1 + 1 + I} \\ - (2 + 2 I) u_{i 1 + 2 + I} - 9 I u_{i 1} + (6 + 6 I) u_{i 1 + 1} - 2 u_{i 1 + 2}), O (∆ x_{1}^{3}) \end{matrix}$

Alternatively, for a second-order quadrant 1 derivative with two variables (23), one has:

(23) $\begin{matrix} \frac{d u (x_{1})}{d x_{1}} = \frac{2 u_{i 1 + I} + (1 - I) u_{i 1 + 1 + I} + (- 3 + 3 I) u_{i 1} - 2 u_{i 1 + 1}}{2 ∆ x_{1}}, O (∆ x_{1}^{3}) \to \\ \frac{d u (x_{1})}{d x_{1}} = (\frac{2 E^{I} + (1 - I) E^{1 + I} + (- 3 + 3 I) - 2 E}{2 ∆ x_{1}}) u_{i 1} \\ \frac{\partial^{2} u (x_{1}, x_{2})}{\partial x_{1} \partial x_{2}} = (\frac{2 E_{x 1}^{I} + (1 - I) E_{x 1}^{1 + I} + (- 3 + 3 I) - 2 E_{x 1}}{2 ∆ x_{1}}) (\frac{2 E_{x 2}^{I} + (1 - I) E_{x 2}^{1 + I} + (- 3 + 3 I) - 2 E_{x 2}}{2 ∆ x_{2}}) u_{i 1, i 2} \\ \frac{\partial^{2}}{\partial x_{1} \partial x_{2}} u (x_{1}, x_{2}) = \frac{1}{2 ∆ x_{1} ∆ x_{2}} (2 u_{i 1 + I, i 2 + I} + (1 - I) u_{i 1 + I, i 2 + 1 + I} \\ + (- 3 + 3 I) u_{i 1 + I, i 2} - 2 I u_{i 1 + I, i 2 + 1} \\ + (1 - I) u_{i 1 + 1 + I, i 2 + I} - I u_{i 1 + 1 + I, i 2 + 1 + I} \\ + 3 I u_{i 1 + 1 + I, i 2} - (1 + I) u_{i 1 + 1 + I, i 2 + 1} \\ + (- 3 + 3 I) u_{i 1, i 2 + I} + 3 I u_{i 1, i 2 + 1 + I} - 9 I u_{i 1, i 2} \\ + (3 + 3 I) u_{i 1, i 2 + 1} - 2 I u_{i 1 + 1, i 2 + I} \\ - (1 + I) u_{i 1 + 1, i 2 + 1 + I} + (3 + 3 I) u_{i 1 + 1, i 2} \\ - 2 u_{i 1 + 1, i 2 + 1}), O (∆ x_{1}^{3} ∆ x_{2}^{3}) \end{matrix}$

For another example, of a third-order quadrant 1 derivative with two independent variables (24), one has:

(24) $\begin{matrix} \frac{d u (x_{1})}{d x_{1}} = \frac{2 u_{i 1 + I} + (1 - I) u_{i 1 + 1 + I} + (- 3 + 3 I) u_{i 1} - 2 u_{i 1 + 1}}{2 ∆ x_{1}}, O (∆ x_{1}^{3}) \to \\ \frac{d u (x_{1})}{d x_{1}} = (\frac{2 E^{I} + (1 - I) E^{1 + I} + (- 3 + 3 I) - 2 E}{2 ∆ x_{1}}) u_{i 1} \\ \frac{\partial^{3} u (x_{1}, x_{2})}{\partial x_{1}^{2} \partial x_{2}} = {(\frac{2 E_{x 1}^{I} + (1 - I) E_{x 1}^{1 + I} + (- 3 + 3 I) - 2 E_{x 1}}{2 ∆ x_{1}})}^{2} (\frac{2 E_{x 2}^{I} + (1 - I) E_{x 2}^{1 + I} + (- 3 + 3 I) - 2 E_{x 2}}{2 ∆ x_{2}}) u_{i 1, i 2} \\ \frac{\partial^{3}}{{\partial x}_{1}^{2} \partial x_{2}} = \frac{1}{4 ∆ x_{1}^{2} ∆ x_{2}} (4 u_{i 1 + 2 I, i 2 + I} + (2 - 2 I) u_{i 1 + 2 I, i 2 + 1 + I} + (- 6 + 6 I) u_{i 1 + 2 I, i 2} \\ - 4 I u_{i 1 + 2 I, i 2 + 1} + (4 - 4 I) u_{i 1 + 1 + 2 I, i 2 + I} - 4 I u_{i 1 + 1 + 2 I, i 2 + 1 + I} \\ + 12 I u_{i 1 + 1 + 2 I, i 2} - (4 + 4 I) u_{i 1 + 1 + 2 I, i 2 + 1} - 2 I u_{i 1 + 2 + 2 I, i 2 + I} \\ - (1 + I) u_{i 1 + 2 + 2 I, i 2 + 1 + I} + (3 + 3 I) u_{i 1 + 2 + 2 I, i 2} - 2 u_{i 1 + 2 + 2 I, i 2 + 1} \\ + (- 12 + 12 I) u_{i 1 + I, i 2 + I} + 12 I u_{i 1 + I, i 2 + 1 + I} - 36 I u_{i 1 + I, i 2} \\ + (12 + 12 I) u_{i 1 + I, i 2 + 1} + 4 I u_{i 1 + 1 + I, i 2 + I} + (2 + 2 I) u_{i 1 + 1 + I, i 2 + 1 + I} \\ - (6 + 6 I) u_{i 1 + 1 + I, i 2} + 4 u_{i 1 + 1 + I, i 2 + 1} - (4 + 4 I) u_{i 1 + 2 + 2 I, i 2 + I} \\ - 4 u_{i 1 + 2 + I, i 2 + 1 + I} + 12 u_{i 1 + 2 + I, i 2} + (- 4 + 4 I) u_{i 1 + 2 + I, i 2 + 1} - 18 u_{i 1, i 2 + I} \\ - (9 + 9 I) u_{i 1, i 2 + 1 + I} + (27 + 27 I) u_{i 1, i 2} - 18 u_{i 1, i 2 + 1} \\ + (12 + 12 I) u_{i 1 + 1, i 2 + I} + 12 u_{i 1 + 1, i 2 + 1 + I} - 36 u_{i 1 + 1, i 2} \\ + (12 - 12 I) u_{i 1 + 1, i 2 + 1} - 4 u_{i 1 + 2, i 2 + I} + (- 2 + 2 I) u_{i 1 + 2, i 2 + 1 + I} \\ + (6 - 6 I) u_{i 1 + 2, i 2} + 4 I u_{i 1 + 2, i 2 + 1}), O (∆ x_{1}^{3} ∆ x_{2}^{3}) \end{matrix}$

Now, for a third-order quadrant 1 partial derivative with three independent variables (25), we have:

(25) $\begin{matrix} \frac{d u (x_{1})}{d x_{1}} = \frac{2 u_{i 1 + I} + (1 - I) u_{i 1 + 1 + I} + (- 3 + 3 I) u_{i 1} - 2 u_{i 1 + 1}}{2 ∆ x_{1}}, O (∆ x_{1}^{3}) \to \\ \frac{d u (x_{1})}{d x_{1}} = (\frac{2 E^{I} + (1 - I) E^{1 + I} + (- 3 + 3 I) - 2 E}{2 ∆ x_{1}}) u_{i 1} \\ \frac{\partial^{3} u (x_{1}, x_{2}, x_{3})}{\partial x_{1} \partial x_{2} \partial x_{3}} = (\frac{2 E_{x 1}^{I} + (1 - I) E_{x 1}^{1 + I} + (- 3 + 3 I) - 2 E_{x 1}}{2 ∆ x_{1}}) (\frac{2 E_{x 2}^{I} + (1 - I) E_{x 2}^{1 + I} + (- 3 + 3 I) - 2 E_{x 2}}{2 ∆ x_{2}}) \\ (\frac{2 E_{x 3}^{I} + (1 - I) E_{x 3}^{1 + I} + (- 3 + 3 I) - 2 E_{x 3}}{2 ∆ x_{3}}) u_{i 1, i 2, i 3} \\ \frac{\partial^{3}}{\partial x_{1} \partial x_{2} \partial x_{3}} = \frac{1}{4 ∆ x_{1} ∆ x_{2} ∆ x_{3}} (4 u_{i 1 + I, i 2 + I, i 3 + I} + (2 - 2 I) u_{i 1 + I, i 2 + I, i 3 + 1 + I} \\ + (- 6 + 6 I) u_{i 1 + I, i 2 + I, i 3} - 4 u_{i 1 + I, i 2 + I, i 3 + 1} + (2 - 2 I) u_{i 1 + I, i 2 + 1 + I, i 3 + I} \\ - 2 I u_{i 1 + I, i 2 + 1 + I, i 3 + 1 + I} + 6 I u_{i 1 + I, i 2 + 1 + I, i 3} - (2 + 2 I) u_{i 1 + I, i 2 + 1 + I, i 3 + 1} \\ + (- 6 + 6 I) u_{i 1 + I, i 2, i 3 + I} + 6 I u_{i 1 + I, i 2, i 3 + 1 + I} - 18 I u_{i 1 + I, i 2, i 3} \\ + (6 + 6 I) u_{i 1 + I, i 2, i 3 + 1} - 4 I u_{i 1 + I, i 2 + 1, i 3 + I} - (2 + 2 I) u_{i 1 + I, i 2 + 1, i 3 + 1 + I} \\ + (6 + 6 I) u_{i 1 + I, i 2 + 1, i 3} - 4 u_{i 1 + I, i 2 + 1, i 3 + 1} + (2 - 2 I) u_{i 1 + 1 + I, i 2 + I, i 3 + I} \\ - 2 I u_{i 1 + 1 + I, i 2 + I, i 3 + 1 + I} + 6 I u_{i 1 + 1 + I, i 2 + I, i 3} - (2 + 2 I) u_{i 1 + 1 + I, i 2 + I, i 3 + 1} \\ - 2 I u_{i 1 + 1 + I, i 2 + 1 + I, i 3 + I} - (1 + I) u_{i 1 + 1 + I, i 2 + 1 + I, i 3 + 1 + I} \\ + (3 + 3 I) u_{i 1 + 1 + I, i 2 + 1 + I, i 3} - 2 u_{i 1 + 1 + I, i 2 + 1 + I, i 3 + 1} + 6 I u_{i 1 + 1 + I, i 2, i 3 + I} \\ + (3 + 3 I) u_{i 1 + 1 + I, i 2, i 3 + 1 + I} - (9 + 9 I) u_{i 1 + 1 + I, i 2, i 3} + 6 u_{i 1 + 1 + I, i 2, i 3 + 1} \\ - (2 + 2 I) u_{i 1 + 1 + I, i 2 + 1, i 3 + I} - 2 u_{i 1 + 1 + I, i 2 + 1, i 3 + 1 + I} + 6 u_{i 1 + 1 + I, i 2 + 1, i 3} \\ + (- 2 + 2 I) u_{i 1 + 1 + I, i 2 + 1, i 3 + 1} + (- 6 + 6 I) u_{i 1, i 2 + I, i 3 + I} + 6 I u_{i 1, i 2 + I, i 3 + 1 + I} \\ - 18 I u_{i 1, i 2 + I, i 3} + (6 + 6 I) u_{i 1, i 2 + I, i 3 + 1} + 6 I u_{i 1, i 2 + 1 + I, i 3 + I} \\ + (3 + 3 I) u_{i 1, i 2 + 1 + I, i 3 + 1 + I} - (9 + 9 I) u_{i 1, i 2 + 1 + I, i 3} + 6 u_{i 1, i 2 + 1 + I, i 3 + 1} \\ - 18 I u_{i 1, i 2, i 3 + I} - (9 + 9 I) u_{i 1, i 2, i 3 + 1 + I} + (27 + 27 I) u_{i 1, i 2, i 3} \\ - 18 u_{i 1, i 2, i 3 + 1} + (6 + 6 I) u_{i 1, i 2 + 1, i 3 + I} + 6 u_{i 1, i 2 + 1, i 3 + 1 + I} - 18 u_{i 1, i 2 + 1, i 3} \\ + (6 - 6 I) u_{i 1, i 2 + 1, i 3 + 1} - 4 I u_{i 1 + 1, i 2 + I, i 3 + I} \\ - (2 + 2 I) u_{i 1 + 1, i 2 + I, i 3 + 1 + I} (6 + 6 I) u_{i 1 + 1, i 2 + I, i 3} - 4 u_{i 1 + 1, i 2 + I, i 3 + 1} \\ - (2 + 2 I) u_{i 1 + 1, i 2 + 1 + I, i 3 + I} - 2 u_{i 1 + 1, i 2 + 1 + I, i 3 + 1 + I} + 6 u_{i 1 + 1, i 2 + 1 + I, i 3} \\ + (- 2 + 2 I) u_{i 1 + 1, i 2 + 1 + I, i 3 + 1} + (6 + 6 I) u_{i 1 + 1, i 2, i 3 + I} + 6 u_{i 1 + 1, i 2, i 3 + 1 + I} \\ - 18 u_{i 1 + 1, i 2, i 3} + (6 - 6 I) u_{i 1 + 1, i 2, i 3 + 1} - 4 u_{i 1 + 1, i 2 + 1, i 3 + I} \\ + (- 2 + 2 I) u_{i 1 + 1, i 2 + 1, i 3 + 1 + I} + (6 - 6 I) u_{i 1 + 1, i 2 + 1, i 3} \\ + 4 I u_{i 1 + 1, i 2 + 1, i 3 + 1}), O (∆ x_{1}^{3} ∆ x_{2}^{3} ∆ x_{3}^{3}) \end{matrix}$

In this way, obtaining a finite difference formula becomes a simple multiplication of polynomials that can be automated with the use of computational tools.

7. Testing the Finite Difference Formulas Obtained and Checking the Error Order

It was observed that all real and complex finite difference formulas obtained could be numerically verified and tested, and even the order of the error could be numerically validated by simply choosing an easily derived function with an exponential and where the Ithprime function returns the ith prime number, where the first prime number is Ithprime(1) = 2, the second prime number is Ithprime(2) = 3, and the third prime number is Ithprime(3) = 5, for example:

$u = e^{[\frac{1000}{(1000 + I t h p r i m e (1))}] x_{1}} e^{[\frac{1000}{(1000 + I t h p r i m e (2))}] x_{2}} e^{[\frac{1000}{(1000 + I t h p r i m e (3))}] x_{3}} \dots$

We chose a value for x₁, x₂, x₃,…; for example, where each x₁, x₂, x₃, … is a different, randomly generated real or complex value (for complex finite difference formulas), and then chose a small value for the step or increment of Δx₁ or Δx₂ or Δx₃…; for example, where each Δx₁, Δx₂, Δx₃, … is a different randomly generated real or complex value (for complex finite difference formulas), with absolute value in the range of 0.1 and 0.01. We then compared the analytically obtained derivative using derivation rules, which is the exact value, with the values of the derivatives obtained using the finite difference formulas, which is the approximate value. We could then calculate the error, which is the absolute value of the difference between the exact value and the approximate value. Then, a smaller step value was selected, for example, do Δx₁ = Δx₁/10 or Δx₂ = Δx₂/10 or Δx₃ = Δx₃/10 or…, ten times smaller. The value of the derivative was recalculated using the finite difference formulas and the new error was obtained. With the values of the two errors with different steps, we could obtain the error order of the approximation by finite differences, making the logarithm of the ratio of errors on the basis of the ratio of the steps and rounding the result of the logarithm to the nearest integer [21]. Thus, we obtained:

(26) $\begin{matrix} S t e p 1 = 0.01 \to E r r o r 1 = |u_{e x a c t} - u_{a p r o x 1}| \\ S t e p 2 = 0.001 \to E r r o r 2 = |u_{e x a c t} - u_{a p r o x 2}| \\ E r r o r O r d e r = l = r o u n d ({l o g}_{(\frac{S t e p 1}{S t e p 2})} (\frac{E r r o r 1}{E r r o r 2})) = \\ r o u n d (\frac{l n (\frac{E r r o r 1}{E r r o r 2})}{l n (\frac{S t e p 1}{S t e p 2})}) \to O ({∆ x}^{l}) \end{matrix}$

Special care must be taken so that Error1 and Error2 are nonzero so as not to produce an error in the expression (26) above. This test checks not only whether the finite difference formula approximates the real value of the derivatives but also the order of their error. It is, therefore, an important test for the verification of the finite difference equations obtained [22]. In this analytical test, the finite difference derivative approximations in the complex domain showed exceptional performance, with a small absolute error and high precision and accuracy. Passing this test, it was suggested to verify whether the approximation formulas of the calculus of derivatives by finite differences are useful in practical problems of engineering and physics [23,24].

8. Approximations of Derivatives Using Complex Finite Differences for Non-Rectangular Grids

An interesting feature that appears in the complex domain is the use of non-rectangular or square grids. The use of this type of grid is important for problems with heterogeneous materials, complex geometries, or with different initial and boundary conditions.

Figure A1, Figure A2, Figure A3, Figure A4, Figure A5, Figure A6, Figure A7 and Figure A8 in Appendix B show examples of finite difference formulas for non-rectangular grids. These equations were deduced in a similar way to the others using complex linear systems obtained from the Taylor series for the grid points, Equations (3)–(5) and (12). Only the first-order derivatives were calculated. Higher-order derivatives and several variables could be obtained using polynomial multiplication with the Shift (E) operator. When using non-rectangular or square grids, care must be taken when applying the finite difference method with successive use in numbering the formula indices. Therefore, symmetric and regular grids are computationally easier to use.

9. Comparison of the Real Solution with the Complex One by Finite Differences in the Theory of Disk Elasticity in Compression

As the main subject of this work, the finite difference formulas obtained were used to solve elasticity theory problems, with fourth-order derivatives, of a ring in compression by a constant force, as shown in Figure 4, where the radius of the inner disk (R1) was 0.02 m, the radius of the outer disk (R2) was 0.04 m, and the constant force (P) was 100 N, with an angle α of 10 degrees and disc thickness of 0.0025 m. This is an interesting problem because in addition to the Dirichlet or a Neumann boundary condition, it presents a boundary condition using second derivatives for the normal and shear stresses [23,25].

In the real case, the problem can be mathematically modeled by Airy’s stress function (ϕ) as:

(27) $\frac{\partial^{4} ϕ}{{\partial x}^{4}} + 2 \frac{\partial^{4} ϕ}{{\partial x}^{2} {\partial y}^{2}} + \frac{\partial^{4} ϕ}{{\partial y}^{4}} = 0$

In polar coordinates (r² = x² + y², θ = arctan(y/x)):

$(\frac{\partial^{2}}{{\partial r}^{2}} + \frac{1}{r} \frac{\partial}{\partial r} + \frac{1}{r^{2}} \frac{\partial^{2}}{{\partial θ}^{2}}) (\frac{\partial^{2} ϕ}{{\partial r}^{2}} + \frac{1}{r} \frac{\partial ϕ}{\partial r} + \frac{1}{r^{2}} \frac{\partial^{2} ϕ}{{\partial θ}^{2}}) = 0$

We assumed that the two-dimensional Cartesian normal (σ_x and σ_y) and shear (τ_xy) stress components could be represented by a stress function, ϕ, such that:

(28) $σ_{x} = \frac{\partial^{2} ϕ}{{\partial y}^{2}} σ_{y} = \frac{\partial^{2} ϕ}{{\partial x}^{2}} τ_{x y} = - \frac{\partial^{2} ϕ}{\partial x \partial y}$

In polar coordinates, with radial stress, σ_r, circumferential stress, σ_θ, and polar shear stress, τ_rθ:

(29) $σ_{r} = \frac{1}{r} \frac{\partial ϕ}{\partial r} + \frac{1}{r^{2}} \frac{\partial^{2} ϕ}{{\partial θ}^{2}} σ_{θ} = \frac{\partial^{2} ϕ}{{\partial r}^{2}} τ_{r θ} = \frac{1}{r^{2}} \frac{\partial ϕ}{\partial θ} - \frac{1}{r} \frac{\partial^{2} ϕ}{\partial r \partial θ}$

The principal stresses yielded:

(30) $σ_{1} = \frac{σ_{x} + σ_{y}}{2} + \sqrt{{(\frac{σ_{x} - σ_{y}}{2})}^{2} + {τ_{x y}}^{2}} σ_{2} = \frac{σ_{x} + σ_{y}}{2} - \sqrt{{(\frac{σ_{x} - σ_{y}}{2})}^{2} + {τ_{x y}}^{2}}$

The maximum shear stresses:

(31) $τ_{m a x} = \{\begin{matrix} \frac{(σ_{1} - σ_{2})}{2}, & σ_{1} σ_{2} \leq 0 \\ \frac{σ_{1}}{2}, & σ_{1} > 0 a n d σ_{2} > 0 \\ \frac{|σ_{2}|}{2} & σ_{1} < 0 a n d σ_{2} < 0 \end{matrix}$

Expressions could be obtained that related the polar stress components to the Cartesian stress components, as follows [26]:

(32) $\{\begin{matrix} σ_{r} = σ_{x} {c o s}^{2} θ + σ_{y} {s i n}^{2} θ + τ_{x y} s i n 2 θ \\ σ_{θ} = σ_{y} {c o s}^{2} θ + σ_{x} {s i n}^{2} θ - τ_{x y} s i n 2 θ \\ τ_{r θ} = {(σ}_{y} - σ_{x}) s i n θ c o s θ {+ τ}_{x y} c o s 2 θ \end{matrix} or \{\begin{matrix} σ_{x} = σ_{r} {c o s}^{2} θ + σ_{θ} {s i n}^{2} θ - τ_{r θ} s i n 2 θ \\ σ_{y} = σ_{r} {s i n}^{2} θ + σ_{θ} {c o s}^{2} θ + τ_{r θ} s i n 2 θ \\ τ_{x y} = {(σ}_{r} - σ_{θ}) s i n θ c o s θ {+ τ}_{r θ} c o s 2 θ \end{matrix}$

The in-plane total stress was determined by the following formula:

(33) $σ = σ_{r} + σ_{θ}$

The equivalent tensile stress or equivalent von Mises stress was defined as follows [26]:

(34) $v o n = \sqrt{\frac{{(σ_{1} - σ_{2})}^{2} + {(σ_{1})}^{2} + {(σ_{2})}^{2}}{2}}$

The principal angle could be computed according to Formula (35):

(35) $β = \frac{1}{2} a r c t a n (\frac{(2 σ_{r} - σ) s i n 2 θ + 2 τ_{r θ} c o s 2 θ}{(2 σ_{r} - σ) c o s 2 θ - 2 τ_{r θ} s i n 2 θ})$

The disk was subjected to in-plane uniform pressure that was distributed on the two diametrically opposite regions of the external circular surface (boundary conditions) [26]:

(36) $\begin{matrix} σ_{r} (R 2, θ) = \{\begin{matrix} - P, θ \in [\frac{π}{2} - α, \frac{π}{2} + α] \cup [\frac{3 π}{2} - α, \frac{3 π}{2} + α] \\ 0, i n o t h e r c a s e s \end{matrix} \\ {o t h e r w i s e, σ}_{r} (R 1, θ) = 0, τ_{r θ} (R 1, θ) = 0, τ_{r θ} (R 2, θ) = 0, θ \in (0, \frac{π}{2}] \end{matrix}$

One point that may go unnoticed is that because the step (Δx) was constant, we could change the values of the indices in the finite difference formulas. To use the finite difference method in all potentialities, in addition to having a bank of finite difference approximations, it is also necessary to manipulate the indices in the finite difference formulas, where one can:

(a). add or subtract the indices of the equations by any value, shifting the indices, which is useful at the border of the geometric domain of the mesh,
(b). use the formulas obtained for total and partial derivatives,
(c). swap the variable with respect to the derivative,
(d). change the orders of the variables of the partial derivatives, changing the indices of equations,
(e). obtain the backward equations based on the forward equations, and vice versa,
(f). with symmetry in the finite difference coefficients of different quadrants, considering the formula of quadrant 1, one can obtain that of quadrants 2, 3, and 4,
(g). observe the symmetry of the numerator coefficients of centered equations.

These characteristics of finite difference formulas are just a few examples that show the versatility and potential of their use in solving differential equation problems. Mainly in implicit formulations of the finite difference method, we arrived at linear systems that can be revolved most of the time by iterative methods or by successive and cyclic application of the equations. Table A11 in Appendix A shows the finite difference formulas utilized in solving the problem [27,28,29].

In the complex case, the problem can be mathematically modeled by Westergaard’s stress function (Z). Westergaard introduced a complex stress function, Z(z), that is related to Airy’s stress function, ϕ, by the equation:

(37) $\emptyset = R e \bar{Z} + y I m \bar{Z}$

Since Z is a complex function, it is clear that:

(38) $Z (z) = R e Z + I m Z$

where z is defined as:

(39) $z = x + y = r e^{θ}$

Note that Z was analytic over the region of interest, and the Cauchy–Riemann conditions led to:

(40) $\nabla^{2} R e Z = \nabla^{2} I m Z = 0$

This result shows that the Westergaard stress function automatically satisfied Equation (27). The bars over the stress function Z in Equation (37) indicate integration [22,30]. Thus,

(41) $\begin{matrix} \frac{d \bar{\bar{Z}}}{d z} = \bar{Z} o r \bar{\bar{Z}} = \int \bar{Z} d z \\ \frac{d \bar{Z}}{d z} = Z o r \bar{Z} = \int Z d z \\ \frac{d Z}{d z} = Z^{'} o r Z = \int Z^{'} d z \end{matrix}$

where the bar and the prime represent integration and differentiation, respectively. Substituting Equation (37) into Equation (28) yielded the Cartesian components of stress in terms of real and imaginary parts of the Westergaard stress function, as:

(42) $σ_{x x} = R e Z - y I m Z^{'} σ_{y y} = R e Z + y I m Z^{'} τ_{x y} = - y R e Z^{'}$

Equation (42) will yield stresses for functions Z(z) that are analytic; however, the stress function must satisfy boundary conditions corresponding to the problem being investigated. The Westergaard stress function, Z, can be used in fracture mechanics to model, for example, the stress field in regions adjacent to the crack tip, and Airy’s stress function, ϕ, cannot model this.

To find the complex approximations of the derivatives by finite differences, the following algorithm was used:

(1). Define the error order of the approximation and thus the number of formula intervals (Figure 1, Figure 2 and Figure 3).
(2). Define whether to use equations from the first, second, third, or fourth quadrants or central equations, thus having the discrete points of the finite difference formulas (depending on the need in the problem mesh).
(3). Assemble the complex linear system using Taylor series (Equations (3)–(6)).
(4). Solve the complex linear system and obtain approximations for the first derivative (Equations (11) and (12)).
(5). Write the first derivative approximations obtained in terms of the Shift operator E, with complex powers (Section 6).
(6). Define the derivative that must be calculated with the number of independent complex and real variables and the order of the derivative of each variable (depending on the differential equation to be solved).
(7). Perform polynomial multiplication in terms of the Shift operator E to obtain the desired finite difference approximation (Section 4 and Section 6).
(8). Perform an algebraic test of the formula obtained (Section 7).
(9). Use the finite difference formula obtained in the practical example.
(10). Apply the finite difference approximation with increasingly smaller increment steps until obtaining the desired precision and accuracy.

Table A12 and Table A13 in Appendix A show the complex finite difference formulas utilized in solving the problem.

The implicit method was used to solve the problem by finite differences, making Δr = 0.001 m and Δθ = 0.01 rad, and then comparing with Δr = 0.0001 m and Δθ = 0.001 rad using 8 decimal places. Although the edges were circular, as the steps were very small and, consequently, the elements as well, the Control–Volume Approach or Volume–Integral Approach (a technique used for an irregularly shaped boundary, where weighting coefficients are applied to account for the nonuniform spacing in the vicinity of the nonrectangular boundary) were not used [1,2,3].

The computational solution of the problem was performed, and the linear system was solved by an iterative method. Figure 5 graphically shows the results, which were compared with the analytical solutions and were successful up to the eighth decimal place of precision (ξ ≤ 10⁻⁸). To gain computational efficiency, only the nodes where x² + y² were between the internal radius R1² (0.02²) and the external radius R2² (0.04²) were activated for the calculation. Each node is a real value of the incognita (ϕ) Airy’s stress function or a complex value of the incognita (Z) Westergaard’s stress function, but it is also a finite difference equation. At the edge of the boundary of the inner radius and outer radius, an increase in the grid of eight radial lines was provided for the stress boundary conditions. Thus, the number of unknowns (nodes) was equal to the number of equations in a large linear system that can be solved by Gauss–Seidel-type iterative techniques [2,31]. For this example of compression rings, the exact value of the analytical solution of the differential equation at each position, x and y, was also calculated. Thus, the error of using the finite difference method was obtained in the real case and in the complex case.

The results obtained by finite differences were compared with the analytical solution up to an error of less than eight decimal places. This example demonstrated the effectiveness of the finite difference method in accurately solving an engineering application involving fourth-order differential equations and circular geometry [26].

Both in the real model with the Airy function and in the complex model with the Westergaard function, the finite difference method was successful in obtaining the stress diagrams. No significant visual differences were observed between the real solution and the complex solution. Table 3 shows a comparison of the errors of the real solution in relation to the errors of the complex solution [23]. It should be noted that the complex solution obtained fewer errors, being considered much better than the real solution.

To date, we do not have a full understanding of why the complex solution is more accurate than the real solution. However, it is believed that the order of error in complex finite difference formulas is greater than that used in real finite difference formulas. Another reason comes from the modulus of the coefficients of the terms in complex finite difference equations being smaller than the real coefficients, reducing rounding errors in mathematical operations.

Finite difference formulas in the complex plane are in general more accurate than finite difference formulas in the real domain [17,18,19,20]. The reasons for this are that complex derivatives depend equally on data from all directions around the point, i.e., stencils must not extend very far, only along a line, and the Cauchy–Riemann equations impose restrictions on the range of functions that need to be considered. Even as stencil sizes and orders of accuracy increase, complex planar finite difference stencils continue to extract their key information from a very small neighborhood of the point of interest. By using standard finite difference approximations along the real axis and increasing the stencil width and precision order to a fixed derivative, the weights converge to a limit, known as the pseudospectral method. This contrasts with traditional pseudospectral methods which, with their slow algebraic decay of weights, rely on distant data, even when approximating local operators as derivatives. It seems that the derivative at a point cannot depend on the direction of approach in the complex plane, resulting in greater stability for the finite difference method. Holding a stencil size fixed and increasing the order p of the derivative $d^{p} / {d z}^{p}$ (as far as the stencil size permits) will typically make the weights (following the leading factor $1 / h^{p}$ ) increase rapidly. The Cauchy–Riemann equations also seemed to smooth the variations of finite differences in the complex domain, which improved the numerical precision and accuracy of numerical solutions of differential equations. The Euler–Maclaurin approximations benefited from a compensatory decrease in the magnitude of the coefficients for successive terms. If a very high-order derivative is to be approximated, algorithms based on Cauchy integrals may be preferable to finite difference approximations. Another option could be to find sets of finite difference weights where some orders of precision have been traded off against reductions in the magnitude of the weights. This has been used successfully in contexts such as numerical quadrature and in cases of Gregory-type and Newton–Cotes-type formulas, respectively.

Another issue that helps in the use of complex finite difference approximations is the square increase in the points used per interval. Thus, a spacing interval uses two points in the real case and four points in the complex case, two intervals uses three points in the real case and nine points in the complex case, three intervals uses four points in the real case and sixteen points in the complex case, and so on. This increases the precision and error order of the finite difference approximation and decreases the absolute value of the coefficient values of the finite difference approximations.

Table 4 shows a comparison between the various finite difference formulas used in numerically solving the stress problem. Note that only the main finite difference approximations for the fourth-order derivatives have been changed, as the initial and boundary conditions that use Dirichlet, Neumann, or second- and third-order derivatives are the same finite difference equations.

Comparing finite difference approximations in the complex field, it was noted that the central equations and more points had advantages, as well as equations that have geometric symmetry between the four quadrants of the complex plane and that try to get closer to circles centered at the origin. In other words, the symmetrical and regular distribution of the stencils helps to obtain more precise and accurate finite difference approximations.

Many other applications and tests of finite difference derivative approximations are provided in the literature [21,32,33,34,35,36,37,38,39,40,41].

10. Conclusions

The finite difference method is an important numerical method for solving both differential equations and computer simulations. For the same order of error, approximations of derivatives by finite differences were obtained with a smaller absolute value of the coefficients in the complex domain than in the real domain, which reduced the propagation of the rounding error. Finite difference formulas in the complex plane are in general more precise and accurate than finite difference formulas in the real domain. Furthermore, a much greater diversity of finite difference equations was obtained in the complex field than in the real domain, in all orders of errors and for the most varied boundary and initial conditions. The work is innovative, with a complex potentiation of numerical operators.

It was possible in the complex domain to obtain finite difference approximations of a high order of error or precision with few points of the function, even reaching equations with an error order of O(hⁿ) and with only n points of the function. The reduced number of points in the finite difference approximation was important, especially in the boundary and initial conditions. Obtaining approximations of derivatives with a high order of error improved the precision of the numerical resolution of differential equations and reduced the computational effort in the finite difference method. The research was successful in obtaining approximations of derivatives by finite differences in the complex domain with a high order of error and with few points in the equations, which increased the applicability of the finite difference method and reduced restrictions on its use.

Due to the complex boundary and initial conditions that appear in structural engineering and solid mechanics problems, a large number of finite difference derivative approximations are necessary to solve the problems with accuracy. This work showed that, with modern microcomputers and a database with several finite difference equations for derivative approximations, many structural analysis problems can be solved with high numerical precision by applying the finite difference method.

This paper showed that although it is outdated, the finite difference method is still very important and cannot be abandoned. It is believed that it is best to use the finite difference method in conjunction with other numerical methods of solving differential equations, such as the finite element method and the boundary element method. The finite differences method presents high stability and easy application in many practical cases of problems involving differential equations, mainly problems involving long time series or sequences and derived from high orders or problems with part of the geometry having a simpler rectangular shape and large areas and volumes with homogeneous media. In these cases, the finite difference method has great advantages.

This article showed how structural mechanics problems can be solved using the finite difference method in the complex domain. Furthermore, it showed how advantageous it is to use complex mathematical models in relation to models in the real domain. Westergaard stress functions are more general and comprehensive than Airy’s stress function, and the finite difference method, which is a numerical method widely used in computer simulations, seemed to work better in the complex plane.

This paper showed how to obtain, in a simple and practical way, finite difference equations that perform approximations of the calculation of derivatives with high orders of error. Correctly using these finite difference equations, we obtained precise and exact solutions of differential equations in practical science and engineering problems. Thus, the finite difference method remains a powerful and robust numerical method for solving differential equation problems that efficiently utilizes recent computational advances in hardware and software.

Since finite differences formulas in the complex plane are generally considerably more accurate than traditional ones that use function values only along the real axis, it is concluded that it may be advantageous to solve differential equation problems using the finite difference method with variables in the complex domain instead of the real domain, which can even generalize many mathematical models of physical and engineering phenomena. The application of the finite difference method in the complex plane can be used to determine whether a variable is better modeled as belonging to the real domain or the complex domain in the study or analysis of a certain phenomenon.

As Supplementary Material of the paper, many of the approximations of derivatives by complex finite differences are presented in a PDF (https://zenodo.org/records/10901007, accessed on 3 April 2024).

Author Contributions

All authors reviewed the manuscript. All authors wrote the main manuscript text. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

Data is contained within the article or Supplementary Material.

Acknowledgments

The authors thank Pontifícia Universidade Católica de Minas Gerais-PUC Minas “Pontifical Catholic University of Minas Gerais” for the support. The authors thank the financial support of the Conselho Nacional de Desenvolvimento Científico e Tecnológico-CNPq “National Counsel of technological and scientific Development” and the Fundação de Amparo à Pesquisa de Minas Gerais-FAPEMIG “Foundation for Research Support of Minas Gerais”. This study was financed in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior—Brasil (CAPES)—Finance Code 001.

Conflicts of Interest

The authors declare no conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Spacing of the independent variable in the complex plane with intervals 01 and 02.

Figure 2. Spacing of the independent variable in the complex plane with interval 03.

Figure 3. Spacing of the independent variable in the complex plane with interval 04.

Figure 4. Ring in compression by a constant force with an angle of 2α.

Figure 5. Results obtained by finite differences.

Table 1

Definition of numeric operators.

Operator	Symbol	Definition
Shift	E	Eu(x_i) = u(x_i+1) or Ey_i = y_i+1
Forward difference	Δ	Δu(x_i) = u(x_i+1) − u(x_i) or Δy_i = y_i+1 − y_i or Δ = E − 1
Backward difference	∇	∇u(x_i) = u(x_i) − u(x_i−1) or ∇y_i = y_i − y_i−1 or ∇ = 1 − E⁻¹
Central difference	δ	δu(x_i) = u(x_i+½) − u(x_i−½) or δy_i = y_i+½ − y_i−½ or δ = E^1/2 − E^−1/2
Average	μ	μu(x_i) = (½)[ u(x_i+½) + u(x_i−½) ] or μy_i = (½)(y_i+½ + y_i−½ ) or μ = (½)(E^½ + E^−½ )
Average central difference	μδ	μδu(x_i) = (½)[ u(x_i+1) − u(x_i−1) ] or μδy_i = (½)(y_i+1 − y_i−1) or μδ = (½)(E¹ − E⁻¹)

Table 2

Successive application of numeric operators (potentiation).

Operator	Symbol	Successive Application
Shift	E	Eⁿu(x_i) = u(x_i+n) or Eⁿy_i = y_i+n or Eⁿ = E(E(E(….E)))
Forward difference	Δ	Δⁿu(x_i) = Δⁿ⁻¹u(x_i+1) − Δⁿ⁻¹u(x_i) or Δⁿy_i = Δⁿ⁻¹y_i+1 − Δⁿ⁻¹y_i or Δⁿ = Δ(Δ(Δ(….Δ)))
Backward difference	∇	∇ⁿu(x_i) = ∇ⁿ⁻¹u(x_i) − ∇ⁿ⁻¹u(x_i−1) or ∇ⁿy_i = ∇ⁿ⁻¹y_i − ∇ⁿ⁻¹y_i−1 or ∇ⁿ = ∇(∇(∇(….∇)))
Central difference	δ	δⁿu(x_i) = δⁿ⁻¹u(x_i+½) − δⁿ⁻¹u(x_i−½) or δⁿy_i = δⁿ⁻¹y_i+½ − δⁿ⁻¹y_i−½ or δ ⁿ = δ(δ(δ(….δ))) δ²u(x_i) = δu(x_i+½) − δu(x_i−½) = [u(x_i+½+½) − u(x_i+½−½)] − [u(x_i−½+½) − u(x_i−½−½)] = [u(x_i+1) − u(x_i)] − [u(x_i) − u(x_i−1)] = u(x_i+1) − u(x_i) − u(x_i) + u(x_i−1) = u(x_i+1) − 2 u(x_i) + u(x_i−1) δ²u(xi) = u(x_i+1) − 2 u(x_i) + u(x_i−1) or δ²y_i = y_i+1 − 2 y_i + y_i−1 δ⁴u(xi) = δ²u(x_i+1) − 2 δ²u(x_i) + δ²u(x_i−1) or δ²y_i = δ²y_i+1 − 2 δ²y_i + δ²y_i−1 δⁿu(xi) = δⁿ⁻²u(x_i+1) − 2δⁿ⁻²u(x_i) + δⁿ⁻²u(x_i−1) or δⁿy_i = δⁿ⁻²y_i+1 − 2δⁿ⁻²y_i + δⁿ⁻²y_i−1
Average central difference	μδ	μδⁿu(x_i) = (½) [δⁿ⁻¹u(x_i+½) − δⁿ⁻¹u(x_i−½)] or μδⁿy_i = (½) (δⁿ⁻¹y_i+½ − δⁿ⁻¹y_i−½ ) or μδⁿ = μδ(δ(δ(….δ))) μδu(x_i) = (½) [u(x_i+1) − u(x_i−1) ] or μδy_i = (½)(y_i+1 − y_i−1) μδ³u(x_i) = (½) [δ²u(x_i+1) − δ²u(x_i−1)] or μδ³y_i = (½)(δ²y_i+1 − δ²y_i−1) μδⁿu(x_i) = (½) [ δⁿ⁻¹u(x_i+1) − δⁿ⁻¹u(x_i−1)] or μδⁿy_i = (½)(δⁿ⁻¹y_i+1 − δⁿ⁻¹y_i−1)

Table 3

Comparison of the real solution with the complex solution of the Theory of Elasticity problem of circular rings in compression.

Step or Increment of Cartesian Coordinates x and y—Δx or Δy (m)	Step or Increment of Polar Coordinates Radial Distance—Δr (m)	Step or Increment of Polar Coordinates Polar Angle—Δθ (rad)	Average Stresses Error between the Exact Value and the Real Solution Due to Finite Differences at the Mesh Nodes—σ or τ (Pa)	Average Stresses Error between the Exact Value and the Complex Solution Due to Finite Differences at the Mesh Nodes—σ or τ (Pa)
h = 0.1	0.14142	0.07854	0.02805	0.006871
h = 0.01	0.014142	0.007854	2.4463 × 10⁻⁶	5.8212 × 10⁻¹¹
h = 0.001	0.0014142	0.0007854	7.7650 × 10⁻¹⁰	8.5163 × 10⁻¹⁹
h = 10⁻⁴	1.14142 × 10⁻⁴	7.854 × 10⁻⁵	1.5339 × 10⁻¹⁴	1.2798 × 10⁻²⁷
h = 10⁻⁵	1.14142 × 10⁻⁵	7.854 × 10⁻⁶	7.4519 × 10⁻¹⁸	6.7628 × 10⁻³⁵
h = 10⁻⁶	1.14142 × 10⁻⁶	7.854 × 10⁻⁷	1.7867 × 10⁻²²	8.6599 × 10⁻⁴³

Table 4

Comparison of the precision of numerical resolution of differential equations between the different grids in the complex plane.

Grid or Set of Points of the Main Stencil Used in the Numerical Resolution of the Stress Calculation Problem:	Average Stresses Error between the Exact Value and the Complex Solution Due to Finite Differences at the Mesh Nodes with Step of Cartesian Coordinates x and y Equal to:
	h = 0.1	h = 0.01	h = 0.001	h = 10⁻⁴	h = 10⁻⁵
Square with Interval 01—Quadrant 01—Table A4a	2.04 × 10²	6.83 × 10⁻¹	9.59 × 10⁻⁴	4.50 × 10⁻⁷	9.74 × 10⁻¹⁰
Square with Interval 01—Quadrant 02—Table A5a	3.88 × 10²	8.91 × 10⁻¹	7.38 × 10⁻⁴	1.75 × 10⁻⁸	2.48 × 10⁻¹⁰
Square with Interval 01—Quadrant 03—Table A6a	3.88 × 10²	3.36 × 10⁻¹	6.95 × 10⁻⁴	4.28 × 10⁻⁷	1.29 × 10⁻¹⁰
Square with Interval 01—Quadrant 04—Table A7a	2.27 × 10	1.26 × 10⁻¹	1.40 × 10⁻⁴	4.74 × 10⁻⁷	1.38 × 10⁻¹⁰
Square with Interval 02—Quadrant 01—Table A4b	3.81 × 10⁻³	3.37 × 10⁻¹¹	4.86 × 10⁻²⁰	5.67 × 10⁻²⁷	7.91 × 10⁻³⁵
Square with Interval 02—Quadrant 02—Table A5b	3.17 × 10⁻³	7.21 × 10⁻¹¹	6.76 × 10⁻¹⁹	7.29 × 10⁻²⁷	4.08 × 10⁻³⁵
Square with Interval 02—Quadrant 03—Table A6b	8.40 × 10⁻³	3.19 × 10⁻¹¹	6.46 × 10⁻¹⁹	5.17 × 10⁻²⁷	2.78 × 10⁻³⁵
Square with Interval 02—Quadrant 04—Table A7b	7.78 × 10⁻³	6.94 × 10⁻¹¹	4.22 × 10⁻¹⁹	2.95 × 10⁻²⁷	7.70 × 10⁻³⁵
Square with Interval 02—Central—Table A8a	3.74 × 10⁻³	4.45 × 10⁻¹¹	3.00 × 10⁻¹⁹	3.47 × 10⁻²⁷	1.45 × 10⁻³⁵
Square with Interval 03—Quadrant 01—Table A4c	5.93× 10⁻¹⁰	4.49 × 10⁻²⁵	7.82 × 10⁻⁴⁰	4.18 × 10⁻⁵⁵	3.33 × 10⁻⁷¹
Square with Interval 03—Quadrant 02—Table A5c	9.28 × 10⁻¹⁰	4.63 × 10⁻²⁵	2.67 × 10⁻⁴⁰	6.06 × 10⁻⁵⁵	8.83 × 10⁻⁷⁰
Square with Interval 03—Quadrant 03—Table A6c	8.08× 10⁻¹⁰	4.49 × 10⁻²⁵	1.47 × 10⁻⁴⁰	6.56 × 10⁻⁵⁷	3.80 × 10⁻⁷⁰
Square with Interval 03—Quadrant 04—Table A7c	9.42 × 10⁻¹⁰	5.81 × 10⁻²⁵	2.63 × 10⁻⁴⁰	7.71 × 10⁻⁵⁵	2.56 × 10⁻⁷⁰
Square with Interval 04—Quadrant 01—Table A4d	2.85 × 10⁻¹⁹	4.90 × 10⁻⁴³	3.28 × 10⁻⁶⁷	3.87 × 10⁻⁹³	9.37 × 10⁻¹¹⁵
Square with Interval 04—Quadrant 02—Table A5d	2.55 × 10⁻¹⁹	2.25 × 10⁻⁴³	6.99 × 10⁻⁶⁹	4.82 × 10⁻⁹¹	7.77 × 10⁻¹¹⁵
Square with Interval 04—Quadrant 03—Table A6d	5.51 × 10⁻¹⁹	5.93 × 10⁻⁴³	2.88 × 10⁻⁶⁷	6.04 × 10⁻⁹¹	3.07 × 10⁻¹¹⁵
Square with Interval 04—Quadrant 04—Table A7d	4.76 × 10⁻¹⁹	2.18 × 10⁻⁴³	5.04 × 10⁻⁶⁷	9.71 × 10⁻⁹²	3.72 × 10⁻¹¹⁵
Square with Interval 04—Central—Table A8b	4.61 × 10⁻¹⁹	1.73 × 10⁻⁴³	2.84 × 10⁻⁶⁷	6.20 × 10⁻⁹²	1.63 × 10⁻¹¹⁶
Square with Interval 06—Central—Table A8c	3.24 × 10⁻⁴³	8.05 × 10⁻⁹²	2.40 × 10⁻¹⁴⁰	2.89 × 10⁻¹⁸⁷	4.55 × 10⁻²³⁵
Triangle 01—Quadrant 01—Figure A1a	2.72 × 10³	2.49 × 10	5.95 × 10⁻¹	4.75 × 10⁻⁴	2.20 × 10⁻⁵
Triangle 01—Quadrant 02—Figure A1b	5.88 × 10³	9.82 × 10	1.92 × 10⁻¹	7.83 × 10⁻³	6.13 × 10⁻⁵
Triangle 01—Quadrant 03—Figure A1c	5.08 × 10³	1.60 × 10	9.20 × 10⁻¹	5.90 × 10⁻³	6.71 × 10⁻⁵
Triangle 01—Quadrant 04—Figure A1d	8.42 × 10³	5.53 × 10	1.41 × 10⁻¹	6.50 × 10⁻³	2.90 × 10⁻⁵
Diamond 01—Central—Figure A1e	2.56 × 10⁰	8.22 × 10⁻⁴	4.12 × 10⁻⁷	2.37 × 10⁻¹¹	2.71 × 10⁻¹⁵
Triangle 02—Quadrant 01—Figure A2a	9.28 × 10⁰	2.52 × 10⁻⁵	9.63 × 10⁻¹⁰	3.15 × 10⁻¹⁵	8.28 × 10⁻²⁰
Triangle 02—Quadrant 02—Figure A2b	8.34 × 10⁰	1.43 × 10⁻⁵	1.60 × 10⁻¹⁰	6.00 × 10⁻¹⁵	8.71 × 10⁻²⁰
Triangle 02—Quadrant 03—Figure A2c	5.20 × 10⁰	3.53 × 10⁰⁵	2.10 × 10⁻¹⁰	4.60 × 10⁻¹⁵	2.78 × 10⁻²¹
Triangle 02—Quadrant 04—Figure A3a	1.81 × 10⁰	6.62 × 10⁻⁵	7.90 × 10⁻¹⁰	3.78 × 10⁻¹⁵	1.78 × 10⁻²⁰
Diamond 02—Central—Figure A3b	2.09 × 10⁻⁶	2.57 × 10⁻¹⁷	3.88 × 10⁻²⁸	1.10 × 10⁻³⁹	2.79 × 10⁻⁵⁰
Hexagon—Central—Figure A4	2.96 × 10⁻¹²	2.63 × 10⁻²⁷	1.18 × 10⁻⁴³	1.59 × 10⁻⁵⁹	2.87 × 10⁻⁷⁵
Octagon—Central—Figure A5	2.00 × 10⁻²³	1.50 × 10⁻⁵¹	1.59 × 10⁻⁷⁹	2.88 × 10⁻¹⁰⁷	3.03 × 10⁻¹³⁵
Another Octagon with horizontal and vertical sides—Central—Figure A6	1.22 × 10⁻³²	7.50 × 10⁻⁶⁸	1.90 × 10⁻¹⁰³	2.38 × 10⁻¹³⁹	8.11 × 10⁻¹⁷⁶
Dodecagon with horizontal and vertical sides—Central—Figure A7 and Figure A8	1.32 × 10⁻⁵⁵	1.51 × 10⁻¹¹⁵	1.10 × 10⁻¹⁷⁵	7.71 × 10⁻²³⁶	1.73 × 10⁻²⁹⁶

Supplementary Materials

The supporting information can be downloaded at: https://zenodo.org/records/10901007, accessed on 3 April 2024.

References

1. Hoffman, J.D.; Hoffman, J.D.; Frankel, S. Numerical Methods for Engineers and Scientists; CRC Press: Boca Raton, FL, USA, 2018; [DOI: https://dx.doi.org/10.1201/9781315274508]

2. Steven, C.C. Applied Numercial Methods with MATLAB for Engineers and Scientist; McGraw Hill: New York, NY, USA, 2013; Volume 53.

3. Chapra, S.C.; Canale, R.P. Numerical Methods for Engineers; 2nd ed. New Age International: New Delhi, India, 2006.

4. Nayak, G.C. Applied Numerical Methods, B. Carnahan, H.A. Lither and J. O. Wilkes, Wiley, New York, 1969. No. of Pages: 604. Price: £6·60. Int. J. Numer. Methods Eng.; 1972; 4, pp. 599-600. [DOI: https://dx.doi.org/10.1002/nme.1620040415]

5. Seiler, M.C.; Seiler, F.A. Numerical Recipes in C: The Art of Scientific Computing. Risk Anal.; 1989; 9, pp. 415-416. [DOI: https://dx.doi.org/10.1111/j.1539-6924.1989.tb01007.x]

6. Dukkipati, R.V. Numerical Methods; Courier Corporation: New Delhi, India, 2010.

7. Chapra, S.C. Métodos Numéricos Para Ingenieros; McGraw Hill: New York, NY, USA, 2011.

8. Chen, Q.; Wang, P.; Zhu, D. A Second-Order Finite-Difference Method for Derivative-Free Optimization. J. Math.; 2024; 2024, 1947996. [DOI: https://dx.doi.org/10.1155/2024/1947996]

9. Chelnokov, Y.N. Quaternion Methods and Models of Regular Celestial Mechanics and Astrodynamics. Appl. Math. Mech.; 2022; 43, pp. 21-80. [DOI: https://dx.doi.org/10.1007/s10483-021-2797-9]

10. Roberts, J.L. A Method for Calculating Meshless Finite Difference Weights. Int. J. Numer. Methods Eng.; 2008; 74, [DOI: https://dx.doi.org/10.1002/nme.2169]

11. Khan, I.R.; Ohba, R. Closed-Form Expressions for the Finite Difference Approximations of First and Higher Derivatives Based on Taylor Series. J. Comput. Appl. Math.; 1999; 107, pp. 179-193. Corrigendum in J. Comput. Appl. Math. 2001, 130, 385 [DOI: https://dx.doi.org/10.1016/S0377-0427(99)00088-6] [DOI: https://dx.doi.org/10.1016/s0377-0427(00)00263-6]

12. Khan, I.R.; Ohba, R.; Hozumi, N. Mathematical Proof of Closed Form Expressions for Finite Difference Approximations Based on Taylor Series. J. Comput. Appl. Math.; 2003; 150, pp. 303-309. [DOI: https://dx.doi.org/10.1016/S0377-0427(02)00667-2]

13. Khan, I.R.; Ohba, R. Taylor Series Based Finite Difference Approximations of Higher-Degree Derivatives. J. Comput. Appl. Math.; 2003; 154, pp. 115-124. [DOI: https://dx.doi.org/10.1016/S0377-0427(02)00816-6]

14. Khan, I.R.; Ohba, R. New Finite Difference Formulas for Numerical Differentiation. J. Comput. Appl. Math.; 2000; 126, pp. 269-276. [DOI: https://dx.doi.org/10.1016/S0377-0427(99)00358-1]

15. Li, J. General Explicit Difference Formulas for Numerical Differentiation. J. Comput. Appl. Math.; 2005; 183, pp. 29-52. [DOI: https://dx.doi.org/10.1016/j.cam.2004.12.026]

16. Feagin, T. High-Order Explicit Runge-Kutta Methods Using m-Symmetry. Neural Parallel Sci. Comput.; 2012; 20, pp. 437-458.

17. Fornberg, B. Infinite-Order Accuracy Limit of Finite Difference Formulas in the Complex Plane. IMA J. Numer. Anal.; 2023; 43, pp. 3055-3072. [DOI: https://dx.doi.org/10.1093/imanum/drac064]

18. Fornberg, B.; Piret, C. Computation of Fractional Derivatives of Analytic Functions. J. Sci. Comput.; 2023; 96, 79. [DOI: https://dx.doi.org/10.1007/s10915-023-02293-4]

19. Fornberg, B. Finite Difference Formulas in the Complex Plane. Numer. Algorithms; 2022; 90, pp. 1305-1326. [DOI: https://dx.doi.org/10.1007/s11075-021-01231-5]

20. Abrahamsen, D.; Fornberg, B. Solving the Korteweg-de Vries Equation with Hermite-Based Finite Differences. Appl. Math. Comput.; 2021; 401, 126101. [DOI: https://dx.doi.org/10.1016/j.amc.2021.126101]

21. Li, T.; Wang, Q.W. Structure Preserving Quaternion Full Orthogonalization Method with Applications. Numer. Linear Algebra Appl.; 2023; 30, e2495. [DOI: https://dx.doi.org/10.1002/nla.2495]

22. SEM. Handbook on Experimental Mechanics. Exp. Tech.; 1989; 13, 21. [DOI: https://dx.doi.org/10.1111/j.1747-1567.1989.tb01038.x]

23. Tokovyy, Y.V.; Hung, K.M.; Ma, C.C. Determination of Stresses and Displacements in a Thin Annular Disk Subjected to Diametral Compression. J. Math. Sci.; 2010; 165, pp. 342-354. [DOI: https://dx.doi.org/10.1007/s10958-010-9803-6]

24. Markides, C.F.; Kourkoulis, S.K. Stresses and Displacements in an Elliptically Perforated Circular Disc under Radial Pressure. Eng. Trans.; 2014; 62, pp. 131-169.

25. Vasil’ev, V.V.; Fedorov, L.V. Stress functions in elasticity theory. Mech. Solids; 2022; 57, pp. 770-778. [DOI: https://dx.doi.org/10.3103/S0025654422040197]

26. Sadd, M.H. Elasticity: Theory, Applications, and Numerics; Elsevier: Amsterdam, The Netherlands, 2020; [DOI: https://dx.doi.org/10.1016/C2017-0-03720-5]

27. Jiang, Y.; Pu, C. The Problem of I–II Combined Plane Crack Solved with Westergaard Stress Function. Mech. Eng.; 2020; 42, pp. 504-507. [DOI: https://dx.doi.org/10.6052/1000-0879-19-476]

28. Dumont, N.A.; Mamani, E.Y. Generalized Westergaard Stress Functions as Fundamental Solutions. CMES Comput. Model. Eng. Sci.; 2011; 78, pp. 109-149.

29. Sun, C.T.; Farris, T.N. On the Completeness of the Westergaard Stress Functions. Int. J. Fract.; 1989; 40, pp. 73-77. [DOI: https://dx.doi.org/10.1007/BF01150867]

30. Leven, M.M. Discussion: “A New Method to ‘Lock-In’ Elastic Effects for Experimental Stress Analysis” (Dally, J.W., Durelli, A.J., and Riley, W.F., 1958, ASME J. Appl. Mech., 25, Pp. 189–195). J. Appl. Mech.; 1959; 26, 152. [DOI: https://dx.doi.org/10.1115/1.4011951]

31. Clough, D.E.; Chapra, S.C. Introduction to Engineering and Scientific Computing with Python; Informa UK Limited: London, UK, 2022; [DOI: https://dx.doi.org/10.1201/9781003256861]

32. Magalhães, P.A.A. New Equations for Phase Evaluation in Measurements with an Arbitrary but Constant Phase Shift between Captured Intensity Signs. Opt. Eng.; 2009; 48, 113602. [DOI: https://dx.doi.org/10.1117/1.3265438]

33. Cunha, P.P.S.R.; da Souza, P.M.; Valente, L.C.; Mello, G.M.V.; de Junior, P.A.A.M. Analysis of Induced Drag and Vortex at the Wing Tip of a Blended Wing Body Aircraft. Int. J. Adv. Eng. Res. Sci.; 2018; 5, pp. 7-9. [DOI: https://dx.doi.org/10.22161/ijaers.5.6.2]

34. Magalhaes, P.A.A.; Magalhaes, C.A. Higher-Order Newton-Cotes Formulas. J. Math. Stat.; 2010; 6, pp. 193-204. [DOI: https://dx.doi.org/10.3844/jmssp.2010.193.204]

35. Almeida Magalhaes, C.; Americo Almeida Magalhaes, P. New Numerical Methods for the Photoelastic Technique with High Accuracy. J. Appl. Phys.; 2012; 112, 083111. [DOI: https://dx.doi.org/10.1063/1.4761979]

36. Lamon, G.P.S.; Magalhães Júnior, P.A.A.; Carneiro, J.R.G. Cole type Pitot tube, discharge factor survey and calibration. Meas. Sens.; 2024; 33, 101152. [DOI: https://dx.doi.org/10.1016/j.measen.2024.101152]

37. López-Pérez, A.; Febrero-Bande, M.; González-Manteiga, W. Parametric estimation of diffusion processes: A review and comparative study. Mathematics; 2021; 9, 859. [DOI: https://dx.doi.org/10.3390/math9080859]

38. Berdyshev, A.; Baigereyev, D.; Boranbek, K. Numerical Method for Fractional-Order Generalization of the Stochastic Stokes–Darcy Model. Mathematics; 2023; 11, 3763. [DOI: https://dx.doi.org/10.3390/math11173763]

39. Pekmen Geridonmez, B.P.; Oztop, H.F. Entropy Generation Due to Magneto-Convection of a Hybrid Nanofluid in the Presence of a Wavy Conducting Wall. Mathematics; 2022; 10, 4663. [DOI: https://dx.doi.org/10.3390/math10244663]

40. Sun, J.; Wang, L.; Gong, D. A Joint Optimization Algorithm Based on the Optimal Shape Parameter–Gaussian Radial Basis Function Surrogate Model and Its Application. Mathematics; 2023; 11, 3169. [DOI: https://dx.doi.org/10.3390/math11143169]

41. Alahmadi, R.A.; Raza, J.; Mushtaq, T.; Abdelmohsen, S.A.M.; RGorji, M.; Hassan, A.M. Optimization of MHD Flow of Radiative Micropolar Nanofluid in a Channel by RSM: Sensitivity Analysis. Mathematics; 2023; 11, 939. [DOI: https://dx.doi.org/10.3390/math11040939]

© 2024 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/). Notwithstanding the ProQuest Terms and Conditions, you may use this content in accordance with the terms of the License.

Numerical Resolution of Differential Equations Using the Finite Difference Method in the Real and Complex Domain

Content area

Abstract

Full text