Computational Advances in the Implementation of the Finite Difference Method Applied to

Solving Problems in the Theory of Elasticity Using Airys Stress Functions

Ana Laura Mendonca Almeida Magalhães1, Pedro Américo Almeida Magalhães Junior*, Pedro
Henrique Mendonça Almeida Magalhães2, Cristina Almeida Magalhães3, Geraldo Pedro da Silva
Lamon1

1
Pontifícia Universidade Católica de Minas Gerais (PUCMINAS), Av: Dom José Gaspar, 500 –
Prédio 10 - Coração Eucarístico - Belo Horizonte – MG – Cep 30535-901, Brazil
*
Pontifícia Universidade Católica de Minas Gerais (PUCMINAS), Av. Brasil, 2023 - Funcionários,
Belo Horizonte - MG, Cep 30140-002, Brazil - Telephone: +55 (31) 99951-3121 – E-mail:
pamerico@[Link] (Corresponding Author)
2
Universidade Federal de Minas Gerais (UFMG), Av. Pres. Antônio Carlos, 6627 - Pampulha,
Belo Horizonte – MG - Cep 31270-901, Brazil
3
Centro Federal de Educação Tecnológica de Minas Gerais (Cefet-MG), Av. Amazonas 7675-
Nova Gameleira - Belo Horizonte - MG – Cep 30510-000, Brazil

Abstract

The paper exposes how recent computational advances in hardware and software favor and
facilitate the use of the finite difference method, especially in more complex practical problems
involving differential equations. Currently, complex problems in the Theory of Elasticity involving
Airys Stress Function with fourth-order partial derivatives, Cartesian and polar coordinates and
exotics and different types of initial and boundary conditions can be solved without
programming using only electronic spreadsheets. The finite difference method is an important
numerical method for solving both ordinary differential equations and partial differential
equations. However, in recent years, this method has been losing ground to other numerical
methods, such as the finite element method and boundary element method. One of the reasons
for this trend is the difficulty of obtaining approximations of derivatives by finite differences in
high precision and ease of use in practical problems. The paper presents a way to deduce finite
difference formulas that approximate the derivatives, mainly high-order and very precise
formulas using numerical operators. As supplementary material of the paper, it follows many of
the approximations of derivatives by finite differences in PDF and text format, following the
table model proposed in the paper to store the finite difference equations. Thus, the article
provides a large form of finite difference equations for approximating derivatives, which can be
used in various applications.

Keywords: Airy stress function; Theory of Elasticity; finite difference methods; error order;
numerical resolution using spreadsheets; computational advances for numerical methods;
simplification of the use of numerical techniques
 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–7]. 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:

𝑑𝑥 ≈ Δ𝑥, 𝑑𝑦 ≈ Δ𝑦

𝑑𝑦 Δ𝑦 𝑑2 𝑦 Δ2 𝑦 𝑑3 𝑦 Δ3 𝑦
≈ , ≈ , ≈ (1)
𝑑𝑥 Δ𝑥 𝑑𝑥 2 Δ𝑥 2 𝑑𝑥 3 Δ𝑥 3

∂𝑢 Δ𝑢 ∂2 𝑢 Δ2 𝑢 ∂3 𝑢 Δ3 𝑢 ∂2 𝑢 Δ2 𝑢
{∂𝑥 ≈ Δ𝑥 , ∂𝑥 2
≈ Δ𝑥 2 , ∂𝑥 3
≈ Δ𝑥3 , ∂𝑥 ∂𝑦
≈ Δ𝑥Δ𝑦

The error is the difference in the absolute value of the exact value and of the approximate value
(Error = |Valueexact – Valueapprox|). In principle, shortening the calculation step reduces the
number of errors. The ratio of step to error decrease is the order of error (e), indicated by O(he),
where h = x = xi+1-xi is the constant step of the independent variables[1,8]. Thus, if a numerical
method is of error order O(h4), it means that if you shorten the step by half, the error decreases
(1/2)4= 1/16, that is, the error decreases approximately 16 times. If you reduce the step by 10
times, the error drops 10 4 = 10000 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[9].

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 has a 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[10]. Numeric operators have as 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. A discrete function (tabulated function) will be represented by
the set of points {(x1, y1), (x2, y2), …, (xn, yn)} or {[x1, f(x1)], [x2, f(x2)], …, [xn, f(xn)]}, where yn= f(xn)
and h=x=xi+1-xi is a constant step of the independent variables (h=x= x2-x1= x3-x2 = x4-x3 = xi+1-
xi = … = xn-xn-1). Table 1 defines the main numeric operators:

Table 1: Definition of numeric operators.

Operator Symbol Definition

Shift E Ef(xi) = f(xi+1) or Eyi = yi+1
Forward  f(xi) = f(xi+1) – f(xi) or yi = yi+1 – yi or  = E - 1
difference
Backward  f(xi) = f(xi) – f(xi-1) or yi = yi – yi-1 or  = 1 – E-1
Difference
Central  f(xi) = f(xi+½) – f(xi-½) or yi = yi+½ – yi-½ or  = E1/2 – E-1/2
Difference
Average  f(xi) = (½)[ f(xi+½) + f(xi-½) ] or yi = (½)(yi+½ + yi-½ ) or
 = (½)(E½ + E-½ )
Average Central  f(xi) = (½)[ f(xi+1) - f(xi-1) ] or yi = (½)(yi+1 - yi-1) or
Difference  = (½)(E1 - E-1)

In addition, the differential operator is usually denoted by D (the differential operator (D) is not
𝑑 𝑑2
a numeric operator), where 𝐷𝑓(𝑥) = 𝑑𝑥 𝑓(𝑥) = 𝑓´(𝑥) 𝑜𝑟 𝐷 2 𝑓(𝑥) = 𝑑𝑥 2 𝑓(𝑥) =
𝑑3 𝑑4
𝑓´´(𝑥) 𝑜𝑟 𝐷 3 𝑓(𝑥) = 𝑓(𝑥) = 𝑓´´´(𝑥) 𝑜𝑟 𝐷 4 𝑓(𝑥) = 4 𝑓(𝑥) =
𝑑𝑥 3 𝑑𝑥
𝑑𝑛
𝑓 (4) (𝑥) 𝑜𝑟 … 𝑜𝑟 𝐷 𝑛 𝑓(𝑥) = 𝑑𝑥 𝑛 𝑓(𝑥) = 𝑓 (𝑛) (𝑥). 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:
Table 2: Successive application of numeric operators.

Operator Symbol Successive Application

Shift E Enf(xi) = f(xi+n) or Enyi = yi+n or En = E(E(E(....E))))
Forward  nf(xi) = n-1f(xi+1) – n-1f(xi) or nyi = n-1yi+1 – n-1yi or
difference n = (((....))))
Backward  nf(xi) = n-1f(xi) – n-1f(xi-1) or nyi = n-1yi – n-1yi-1 or
Difference n = (((....))))
Central  nf(xi) = n-1f(xi+½) – n-1f(xi-½) or nyi = n-1yi+½ – n-1yi-½ or
Difference  n = (((....))))
2f(xi) = f(xi+½) – f(xi-½) = [f(xi+½+½)- f(xi+½-½)] – [f(xi-½+½)- f(xi-½-½)] =
[f(xi+1)- f(xi)] – [f(xi)- f(xi-1)] = f(xi+1) - f(xi) – f(xi) + f(xi-1) = f(xi+1) – 2
f(xi) + f(xi-1)
2f(xi) = f(xi+1) -2 f(xi) + f(xi-1) or 2yi = yi+1 -2 yi + yi-1
4f(xi) = 2f(xi+1) -2 2f(xi) + 2f(xi-1) or 2yi = 2yi+1 -2 2yi + 2yi-1
nf(xi)=n-2f(xi+1)-2n-2f(xi)+n-2f(xi-1) or nyi=n-2yi+1-2n-2yi+n-2yi-1
Average Central  nf(xi) = (½)[ n-1f(xi+½) – n-1f(xi-½) ] or
Difference nyi = (½)(n-1yi+½ – n-1yi-½ ) or
 n = (((....))))
f(xi) = (½)[ f(xi+1) - f(xi-1) ] or yi = (½)(yi+1 - yi-1)
3f(xi) = (½)[ 2f(xi+1) - 2f(xi-1) ] or 3yi = (½)(2yi+1 - 2yi-1)
nf(xi)= (½)[ n-1f(xi+1) - n-1f(xi-1) ] or nyi= (½)(n-1yi+1 - n-1yi-1)

The successive applications of the numerical operators are calculated using the finite
difference tables shown in the tables below. Table 3 shows the table of forward differences,
Table 4 shows the table of backward differences and Table 5 shows the table of central
differences, where the powers of the even values  and the powers of the odd values . This
step is performed so that we only have integer indices of the image of the discrete functions
and the powers can be directly calculated from the given point table.

Table 3: Table of Forward Differences

i xi yi yi 2yi 3yi 4yi

1 x1 y1 y1 = y2 - y1 2y1 = y2 - y1 3y1 =2y2 - 2y1 4y1 =3y2 - 3y1
2 x2 y2 y2 = y3 – y2 2y2 = y3 –y2 3y2 = 2y3 –2y2 ----------
3 x3 y3 y3 = y4 - y3 2y3 = y4 - y3 ----------- ----------
4 x4 y4 y4 = y5 - y4 ----------- ----------- -----------
5 x5 y5 ---------- ----------- ------------ -----------
Table 4: Table of Backward Differences

i xi yi yi 2yi 3yi 4yi

1 x1 y1 ---------- ----------- ----------- ----------
2 x2 y2 y1 = y2 - y1 ----------- ------------ ----------
3 x3 y3 y2 = y3 – y2  y1 = y2 - y1
2 ----------- ----------
4 x4 y4 y3 = y4 - y3 2y2 = y3 –y2  y1 =2y2 - 2y1
3 -----------
5 x5 y5 y4 = y5 - y4 2y3 = y4 - y3 3y2 = 2y3 –2y2 4y1 =3y2 - 3y1

Table 5: Table of Central Differences

i xi yi yi 2yi 3yi 4yi

1 x1 y1 ---------- ----------- ----------- ----------
2 x2 y2 y2=(½)(y3 - y1)  y2=y3-2y2 +y1
2 ----------- ----------
3 x3 y3 y3=(½)(y4 – y2) 2y3=y4-2y3 +y2  y3=(½)(2y4 – 2y2)
3
 y3= y4 –22y3+2y2
4 2

4 x4 y4 y4=(½)(y5 - y3) 2y4=y5 -2y4 +y3 ----------- ----------

5 x5 y5 ---------- ----------- ----------- ----------

3) Review on Finite Difference Interpolation

The obtaining of the Forward Gregory-Newton interpolation method is showed[1,8,11].

Consider the x value of the abscissa where you want to interpolate, f(x) value to be interpolated,
xi any tabulated value before x, h = x (constant) spacing of the abscissa and z a new variable,
so that one has x = xi + z h where z indicates the fractional shift to, from xi, to achieve x. To
determine the value to be interpolated f(x) one takes y = f(x) = f(xi + z h) = Ez f(xi) = Ezyi = (1 + )z
∞
yi. One can then make an expansion using Newton's binomial (1 + ∆) 𝑧 = ∑ (𝑘𝑧 )∆𝑘 1𝑛−𝑘 ,
𝑘=0
such as this:

𝑧(𝑧−1) 2 𝑧(𝑧−1)(𝑧−2) 3 𝑧(𝑧−1)(𝑧−2)(𝑧−3) 4

𝑓(𝑥) = 𝑦𝑖 + 𝑧𝛥𝑦𝑖 + 2!
𝛥 𝑦𝑖 + 3!
𝛥 𝑦𝑖 + 4!
𝛥 𝑦𝑖 +. .. (2)

Alternatively,
∆𝑖 𝑦𝑖 𝑖−1 𝑥−𝑥𝑖
𝑓(𝑥) = 𝑦𝑖 + ∑∞
𝑖=1 𝑖!
∏𝑗=0(𝑧 − 𝑗) , 𝑧= ℎ
(3)

One then has in (2) and (3) the interpolation by Forward Gregory-Newton (at the beginning of
the table). The number of terms of the Forward Gregory–Newton series with the degree of the
power of  indicates the degree of the interpolating polynomial. Thus, if you use the series up
to the power one of , you have an interpolating polynomial of degree one; if you use the series
up to the power two of , you have an interpolator polynomial of degree two, etc.
Once the forward Gregory–Newton interpolator polynomial is obtained, one makes the
𝑑𝑓(𝑥) 𝑑𝑓(𝑥) 𝑑𝑧 𝑑𝑧 1
derivative of it by x using the chain rule, 𝑑𝑥
= 𝑑𝑧 𝑑𝑥
, where 𝑑𝑥
= ℎ, and to obtain the
derivative at a tabulated point, it is necessary to make z=0 at x = xi + z h, leaving x = xi. Thus, one
obtains (4) and (5):

1 ∆2 𝑦𝑖 ∆3 𝑦𝑖 ∆4 𝑦𝑖 ∆5 𝑦𝑖 ∆6 𝑦𝑖 ∆7 𝑦𝑖
𝑓´(𝑥𝑖 ) = (∆𝑦𝑖 − + − + − + − ⋯) (4)
ℎ 2 3 4 5 6 7

Alternatively,
1 ∆𝑘 𝑦𝑖
𝑓´(𝑥𝑖 ) = ℎ ∑∞
𝑘=1(−1)
(𝑘+1)
𝑘
(5)

An alternative obtainment of the Backward Gregory-Newton interpolation method is also

proposed. Consider the x value of the abscissa where you want to interpolate, f(x) value to be
interpolated, xi any value tabulated after x, h = x (constant) spacing of the abscissa and z a new
variable, so that one has x = xi – z h, where z indicates the fractional shift to, from xi, to achieve
x. To determine the value to be interpolated f(x) one takes y = f(x) = f(xi - z h) = E-z f(xi) = E-zyi = (1
- )z yi. One can then make an expansion using Newton's binomial (1 − ) 𝑧 =
∞
∑ (𝑘𝑧 )(−)𝑘 1𝑛−𝑘 , such as:
𝑘=0

𝑧(𝑧−1) 2 𝑧(𝑧−1)(𝑧−2) 3 𝑧(𝑧−1)(𝑧−2)(𝑧−3) 4

𝑓(𝑥) = 𝑦𝑖 − 𝑧∇𝑦𝑖 + 2!
∇ 𝑦𝑖 − 3!
∇ 𝑦𝑖 + 4!
∇ 𝑦𝑖 −⋯ (6)

Alternatively,
∇𝑖 𝑦𝑖 𝑖−1 𝑥𝑖 −𝑥
𝑓(𝑥) = 𝑦𝑖 + ∑∞ 𝑖
𝑖=1(−1) . 𝑖!
∏𝑗=0(𝑧 − 𝑗) , 𝑧= ℎ
(7)

We then have in (6) and (7) the interpolation by Backward Gregory-Newton (at the end of the
table). The number of terms of the backward Gregory–Newton series with the degree of the
power of  indicates the degree of the interpolating polynomial. Thus, if you use the series up
to the power one of , you have an interpolator polynomial of degree one; if you use the series
up to the power two of , you have an interpolator polynomial of degree two, etc.

Once the backward Gregory–Newton interpolator polynomial is obtained, one makes the
𝑑𝑓(𝑥) 𝑑𝑓(𝑥) 𝑑𝑧 𝑑𝑧 1
derivative of it by x using the chain rule, = where = − , and to obtain the
𝑑𝑥 𝑑𝑧 𝑑𝑥 𝑑𝑥 ℎ
derivative at a tabulated point, it is necessary to make z=0 at x = xi – z h, leaving x = xi,. Thus, one
obtains (8) and (9):

1 ∇2 𝑦𝑖 ∇3 𝑦𝑖 ∇4 𝑦𝑖 ∇5 𝑦𝑖 ∇6 𝑦𝑖 ∇7 𝑦𝑖
𝑓´(𝑥𝑖 ) = ℎ (∇𝑦𝑖 + 2
+ 3
+ 4
+ 5
+ 6
+ 7
+ ⋯) (8)

Alternatively,
1 ∇𝑘 𝑦𝑖
𝑓´(𝑥𝑖 ) = ℎ ∑∞
𝑘=1 𝑘
(9)
An alternative obtainment of the Stirling interpolation method is again proposed. Consider the
x value of the abscissa where you want to interpolate, f(x) value to be interpolated, xi a value
any tabulated before x, h = x (constant) spacing of the abscissa and z a new variable, so that
one has x = xi + z h where z indicates the fractional shift to, from xi, to achieve x. You can use
the relations between operators: [1 +  + 2/2]= [1 + (½)(E1 - E-1) + (½)(E1/2 – E-1/2)2] = [ 1 + ½E -
½E-1 + (½)(E – 2 + E-1) ] = [ 1 + ½E - ½E-1 + ½E – 1 + ½E-1 ] = E. To determine the value to be
interpolated f(x) one takes: y = f(x) = f(xi + z h) = Ez f(xi) = Ezyi = (1 +  + 2/2)z yi ; one can then
𝑧 ∞
make an expansion using Newton's binomial (1 +  + 𝟐 /2) = ∑ (𝑘𝑧 )( + 𝟐 /
𝑘=0
2)𝑘 1𝑛−𝑘 , like this:

𝑧2 2 𝑧(𝑧 2 −12 ) 𝑧 2 (𝑧 2 −12 ) 𝑧(𝑧 2 −1)(𝑧 2 −22 )

𝑓(𝑥) = 𝑦𝑖 + 𝑧𝜇𝛿𝑦𝑖 + 𝛿 𝑦𝑖 + 𝜇𝛿 3 𝑦𝑖 + 𝛿 4 𝑦𝑖 + 𝜇𝛿 5 𝑦𝑖 +
2! 3! 4! 5!
𝑧 2 (𝑧 2 −12 )(𝑧 2 −22 ) 𝑧(𝑧 2 −12 )(𝑧 2 −22 )(𝑧 2 −32 ) 𝑧 2 (𝑧 2 −12 )(𝑧 2 −4)(𝑧 2 −32 )
𝛿 6 𝑦𝑖 + 𝜇𝛿 7 𝑦𝑖 + 𝛿 8 𝑦𝑖 + ⋯ (10)
6! 7! 8!

Alternatively,
𝑧+𝑘−1 𝑧 𝑥−𝑥𝑖
𝑓(𝑥) = 𝑦𝑖 + ∑∞
𝑘=1 ( ) (𝜇𝛿 2𝑘−1 𝑦𝑖 + 𝛿 2𝑘 𝑦𝑖 ) , 𝑧= (11)
2𝑘 − 1 2𝑘 ℎ

Alternatively,
𝜇𝛿 (2𝑖−1) 𝑦 𝛿 (2𝑖) 𝑦𝑖 𝑥−𝑥𝑖
𝑓(𝑥) = 𝑦𝑖 + ∑∞ 𝑖
𝑖=1 ( 𝑧 (2𝑖−1)! + (2𝑖)!
) ∏𝑖−1
𝑗=0(𝑧
2
− 𝑗2) , 𝑧= ℎ
(12)

One then has in (10), (11) and (12) the interpolation by Stirling (in the middle of the table). The
number of terms of the Stirling series with the degree of the power of  indicates the degree of
the interpolating polynomial. Thus, if you use the series up to the power one of , you have an
interpolator polynomial of degree one; if you use the series up to the power two of , you have
an interpolator polynomial of degree two, etc.

Once the Stirling interpolator polynomial is obtained, one makes the derivative of it by x using
𝑑𝑓(𝑥) 𝑑𝑓(𝑥) 𝑑𝑧 𝑑𝑧 1
the chain rule, 𝑑𝑥
= 𝑑𝑧 𝑑𝑥
where 𝑑𝑥 = ℎ, and to obtain the derivative at a tabulated point,
it is necessary to make z=0 at x = xi + z h, leaving x = xi, and thus one obtains (13) and (14):

1 𝜇𝛿 3 𝑦𝑖 𝜇𝛿 5 𝑦𝑖 𝜇𝛿 7 𝑦𝑖 𝜇𝛿 9 𝑦𝑖 𝜇𝛿 11 𝑦𝑖 𝜇𝛿 13 𝑦𝑖
𝑓´(𝑥𝑖 ) = ℎ (𝜇𝛿𝑦𝑖 − 6
+ 30
− 140
+ 630
− 2772
+ 12012
−⋯) (13)

Alternatively,
1 𝜇𝛿 (2𝑘−1) 𝑦𝑖 𝑘−1
𝑓´(𝑥𝑖 ) = ℎ ∑∞
𝑘=1(−1)
(𝑘+1)
(2𝑘−1)!
∏𝑗=0 (𝑗 + 1)2 (14)
 4) Obtaining Equations by Numerical Derivation Using Finite Differences

Once the forward, backward, and central approximate formulas for first derivatives are
obtained, these can be algebraically manipulated using the relationships between the operators
and the shift operator (E) and with notable products. In the above formulas for the first
derivative f (́ xi), the powers of ,  and  define the order of the error in the calculation of the
derivatives. Thus, if you expand the series to power one, you have a derivative formula with an
error order equal to one; if you expand the series to power two, you have a drift formula with
an error order equal to two. For these Gregory-Newton and Stirling first derivative formulas,
there is a simpler way to determine the order of error of the formulas of the first derivative,
which is to subtract the highest index of yi from the lowest index of yi. The number resulting
from the subtraction indicates the order of the error. Therefore, for example,

1 1 1 1 𝑦𝑖+1 −𝑦𝑖
𝑓´(𝑥𝑖 ) = ℎ (∆𝑦𝑖 ) = ℎ (∆)𝑦𝑖 = ℎ (𝐸 − 1)𝑦𝑖 = ℎ (𝐸𝑦𝑖 − 𝑦𝑖 ) = ℎ
(15)

Error=O(h1-0) = O(h)

Alternatively,
2
1 ∇2 𝑦𝑖 1 ∇2 1 (1−𝐸 −1 )
𝑓´(𝑥𝑖 ) = (∇𝑦𝑖 + ) = (∇ + ) 𝑦𝑖 = ((1 − 𝐸 −1 ) + ) 𝑦𝑖 =
ℎ 2 ℎ 2 ℎ 2

1 1 1
((1 − 𝐸 −1 ) + (1 − 2𝐸 −1 + 𝐸 −2 )) 𝑦𝑖 = ((2 − 2𝐸 −1 ) + (1 − 2𝐸 −1 + 𝐸 −2 ))𝑦𝑖 (16)
ℎ 2 2ℎ
1 1 3𝑦𝑖 − 4𝑦𝑖−1 + 𝑦𝑖−2
= (3 − 4𝐸 −1 + 𝐸 −2 )𝑦𝑖 = (3𝑦𝑖 − 4𝐸 −1 𝑦𝑖 + 𝐸 −2 𝑦𝑖 ) =
2ℎ 2ℎ 2ℎ
Error = O(h 0-(-2)) = O(h2)

Alternatively,
1 𝜇𝛿 3 𝑦𝑖 1 𝜇𝛿 3 1 𝜇𝛿 𝛿 2
𝑓´(𝑥𝑖 ) = (𝜇𝛿𝑦𝑖 − ) = (𝜇𝛿 − ) 𝑦𝑖 = (𝜇𝛿 − ) 𝑦𝑖
ℎ 6 ℎ 6 ℎ 6
1 1 1 1
= ℎ (2 (𝐸 − 𝐸 −1 ) − 6 (2 (𝐸 − 𝐸 −1 )(𝐸 − 2 + 𝐸 −1 ))) 𝑦𝑖
1 1 1
= ℎ (2 (𝐸 − 𝐸 −1 ) − 12 (𝐸 2 − 2𝐸 + 1 − 1 + 2𝐸 −1 − 𝐸 −2 )) 𝑦𝑖 (17)
1
= 12ℎ (6𝐸 − 6𝐸 −1 − (𝐸 2 − 2𝐸 + 2𝐸 −1 − 𝐸 −2 ))𝑦𝑖
1 1
= (−𝐸 2 + 8𝐸 − 8𝐸 −1 + 𝐸 −2 )𝑦𝑖 = (−𝐸 2 𝑦𝑖 + 8𝐸𝑦𝑖 − 8𝐸 −1 𝑦𝑖 + 𝐸 −2 𝑦𝑖 )
12ℎ 12ℎ
−𝑦𝑖+2 + 8𝑦𝑖+1 − 8𝑦𝑖−1 + 𝑦𝑖−2
=
12ℎ
Error = O(h2-(-2)) = O(h4)

In (15), (16) and (17) show forward, backward, and central transformations, respectively, of the
formulas of derivatives using numerical operators in finite difference equations. With the use of
numerical operators, to obtain finite difference formulas for approximation of derivatives, it
becomes a work of the relationship between two operators and purely algebraic manipulation,
which can be automatically achieved by mathematical algebra software for any error order.
Tables 6 through 13 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 x1 or x2 or x3 or... for dependent variable, as well as the index i by i1 or i2 or i3 or... for the
indices of each dependent variable.

Table 6: Finite difference formulas with error orders from O(h) through O(h4).
Table 7: Finite difference formulas with error orders from O(h5) to O(h8).
Table 8: Finite difference formulas with error orders from O(h9) through O(h12).
Table 9: Finite difference formulas with error orders from O(h13) to O(h15).
Table 10: Finite difference formulas with error orders from O(h16) to O(h17).
Table 11: Finite difference formulas with error orders from O(h18) to O(h19).
Table 12: Finite difference formulas with error orders of O(h20).

Table 13: Finite difference formulas with error orders of O(h30).

It is observed that in Table 13, 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)[12].

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 (18), one has:
𝑑𝑢(𝑥1) −𝑢𝑖1−1 +𝑢𝑖1+1 2) 𝑑𝑢(𝑥1) −𝐸 −1 +𝐸
𝑑𝑥1
= 2∆𝑥1
, 𝑂(∆𝑥1 → 𝑑𝑥1
= ( 2∆𝑥1
) 𝑢𝑖1
2
𝑑 2 𝑢(𝑥1) −𝐸 −1 +𝐸 𝐸 −2 −2+𝐸 2 𝐸 −2 𝑢𝑖1 −2𝑢𝑖1 +𝐸 2 𝑢𝑖1
𝑑𝑥12
= ( 2∆𝑥1 ) 𝑢𝑖1 = ( 4∆𝑥12 ) 𝑢𝑖1 = ( 4∆𝑥12
) =
𝑢𝑖1−2 −2𝑢𝑖1 +𝑢𝑖1+2
( 4∆𝑥12
) , 𝑂(∆𝑥12 ) (18)

Alternatively, for a second order centered third derivative (19), one has:
𝑑𝑢(𝑥1) −𝑢𝑖1−1 +𝑢𝑖1+1 2) 𝑑𝑢(𝑥1) −𝐸 −1 +𝐸
𝑑𝑥1
= 2∆𝑥1
, 𝑂(∆𝑥1 → 𝑑𝑥1
= ( 2∆𝑥1
) 𝑢𝑖1
3
𝑑 3 𝑢(𝑥1) −𝐸 −1 +𝐸 −𝐸 −3 +3𝐸 −1 −3𝐸+𝐸 3 −𝐸 −3 𝑢𝑖1 +3𝐸 −1 𝑢𝑖1 −3𝐸𝑢𝑖1 +𝐸 3𝑢𝑖1
𝑑𝑥1 3 = ( 2∆𝑥1
) 𝑢 𝑖1 = ( 8∆𝑥1 3 ) 𝑢 𝑖1 = ( 8∆𝑥13
) =
−𝑢 +3𝑢 −3𝑢 +𝑢
( 𝑖1−3 𝑖1−1 3 𝑖1+1 𝑖1+3 ) , 𝑂(∆𝑥12 ) (19)
8∆𝑥1

In another example, for a second order centered partial derivative with second-order error (20),
one has:
𝑑𝑢(𝑥1) −𝑢𝑖1−1 +𝑢𝑖1+1 2) 𝑑𝑢(𝑥1) −𝐸𝑥1 −1 +𝐸𝑥1
𝑑𝑥1
= 2∆𝑥1
, 𝑂(∆𝑥1 → 𝑑𝑥1
= ( 2∆𝑥1
) 𝑢𝑖1
−1 −1
2
𝜕 𝑢(𝑥1,𝑥2) −𝐸𝑥1 +𝐸𝑥1 −𝐸𝑥2 +𝐸𝑥2 𝐸𝑥1 𝐸𝑥2 −𝐸𝑥1 −1 𝐸𝑥2 −𝐸𝑥1 𝐸𝑥2 −1 +𝐸𝑥1 𝐸𝑥2
−1 −1

𝜕𝑥1𝜕𝑥2
= ( 2∆𝑥1
) ( 2∆𝑥2
) 𝑢 𝑖1,𝑖2 = ( 4∆𝑥1∆𝑥2
) 𝑢𝑖1,𝑖2 =
−1 −1 −1 −1
𝐸𝑥1 𝐸𝑥2 𝑢𝑖1,𝑖2 −𝐸𝑥1 𝐸𝑥2 𝑢𝑖1,𝑖2 −𝐸𝑥1 𝐸𝑥2 𝑢𝑖1,𝑖2 +𝐸𝑥1 𝐸𝑥2 𝑢𝑖1,𝑖2
4∆𝑥1∆𝑥2
=
𝑢𝑖1−1,𝑖2−1 −𝑢𝑖1−1,𝑖2+1 −𝑢𝑖1+1,𝑖2−1 +𝑢𝑖1+1,𝑖2+1
4∆𝑥1∆𝑥2
, 𝑂(∆𝑥12 ∆𝑥22 ) (20)

For another example, of a third order centered derivative with two independent variables (21),
one has:
𝑑𝑢(𝑥1) −𝑢𝑖1−1 +𝑢𝑖1+1 𝑑𝑢(𝑥1) −𝐸𝑥1 −1 +𝐸𝑥1
𝑑𝑥1
= 2∆𝑥1
, 𝑂(∆𝑥12 ) → 𝑑𝑥1
=( 2∆𝑥1
) 𝑢𝑖1
2
𝜕3 𝑢(𝑥1,𝑥2) −𝐸𝑥1 −1 +𝐸𝑥1 −𝐸𝑥2 −1 +𝐸𝑥2
2
𝜕𝑥1 𝜕𝑥2
= ( 2∆𝑥1
) ( 2∆𝑥2
) 𝑢𝑖1,𝑖2 =
−𝐸 −2 𝐸 −1 +𝐸𝑥1 𝐸𝑥2 +2𝐸𝑥2 −2𝐸𝑥2 −𝐸𝑥1 𝐸𝑥2 −1 +𝐸𝑥1 2 𝐸𝑥2
−2 −1 2
( 𝑥1 𝑥2 8∆𝑥12 ∆𝑥2
) 𝑢𝑖1,𝑖2 =
−𝐸𝑥1 𝐸𝑥2 𝑢𝑖1,𝑖2 +𝐸𝑥1 𝐸𝑥2 𝑢𝑖1,𝑖2 +2𝐸𝑥2 𝑢𝑖1,𝑖2 −2𝐸𝑥2 𝑢𝑖1,𝑖2 −𝐸𝑥1 2 𝐸𝑥2 −1 𝑢𝑖1,𝑖2 +𝐸𝑥1 2 𝐸𝑥2 𝑢𝑖1,𝑖2
−2 −1 −2 −1
=
8∆𝑥12 ∆𝑥2
−𝑢𝑖1−2,𝑖2−1 +𝑢𝑖1−2,𝑖2+1 +2𝑢𝑖1,𝑖2−1 −2𝑢𝑖1,𝑖2+1 −𝑢𝑖1+2,𝑖2−1 +𝑢𝑖1+2,𝑖2+1
8∆𝑥12 ∆𝑥2
, 𝑂(∆𝑥12 ∆𝑥22 ) (21)

Now, for a third order centered partial derivative with second-order error and three
independent variables (22), we have:
𝑑𝑢(𝑥1) −𝑢𝑖1−1 +𝑢𝑖1+1 2) 𝑑𝑢(𝑥1) −𝐸𝑥1 −1 +𝐸𝑥1
𝑑𝑥1
= 2∆𝑥1
, 𝑂(∆𝑥1 → 𝑑𝑥1
= ( 2∆𝑥1
) 𝑢𝑖1
3 −1 −1 −1
𝜕 𝑢(𝑥1,𝑥2,𝑥3) −𝐸 +𝐸 −𝐸 +𝐸 −𝐸 +𝐸
𝜕𝑥1𝜕𝑥2𝜕𝑥3
= ( 𝑥12∆𝑥1 𝑥1 ) ( 𝑥22∆𝑥2 𝑥2 ) ( 𝑥32∆𝑥3 𝑥3 ) 𝑢𝑖1,𝑖2,𝑖3 =
𝐸 −1 𝐸 −1 𝐸 −1 +𝐸𝑥1 −1 𝐸𝑥2 −1 𝐸𝑥3 +𝐸𝑥1 −1 𝐸𝑥2 𝐸𝑥3 −1 −𝐸𝑥1 −1 𝐸𝑥2 𝐸𝑥3 +𝐸𝑥1 𝐸𝑥2 −1 𝐸𝑥3 −1 −𝐸𝑥1 𝐸𝑥2 −1 𝐸𝑥3 −𝐸𝑥1 𝐸𝑥2 𝐸𝑥3 −1 +𝐸𝑥1 𝐸𝑥2 𝐸𝑥3
( 𝑥1 𝑥2 𝑥3 8∆𝑥1∆𝑥2∆𝑥3
)
𝑢𝑖1,𝑖2,𝑖3 =
−𝑢𝑖1−1,𝑖2−1,𝑖3−1 +𝑢𝑖1−1,𝑖2−1,𝑖3+1 +𝑢𝑖1−1,𝑖2+1,𝑖3−1 −𝑢𝑖1−1,𝑖2+1,𝑖3+1 +𝑢𝑖1+1,𝑖2−1,𝑖3−1 −𝑢𝑖1+1,𝑖2−1,𝑖3+1 −𝑢𝑖1+1,𝑖2+1,𝑖3−1 +𝑢𝑖1+1,𝑖2+1,𝑖3+1
8∆𝑥1∆𝑥2∆𝑥3
, 𝑂(∆𝑥12 ∆𝑥22 ∆𝑥32 ) (22)
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 uip+n = Enuip 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 Ex1n1Ex2n2Ex3n3... Expnpui =
ui1+n1,i2+n2,i3+n3,...,ip+np. Using computer algebra software makes it even easier to obtain finite
difference approximations of successive derivatives of any order and of any number of
independent variables. As supplementary material to the paper, many tried and ready-to-use
approximations are shown.

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 use a formula with order error two for x1 and with order error
three for x2. 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) Comparison with Techniques for Calculating Finite Difference Coefficients using Taylor
Expansion

Finite difference equations enable you to take derivatives of any order at any point using any
given sufficiently large selection of points. By inputting the locations of your sampled points
below, you will generate a finite difference equation which will approximate the derivative at
any desired location. To take a numerical derivative, you calculate the slope at a point using the
values and relative locations of surrounding points. The locations of these sampled points are
collectively called the finite difference stencil[13–23].

This calculator accepts as input any finite difference stencil and desired derivative order and
dynamically calculates the coefficients for the finite difference equation. You may be familiar
𝜕𝑓 𝑓(𝑥)−𝑓(𝑥−ℎ)
with the backward difference derivative 𝜕𝑥
= ℎ
. This is a special case of a finite
difference equation (where f(x) – f(x-h) is the finite difference and h is the spacing between the
points) and can be displayed below by entering the finite difference stencil {-1,0} for Locations
of Sampled Points and 1 for Derivative Order. The finite difference coefficients calculator can be
used generally for any finite difference stencil and any derivative order. Notable cases include
the forward difference derivative, {0,1} and 1, the second-order central difference, {-1,0,1} and
2, and the fourth-order five-point stencil, {-2,-1,0,1,2} and 4.

For example, to determine the coefficients for difference equation:

𝜕4 𝑓
𝜕𝑥 4
≈ 𝐴𝑓(𝑥 − 2ℎ) + 𝐵𝑓(𝑥 − ℎ) + 𝐶𝑓(𝑥) + 𝐷𝑓(𝑥 + ℎ) + 𝐸𝑓(𝑥 + 2ℎ) (23)
The formula for Taylor expansion is
𝑓´(𝑎) 𝑓´´(𝑎) 𝑓(𝑛) (𝑎)
𝑓(𝑏) = 𝑓(𝑎) + 1!
(𝑏 − 𝑎) + 2!
(𝑏 − 𝑎)2 + ⋯ + 𝑛!
(𝑏 − 𝑎)𝑛 + ⋯ (24)

We began by expanding the first term, A f(x-2h), with b=x-2h and a=x, getting
1 1 1
𝐴𝑓(𝑥) + 𝐴𝑓´(𝑥)(−2ℎ) + 𝐴 2 𝑓´´(𝑥)(−2ℎ)2 + 𝐴 6 𝑓´´´(𝑥)(−2ℎ)3 + 𝐴 24 𝑓´´´´(𝑥)(−2ℎ)4 +
1
𝐴 120 𝑓´´´´´(𝑥)(−2ℎ)5 + ⋯ (25)

We then did the same for all the other terms in the numerator:
1 1 1
+𝐵𝑓(𝑥) + 𝐵𝑓´(𝑥)(−ℎ) + 𝐵 2 𝑓´´(𝑥)(−ℎ)2 + 𝐵 6 𝑓´´´(𝑥)(−ℎ)3 + 𝐵 24 𝑓´´´´(𝑥)(−ℎ)4 +
1
𝐵 𝑓´´´´´(𝑥)(−ℎ)5 +⋯
120
+𝐶𝑓(𝑥)
1 1 1
+𝐷𝑓(𝑥) + 𝐷𝑓´(𝑥)(ℎ) + 𝐷 2 𝑓´´(𝑥)(ℎ)2 + 𝐷 6 𝑓´´´(𝑥)(ℎ)3 + 𝐷 24 𝑓´´´´(𝑥)(ℎ)4 +
1
𝐷 120 𝑓´´´´´(𝑥)(ℎ)5 + ⋯ (26)
1 1 1
+𝐸𝑓(𝑥) + 𝐸𝑓´(𝑥)(2ℎ) + 𝐸 2 𝑓´´(𝑥)(2ℎ)2 + 𝐸 6 𝑓´´´(𝑥)(2ℎ)3 + 𝐸 24 𝑓´´´´(𝑥)(2ℎ)4 +
1
𝐸 120 𝑓´´´´´(𝑥)(2ℎ)5 +⋯

Now since we needed to choose A, B, C, D and E such that the summation of all of these terms
results in the cancellation of all f(x), f´(x), f´´(x), and f´´´(x) terms and results in the f´´´´(x)
coefficients summing to one, this requires us to fulfill the equations:
𝐴+𝐵+𝐶+𝐷+𝐸 =0
−2𝐴 − 𝐵 + 𝐷 + 2𝐸 = 0
4𝐴 + 𝐵 + 𝐷 + 4𝐸 = 0 (27)
−8𝐴 − 𝐵 + 𝐷 + 8𝐸 = 0
24
{16𝐴 + 𝐵 + 𝐷 + 16𝐸 = ⁄ℎ4
These can be rewritten in general form as
4!
(−2)𝑛 𝐴 + (−1)𝑛 𝐵 + 𝐶𝜑(𝑛) + (1)𝑛 𝐷 + (2)𝑛 𝐸 = 𝜑(𝑛 − 4) 𝑓𝑜𝑟 0 ≤ 𝑛 ≤ 4 (28)
ℎ4

Or in matrix form as

1 1 1 1 1   A 0
 −2 −1 0  
1 2  B  0
 1  
4 1 0 1 4  C  = 4  0  (29)
   h  
 −8 −1 0 1 8   D 0
16 1 0 1 16   E   24 

Now this matrix is invertible, and the solution is:

−1
 A 1 1 1 1 1 0 1
B  2  0  
  1  −2 −1 0 1   1  −4 
C  = 4  4 1 0 1 4 0= 46 (30)
  h     h  
D  −8 −1 0 1 8 0  −4 
 E  16 1 0 1 16   24   1 
Substituting these values back into the original equation, we have the difference equation:
𝜕4 𝑓 𝑓(𝑥−2ℎ)−4𝑓(𝑥−ℎ)+6𝑓(𝑥)−4𝑓(𝑥+ℎ)+𝑓(𝑥+2ℎ)
𝜕𝑥 4
≈ ℎ4
(31)

Carefully considering the general form above, it is noticed that this equation can be generalized
to obtain the finite difference equation from any finite difference stencil given the desired
derivative order (provided the desired finite difference equation exists for the given derivative
order). Given a stencil s of length N and derivative order d < N, the coefficients c’s are given by
the finite difference coefficients equation:
𝑑!
𝑠1𝑛 𝑐1 + 𝑠2𝑛 𝑐2 + ⋯ + 𝑠𝑁𝑛 𝑐𝑁 = 𝜑(𝑛 − 𝑑) 𝑓𝑜𝑟 0 ≤ 𝑛 ≤ 𝑁 − 1 (32)
ℎ𝑑

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

0
𝑐1 𝑠10 ⋯ 𝑠𝑁0 ⋮
1
[ ⋮ ] = ℎ𝑑 [ ⋮ ⋱ ⋮ ] 𝑑! (33)
𝑐𝑁 𝑠1𝑁−1 ⋯ 𝑠𝑁𝑁−1 ⋮
[0]
Each time s and d are altered in the calculator above, this matrix equation is solved and the
resulting finite difference equation is displayed.

It can be noted that this technique for calculating finite difference coefficients consumes much
more computational processing time than the technique described in the previous section. In
other words, the technique using numerical operators obtains the same results with much fewer
mathematical operations.

6) Changes in Formula Indices for Shift of Points and Variables

One point that may go unnoticed is that because the step (x) is constant, one can change the
values of the indices in the finite difference formulas. Thus, for example, for Formula (34):

𝑑 2 𝑢(𝑥1) 𝑢𝑖1−4 −16𝑢𝑖1−3 +64𝑢𝑖1−2 +16𝑢𝑖1−1 −130𝑢𝑖1 +16𝑢𝑖1+1 +64𝑢𝑖1+2 −16𝑢𝑖1+3 +𝑢𝑖1+4
𝑑𝑥12
= 144∆𝑥12
, 𝑂(∆𝑥14 ) (34)

One can swap i1 for i1+4 without changing the approximation by finite differences or the order
of error, leaving (24):

𝑑 2 𝑢(𝑥1) 𝑢𝑖1 −16𝑢𝑖1+1 +64𝑢𝑖1+2 +16𝑢𝑖1+3 −130𝑢𝑖1+4 +16𝑢𝑖1+5 +64𝑢𝑖1+6 −16𝑢𝑖1+7 +𝑢𝑖1+8
𝑑𝑥12
= 144∆𝑥12
, 𝑂(∆𝑥14 ) (35)

In another example (36), one can exchange the derivatives for partial derivatives, as follows:

𝑑 3 𝑢(𝑥1) −𝑢𝑖1−3 +3𝑢𝑖1−1 −3𝑢𝑖1+1 +𝑢𝑖1+3

𝑑𝑥13
= 8∆𝑥13
, 𝑂(∆𝑥12 ) (36)

It can be exchanged for (37) without changing the formula coefficients or the order of the
error:
𝜕3 𝑢(𝑥1,𝑥2) −𝑢𝑖1−3,𝑖2 +3𝑢𝑖1−1,𝑖2 −3𝑢𝑖1+1,𝑖2 +𝑢𝑖1+3,𝑖2
𝜕𝑥13
= 8∆𝑥13
, 𝑂(∆𝑥12 ) (37)

One can also swap the variable with respect to the derivative, for example (38):

𝜕3 𝑢(𝑥1,𝑥2) −𝑢𝑖1,𝑖2−3 +3𝑢𝑖1,𝑖2−1 −3𝑢𝑖1,𝑖2+1 +𝑢𝑖1,𝑖2+3

𝜕𝑥23
= 8∆𝑥23
, 𝑂(∆𝑥22 ) (38)

In another case, one can change the order of the variables in the partial derivatives (39):

𝜕3 𝑢(𝑥1,𝑥2) −𝑢𝑖1−2,𝑖2−1 +𝑢𝑖1−2,𝑖2+1 +2𝑢𝑖1,𝑖2−1 −2𝑢𝑖1,𝑖2+1 −𝑢𝑖1+2,𝑖2−1 +𝑢𝑖1+2,𝑖2+1

𝜕𝑥12 𝜕𝑥2
= 8∆𝑥12 ∆𝑥2
, 𝑂(∆𝑥12 ∆𝑥22 ) (39)

Swapping the indices becomes between i1 and i2 in (40):

𝜕3 𝑢(𝑥1,𝑥2) −𝑢𝑖1−1,𝑖2−2 +𝑢𝑖1+1,𝑖2−2 +2𝑢𝑖1−1,𝑖2 −2𝑢𝑖1+1,𝑖2 −𝑢𝑖1−1,𝑖2+2 +𝑢𝑖1+1,𝑖2+2

𝜕𝑥1𝜕𝑥22
= 8∆𝑥1∆𝑥22
, 𝑂(∆𝑥12 ∆𝑥22 ) (40)

Ordering the terms, one has (41):

𝜕3 𝑢(𝑥1,𝑥2) −𝑢𝑖1−1,𝑖2−2 +2𝑢𝑖1−1,𝑖2 −𝑢𝑖1−1,𝑖2+2 +𝑢𝑖1+1,𝑖2−2 −2𝑢𝑖1+1,𝑖2 +𝑢𝑖1+1,𝑖2+2

𝜕𝑥1𝜕𝑥22
= 8∆𝑥1∆𝑥22
, 𝑂(∆𝑥12 ∆𝑥22 ) (41)

Thus, the forward equations can become backward and vice versa, using symmetry as shown in
(42), (43) and (44). In (42) one has the coefficients of the numerators (-1 and 1) in the forward
equation and (-1 and 1) in the backward equation. In (42), the numerator one always has the
variation of u (u), that is, the value of the u of the next point minus the value of the u of the
previous point:
𝑑 −𝑢𝑖1 +𝑢𝑖1+1
𝑢(𝑥1) = , 𝑂( ∆𝑥1 ), 𝐹𝑜𝑟𝑤𝑎𝑟𝑑
𝑑𝑥1 ∆𝑥1
𝑑 −𝑢𝑖1−1 +𝑢𝑖1
𝑑𝑥1
𝑢(𝑥1) = ∆𝑥1
, 𝑂( ∆𝑥1 ), 𝐵𝑎𝑐𝑘𝑤𝑎𝑟𝑑 (42)

In (43), one has the coefficients of the numerators (-3, 4 and -1) in the forward equation and the
opposite in the backward equation with the coefficients of the numerators being (1, -4 and 3):
𝑑 −3𝑢𝑖1 +4𝑢𝑖1+1 −𝑢𝑖1+2
𝑑𝑥1
𝑢(𝑥1) = 2∆𝑥1
, 𝑂( ∆𝑥12 ), 𝐹𝑜𝑟𝑤𝑎𝑟𝑑
𝑑 𝑢𝑖1−2 −4𝑢𝑖1−1 +3𝑢𝑖1
𝑢(𝑥1) = , 𝑂( ∆𝑥12 ), 𝐵𝑎𝑐𝑘𝑤𝑎𝑟𝑑 (43)
𝑑𝑥1 2∆𝑥1

In (44), the coefficients of the numerators are (-11, 18, -9 and 2) in the forward equation and
the opposite in the backward equation, where the coefficients of the numerators are (-2, 9, -18,
11). Thus, having the forward formula, one can obtain the backward formula of finite differences
by symmetry. The indices change from (i1, i1+1, i1+2 and i1+3) in the forward equation to (i1-3,
i1-2, i1-1 and i1) in the backward equation:

(44)
The sign exchange occurs in the numerators above because the step or addition in the
denominator also changes from h or x or x1 in the forward equation to -h or -x or -x1 in
the backward equation. Instead of using the front points, greater than x1i1, the backward
formula uses the back points, smaller than x1i1. These exchanges of indices and variables save
the deduction of new finite difference formulas and simplify the use of finite difference formulas
at the boundaries of the geometric domains of the mesh.
Symmetry is also observed in the coefficients of the numerators of the central finite difference
equations, as shown here:

(45)

where in (45), the first half numerator coefficients (3, -32, 168 and -672) are symmetrical to the
second half coefficients (672, -168, 32 and 3), just by changing the sign. In another example:

(46)

where in (46), the first half numerator coefficients (1, -16, 64 and 16) are symmetric to the
second half coefficients (16, 64, -16 and 1). In this last example:

(47)

where in (47), the first half numerator coefficients (1, -8, 8, -1, -8, 64, -64 and 8) are perfectly
symmetrical to the second half coefficients (8, -64, 64, -8, -1, 8, -8 and 1).
This result shows how robust and practical the approximation equations of derivatives by finite
differences are, 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, (34) and (35),
b) use the formulas obtained for total and partial derivatives, (36) and (37),
c) swap the variable with respect to the derivative, (38) and (39),
d) change the orders of the variables of the partial derivatives, changing the indices of
Equations (40) and (41),
e) obtain the backward equations based on the forward equations and vice versa, (42),
(43) and (44),
f) observe the symmetry of the numerator coefficients of centered equations, (45),
(46) and (47).

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, one arrives at linear systems that can be revolved
most of the time by iterative methods or by successive and cyclic application of the equations
as utilized in spreadsheets.
 7) New Way to Write Finite Difference Formulas

Often, one has to use many different finite difference formulas to solve a given practical problem
involving differential equations. In addition, one has many approximations of finite differences
that can be applied, and these formulas have many coefficients and variations. The idea is to be
able to assemble a database with many approximations of derivatives by finite differences,
which will be used in the solution of differential equations with different boundary conditions.
Thus, with an alternative to represent, save, and load the approximations by finite differences,
we propose a table that can also be stored in a database and even in text files with the columns
separated by semicolons (;), thus given to a formula of general finite difference (48):

Number equation: neq (Integer sequential number)

𝜕𝑚 𝑢(𝑥1,𝑥2,𝑥3,…,𝑥𝑣) 1
𝜕𝑥1𝑚1 𝜕𝑥2𝑚2 𝜕𝑥3𝑚3 …𝜕𝑥𝑣 𝑚𝑣
= (∆𝑥1𝑚1 ∆𝑥2𝑚2 ∆𝑥3𝑚3 …∆𝑥𝑣 𝑚𝑣 ).
𝑎1 𝑎
( 𝑢
𝑏1 𝑖1+𝑛11,𝑖2+𝑛12,𝑖3+𝑛13,…,𝑖𝑣+𝑛1𝑣
+ 𝑏2 𝑢𝑖1+𝑛21,𝑖2+𝑛22,𝑖3+𝑛23,…,𝑖𝑣+𝑛2𝑣 +
2
𝑎3 𝑎
𝑢
𝑏3 𝑖1+𝑛31,𝑖2+𝑛32,𝑖3+𝑛33,…,𝑖𝑣+𝑛3𝑣
+ ⋯ + 𝑏𝑡 𝑢𝑖1+𝑛𝑡1,𝑖2+𝑛𝑡2,𝑖3+𝑛𝑡3,…,𝑖𝑣+𝑛𝑡𝑣 ) (48)
𝑡

𝑂(∆𝑥1𝑜1 ∆𝑥2𝑜2 ∆𝑥3𝑜3 … ∆𝑥𝑣 𝑜𝑣 )

where  = number of independent variables (x1, x2, x3,..., x);

t = number of terms that are summed in the finite difference formula;

m = total order of the derivative to be approximated;

m1, m2, m3, ..., m = partial order of the differentials of each independent variable (x1, x2,
x3, ..., x);

a1, a2, a3, ..., at = numerators of each term of the finite difference formula;

b1, b2, b3, ..., bt = denominators of each term of the finite difference formula;

n11, n12, n13, ..., n1 = indices of the independent variables of the first term of the finite
difference formula;

n21, n22, n23, ..., n2 = indices of the independent variables of the second term of the finite
difference formula;

n31, n32, n33, ..., n3 = indices of the independent variables of the third term of the finite
difference formula, etc., until;

nt1, nt2, nt3, ..., nt = indices of the independent variables of the last term of the finite
difference formula;

o1, o2, o3, ..., o = error order of the finite difference equation of each independent
variable (x1, x2, x3, ..., x);

i1 is the index variable for the first independent variable x1, i2 is the index variable for the
second independent variable x2, i3 is the index variable for the third independent variable x3,
etc. until i is the index variable for the last independent variable x, which is the dependent
variable to be derived.
This finite difference approximation is stored in the table as follows:

neq  o1 o2 o3 … o
t m m1 m2 m3 ... m
a1 b1 n11 n12 n13 ... n1
a2 b2 n21 n22 n23 ... n2
a3 b3 n31 n32 n33 ... n3
... ... ... ... ... ... ...
at bt nt1 nt2 nt3 ... nt

where all the values in the table above are integers and the rows are ordered in forward order
by the values of the indices n1, n2, n3, …., n.

For instance:

a) (49)
Is:

1 1 1
2 1 1
-1 1 0
1 1 1

b) (50)
Is:

307 2 2 2
4 2 1 1
1 4 -1 -1
-1 4 -1 1
-1 4 1 -1
1 4 1 1
 c)

(51)

Is

446 3 2 2 2
8 3 1 1 1
-1 8 -1 -1 -1
1 8 -1 -1 1
1 8 -1 1 -1
-1 8 -1 1 1
1 8 1 -1 -1
-1 8 1 -1 1
-1 8 1 1 -1
1 8 1 1 1

d) (52)
Is

458 1 1
5 4 4
1 1 -4
-4 1 -3
6 1 -2
-4 1 -1
1 1 0

8) Testing the Finite Difference Formulas Obtained and Checking the Error Order

It is observed that all finite difference formulas obtained can be numerically verified and tested,
and even the order of the error can be numerically validated by simply choosing an easily derived
function with an exponential, for example, u = e(1001/1000)x1e(1000/1003)x2e(1000/1005)x3... choose a value
for x1, x2, x3..., for example, one does x1=x2=x3=...=1.0, and then chooses a small value for the
step or increment of x1 or x2 or x3... ; for example, one does x1=0.01 or x2=0.01 or
x3=0.01... ; and then compares 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. One can 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 is selected, for example, do x1=0.001 or x2=0.001 or x3=0.001 or..., ten times smaller;
the value of the derivative is recalculated using the finite difference formulas and the new error
is obtained; with the values of the two errors with different steps, one can 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[24].
Thus, we have:

𝑆𝑡𝑒𝑝1 = 0.01 → 𝐸𝑟𝑟𝑜𝑟1 = |𝑢𝑒𝑥𝑎𝑐𝑡 − 𝑢𝑎𝑝𝑟𝑜𝑥1 |

𝑆𝑡𝑒𝑝2 = 0.001 → 𝐸𝑟𝑟𝑜𝑟2 = |𝑢𝑒𝑥𝑎𝑐𝑡 − 𝑢𝑎𝑝𝑟𝑜𝑥2 | (53)

𝐸𝑟𝑟𝑜𝑟 𝑂𝑟𝑑𝑒𝑟 = 𝑙 = 𝑟𝑜𝑢𝑛𝑑 (𝑙𝑜𝑔(𝑆𝑡𝑒𝑝1⁄ (𝐸𝑟𝑟𝑜𝑟1⁄𝐸𝑟𝑟𝑜𝑟2)) =

𝑆𝑡𝑒𝑝2)

𝑙𝑛(𝐸𝑟𝑟𝑜𝑟1⁄𝐸𝑟𝑟𝑜𝑟2)
𝑟𝑜𝑢𝑛𝑑 ( 𝑆𝑡𝑒𝑝1 ) → 𝑂(∆𝑥 𝑙 )
𝑙𝑛( ⁄𝑆𝑡𝑒𝑝2)

Special care must be taken so that Error1 and Error2 are nonzero so as not to produce an error
in the expression (53) 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[25]. 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[26,27].

9) Application in Theory of Disk Elasticity in Compression Using Electronic Spreadsheets

As the main subject of this work, the finite difference formulas obtained are used to solve
elasticity theory problems, with fourth-order derivatives, of a ring in compression by a constant
force, as shown in Figure 1, where the radius of the inner disk (R1) is 0.02 m, the radius of the
outer disk (R2) is 0.04 m, and the constant force (P) is 100 N, with an angle  of 10 degrees and
disc thickness of 0.0025 m[26,28]. 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.

The problem can be mathematically modeled by Airy's stress function () as:
𝜕4  𝜕4  𝜕4 
𝜕𝑥 4
+ 2 𝜕𝑥 2 𝜕𝑦2 + 𝜕𝑦4 = 0 (54)

In polar coordinates (r2 = x2 + y2,  = arctan(y/x)):

𝜕2 1 𝜕 1 𝜕2 𝜕2  1 𝜕 1 𝜕2 
( + + )( 2 + + ) =0
𝜕𝑟 2 𝑟 𝜕𝑟 2
𝑟 𝜕𝜃 2 𝜕𝑟 𝑟 𝜕𝑟 𝑟 2 𝜕𝜃2

Assume that the two-dimensional cartesian normal (x and y) and shear (xy) stresses
components can be represented by a stress function  such that:
𝜕2  𝜕2  𝜕2 
𝜎𝑥 = 𝜕𝑦2 𝜎𝑦 = 𝜕𝑥 2 𝜏𝑥𝑦 = − 𝜕𝑥𝜕𝑦 (55)
In polar coordinates, with radial stress r, circumferential stress  and polar shear stress r:
1 𝜕 1 𝜕2  𝜕2  1 𝜕 1 𝜕2 
𝜎𝑟 = 𝑟 𝜕𝑟 + 𝑟2 𝜕𝜃2 𝜎𝜃 = 𝜕𝑟2 𝜏𝑟𝜃 = 𝑟2 𝜕𝜃 − 𝑟 𝜕𝑟𝜕𝜃 (56)

The principal stresses yield:

𝜎𝑥 +𝜎𝑦 𝜎𝑥 −𝜎𝑦 2 𝜎𝑥 +𝜎𝑦 𝜎𝑥 −𝜎𝑦 2

𝜎1 = 2
+ √( 2
) + 𝜏𝑥𝑦 2 𝜎2 = 2
− √( 2
) + 𝜏𝑥𝑦 2 (57)

The maximum shear stresses:

(𝜎1 − 𝜎2 )⁄
2, 𝜎1 𝜎2 ≤ 0
𝜎1
𝜏𝑚𝑎𝑥 = ⁄2 , 𝜎1 > 0 𝑎𝑛𝑑 𝜎2 > 0 (58)
|𝜎2 |⁄
{ 2 𝜎1 < 0 𝑎𝑛𝑑 𝜎2 < 0

Expressions can be obtained that relate the polar stress components to the Cartesian stress
components as follows[29]:

𝜎𝑟 = 𝜎𝑥 𝑐𝑜𝑠 2 𝜃 + 𝜎𝑦 𝑠𝑖𝑛2 𝜃 + 𝜏𝑥𝑦 𝑠𝑖𝑛2𝜃 𝜎𝑥 = 𝜎𝑟 𝑐𝑜𝑠 2 𝜃 + 𝜎𝜃 𝑠𝑖𝑛2 𝜃 − 𝜏𝑟𝜃 𝑠𝑖𝑛2𝜃

{𝜎𝜃 = 𝜎𝑦 𝑐𝑜𝑠 2 𝜃 + 𝜎𝑥 𝑠𝑖𝑛2 𝜃 − 𝜏𝑥𝑦 𝑠𝑖𝑛2𝜃 or {𝜎𝑦 = 𝜎𝑟 𝑠𝑖𝑛2 𝜃 + 𝜎𝜃 𝑐𝑜𝑠 2 𝜃 + 𝜏𝑟𝜃 𝑠𝑖𝑛2𝜃 (59)
𝜏𝑟𝜃 = (𝜎𝑦 − 𝜎𝑥 )𝑠𝑖𝑛𝜃𝑐𝑜𝑠𝜃+𝜏𝑥𝑦 𝑐𝑜𝑠2𝜃 𝜏𝑥𝑦 = (𝜎𝑟 − 𝜎𝜃 )𝑠𝑖𝑛𝜃𝑐𝑜𝑠𝜃+𝜏𝑟𝜃 𝑐𝑜𝑠2𝜃

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

𝜎 = 𝜎𝑟 + 𝜎𝜃 (60)

The equivalent tensile stress or equivalent von Mises stress is defined as follows[29]:

(𝜎1 −𝜎2 )2 +(𝜎1 )2 +(𝜎2 )2

𝑣𝑜𝑛 = √ 2
(61)

The principal angle can be computed according to the formula (62):

1 (2𝜎 −𝜎)𝑠𝑖𝑛2𝜃+2𝜏 𝑐𝑜𝑠2𝜃
𝛽 = 2 𝑎𝑟𝑐𝑡𝑎𝑛 ((2𝜎𝑟−𝜎)𝑐𝑜𝑠2𝜃−2𝜏𝑟𝜃 𝑠𝑖𝑛2𝜃) (62)
𝑟 𝑟𝜃

The disk is subjected to in-plane uniform pressure that is distributed on the two diametrically
opposite regions of the external circular surface (boundary conditions)[29]:
𝜋 𝜋 3𝜋 3𝜋
−𝑃, 𝜃 ∈ [ 2 − 𝛼, 2 + 𝛼] ∪ [ 2 − 𝛼, 2 + 𝛼]
𝜎𝑟 (𝑅2, 𝜃) = { (63)
0, 𝑖𝑛 𝑜𝑡ℎ𝑒𝑟 𝑐𝑎𝑠𝑒𝑠
𝜋
𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒, 𝜎𝑟 (𝑅1, 𝜃) = 0, 𝜏𝑟𝜃 (𝑅1, 𝜃) = 0, 𝜏𝑟𝜃 (𝑅2, 𝜃) = 0, 𝜃 ∈ (0, ]
2
 Figure 1: Ring in compression by a constant force with an angle of 2

Tables 14 show the finite difference formulas utilized in solving the problem:
Table 14: Fourth-order approximations used in solving the elasticity problem.
The implicit method was used to solve the problem by finite differences, making r=0.001m,
=0.01rad and then comparing with r=0.0001m, =0.001rad 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 (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) was not used[1–3]. The
border of the mesh that was circular was approached by a series of very small rectangles that
formed a kind of ladder on the border, where each element is a step of the ladder, as shown in
Figure 2 (there are several techniques to improve this, but it was not necessary to use them
here).

Figure 2: Boundary of the spreadsheet with borders on each node in the form of rungs of a
ladder.

The computational solution of the problem was performed using electronic spreadsheets
(Microsoft Excel®), and the linear system was solved by an iterative method of execution of cells
in the spreadsheets. Figure 3 graphically shows the results, which were compared with the
analytical solutions and was successful up to the eighth decimal place of precision ( ≤ 10-8). To
gain computational efficiency in the spreadsheet, only the cells of the spreadsheet where x2+y2
are between the internal radius R12 (0.022) and the external radius R22 (0.04 2) were activated
for the calculation. Thus, the values of x are placed in the first row of the worksheet, and the
values of y are placed in the first column of the worksheet, and only the cells of the spreadsheet
are filled if the position of the cell is such that the value in the first row squared plus the value
of the first column squared (x2+y2) is between the inner radius R12 (0.022) and the outer radius
R22 (0.042). Each cell or node is a value of the incognita () Airy’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 given for the stress boundary conditions.

We chose to use an electronic spreadsheet to show one more of the advantages of the finite
difference method, which is the ease of computational implementation; it can be used with
simplicity even with electronic spreadsheets, without the need for programming. Thus, this is
an important option for the quick and simple numerical solution of differential equations. Few
numerical methods can be applied to solve practical engineering problems without having to
build software programs for their implementation.
 Radial Stress r (MPa) Circumferential Stress  (MPa) Polar Shear Stress r (MPa)

Principal Stress 1 (MPa) Principal Stress 2 (MPa) Maximum Shear Stress max (MPa)

Cartesian Stress x (MPa) Cartesian Stress y (MPa) Cartesian Shear Stress xy (MPa)

Von Mises Stress von (MPa) Principal Angle  (rad) In-plane Total Stress  (MPa)

Figure 3: Results obtained by finite differences.

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

Many other applications and tests of finite difference derivative approximations are provided in
the literature[4–7,24,30–92].

10) Implementation Computational Issues and Advantages

The finite difference method has a simpler mathematical formulation that recalls the beginnings
of differential and integral calculus, where the idea of limit was not yet well developed, and the
calculation of differentials was performed by numerical approximations using small additions.
Thus, the finite differences method is a powerful and efficient, high-precision numerical method
for the practical resolution of differential equations. In addition, the method has been gaining
new uses and new advantages with the evolution of computational technology, such as
processors with many high-speed cores that allow multiprocessing of tasks, machines with a
large amount of RAM (random-access memory) that allows many small steps and many nodes
and elements, graphical interfaces that facilitate geometric manipulation and visualization of
results and large databases that facilitate waiting and storage. The finite difference formulas and
the numerical results simplify the work with a significant volume of data and numbers.

One of the main advantages of finite difference methods is that in their explicit form, they do
not require the resolution of a linear system. Even in its implicit form, the system can often be
solved by iterative methods such as Gauss‒Seidel and relaxation and simply with successive and
repeated application of the formulas of finite differences obtained, as is done in many
spreadsheet programs with Microsoft Excel, enabling this option of iterative calculation of cells
in the worksheets. The need to solve large linear systems previously enabled the use of the finite
element method before the advent of computers with considerable RAM (random-access
memory), and it is still a problem in solving problems with numerous elements.

Another issue is that some of the formulas shown in this work have a very high error order,
which numerically makes the resolution very accurate, even if there is a large spacing of the
variables[93]. Steps and spacing can be used easily, much smaller than, for example, the pixels
of a high-resolution image, and with high speed of resolution in practical problems.

The use of processors with many cores, with large amounts of RAM and the development of
high-speed storage devices such as SSD (solid-state drive) and operating systems that
automatically swap between RAM and SSD, have significantly reduced the problem of the step
being constant (mesh of the same size) of the finite differences method. Meshes of straight or
rectangular segments or parallelepipeds with constant size, used in the finite difference method,
are faster and easier to construct or modify; much processing time is spent building complex
meshes in the finite element method. In addition, topological optimization techniques benefit
from the use of rectangular meshes or constant-sized parallelepipeds. Therefore, the advent of
new computational technologies favored all numerical methods, especially the finite difference
method.
Today, too, even with the finite element method, one has very small meshes in the entire
domain of solving a problem, similar to the finite difference method. Practical problems of
differential equations involving high-order derivatives or derivatives with long series and time
sequences also benefit from the use of finite differences.

11) Conclusion

Following 35 years of articles on the finite element method applied to structural problems, it
can be noted how formulas are deduced using computational algebra software, processors with
multiple cores, large amounts of RAM, graphical interfaces, and equation databases, made the
use and implementation of numerical methods applied to engineering simple and quick. This
work shows that, as modern microcomputers and a database with several finite difference
equations for derivative approximations, the vast majority of structural analysis problems can
be solved with high numerical precision by applying the finite difference method.

The 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.

Despite all the potentiality and applicability of the finite difference method, unfortunately, the
great limitation of the method persists: the need for the step to be constant (x1, x2, etc.).
However, it is believed that with the innovative approach and the line followed in this paper,
similar to the forward, backward and centered difference operators and the Gregory-Newton
and Stirling interpolation, one can create numerical operators of forward, backward and
centered divided differences and interpolation variant using these operators of divided
differences (variants of Newton's Formula for Interpolation). With this development of divided
differences, one would not be limited to the constant step in the independent variables.

The paper shows 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, one obtains 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.

As supplementary material of the paper, it follows many of the approximations of derivatives by

finite differences in PDF and text format, following the table model proposed in the paper to
store the finite difference equations.
Competing interests: The authors declare no competing interests.

Use of AI tools declaration: The authors declare they have not used Artificial Intelligence (AI)
tools in the creation of this article.

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.

References

[1] Hoffman JD, Hoffman JD, Frankel S. Numerical Methods for Engineers and Scientists.
2018. [Link]

[2] Steven CC. Applied Numercial Methods with MATLAB for Engineers and Scientist. vol.
53. 2013.

[3] Chapra SC, Canale RP. Numerical methods for engineers, second edition. Numerical
Methods for Engineers, Second Edition 2006.

[4] Dyatlovitsky LI, Cederbaum G, Jakobson L. Refined version of the finite differences
method in the theory of elasticity. Comput Struct 1995;57.
[Link]

[5] Papakaliatakis G, Simos TE. A finite difference method for the numerical solution of
fourth-order differential equations with engineering applications. Comput Struct
1997;65. [Link]

[6] Tuominen P, Jaako T. Generation of beam elements using the finite difference method.
Comput Struct 1992;44. [Link]

[7] Cocchi GM. Finite difference method with arbitrary grids in the elastic-static analysis of
three-dimensional continua. Comput Struct 2000;75. [Link]
7949(99)00087-5.

[8] Nayak GC. 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.
[Link]

[9] Magalhaes PAA, Junior PAAM, Magalhaes CA, Magalhaes ALMA. New Formulas of
Numerical Quadrature Using Spline Interpolation. Archives of Computational Methods
in Engineering 2020. [Link]

[10] Dukkipati R V. NUMERICAL METHODS. New Delhi: 2010.

[11] Steven C. Chapra. Applied Numerical Methods W/MATLAB: for Engineers & Scientists.
2011.

[12] Feagin T. High-order explicit Runge-Kutta methods using m-symmetry. Neural, Parallel
and Scientific Computations 2012;20.

[13] Verma AK, Dey SS. Integrated analysis of curved bridge superstructures by variational
finite difference method. Comput Struct 1991;38. [Link]
7949(91)90011-A.

[14] Ng SF, Bencharif N. A finite difference computer program for the modelling of thick
rectangular plates. Comput Struct 1989;33. [Link]
7949(89)90436-7.

[15] Ihlenburg FF. Plate bending analysis with variational finite difference methods on
general grid. Comput Struct 1993;48. [Link]

[16] Salehi M, Turvey GJ. Elastic large deflection response of annular sector plates-a
comparison of dr finite-difference, finite element and other numerical solutions.
Comput Struct 1991;40. [Link]

[17] Behravesh A, Davaran A, Kaveh A. A finite difference scheme with variable rectilinear
mesh for solving multi-harmonic partial differential equations. Comput Struct 1992;44.
[Link]

[18] Kaunda MAE. Forward-backward-difference time-integrating schemes with higher order

derivatives for non-linear finite element analysis of solids and structures. Comput Struct
2015;153. [Link]

[19] Arad M, Segev R, Ben-Dor G. Accuracy increase of finite difference calculations on

arbitrary meshes by means of differentiation of the partial differential equations and
their boundary conditions. Comput Struct 1997;64. [Link]
7949(96)00129-0.

[20] Lim KM, Li H. A coupled boundary element/finite difference method for fluid-structure
interaction with application to dynamic analysis of outer hair cells. Comput Struct
2007;85. [Link]

[21] Ahmed SR, Hossain MZ, Uddin MW. A general mathematical formulation for finite-
difference solution of mixed-boundary-value problems of anisotropic materials.
Comput Struct 2005;83. [Link]

[22] Jia X, Mang HA. Assessment of solutions from the consistently linearized eigenproblem
by means of finite difference approximations. Comput Struct 2015;151.
[Link]

[23] Zingoni A. A group-theoretic finite-difference formulation for plate eigenvalue

problems. Comput Struct 2012;112–113.
[Link]

[24] Lamon G, Magalhães Jr. P. Flow Measurement Ultrasonic Technology. SSRN Electronic
Journal 2022. [Link]
[25] HANDBOOK ON EXPERIMENTAL MECHANICS. Exp Tech 1989;13.
[Link]

[26] Tokovyy Y V., Hung KM, Ma CC. Determination of stresses and displacements in a thin
annular disk subjected to diametral compression. Journal of Mathematical Sciences
2010;165. [Link]

[27] Markides CF, Kourkoulis SK. Stresses and displacements in an elliptically perforated
circular disc under radial pressure. Engineering Transactions 2014;62.

[28] Vasil’ev V V., Fedorov L V. STRESS FUNCTIONS IN ELASTICITY THEORY. Mechanics of

Solids 2022;57. [Link]

[29] Sadd MH. Elasticity: Theory, Applications, and Numerics. 2020.

[Link]

[30] Magalhães PAA. New equations for phase evaluation in measurements with an
arbitrary but constant phase shift between captured intensity signs. Optical Engineering
2009;48. [Link]

[31] Cunha PPSR da, Souza PM, Valente LC, Mello GMV de, Junior PAAM. Analysis of induced
drag and vortex at the wing tip of a Blended Wing Body aircraft. International Journal of
Advanced Engineering Research and Science 2018;5.
[Link]

[32] Magalhaes PAA, Magalhaes CA. Higher-order Newton-Cotes formulas. J Math Stat
2010;6. [Link]

[33] Almeida Magalhaes C, Americo Almeida Magalhaes P. New numerical methods for the
photoelastic technique with high accuracy. J Appl Phys 2012;112.
[Link]

[34] Lamon G, Magalhães Jr. P, Carneiro J. Measurement of Sewage or Industrial Effluent.

SSRN Electronic Journal 2022. [Link]

[35] Lamon G, Magalhães Jr. P, Carneiro J. Cole Type Pitot Tube, Discharge Factor Survey
and Calibration. SSRN Electronic Journal 2022. [Link]

[36] Magalhaes PAA, Junior PAAM, Magalhaes CA, Magalhaes ALMA. New Formulas of
Numerical Quadrature Using Spline Interpolation. Archives of Computational Methods
in Engineering 2021;28. [Link]

[37] Almeida P de, Magalhães Jr PAA. Study of Inelastic Behavior of Reticular Structures of
Steel Bars - Phenomena of Plastification. International Journal of Advanced Engineering
Research and Science 2017;4. [Link]

[38] Salomão EA, Almeida PA, Junior M. Analysis of Turbocharged Engine Driven by Pulses
with Split Exhaust System and Distinct Discharge Valves. International Journal of
Advanced Engineering Research and Science 2018;5.
[Link]

[39] Duarte OC, Junior PAAM. Application of Cubic Spline Interpolation to Fit the Stress-
Strain Curve to SAE 1020 Steel. International Journal of Advanced Engineering Research
and Science 2017;4. [Link]
[40] Veloso V, Magalhães PAA, Landre J. A numerical model to study the effects of
aluminum foam filler on the dynamic behavior of a steel tubular energy absorber using
a multi-node displacement evaluation procedure. J Comput Nonlinear Dyn 2014;9.
[Link]

[41] Lima T, Andrade de W, Américo Almeida Magalhães P. Analysis of Stress in Nozzle/Shell

of Cylindrical Pressure Vessel under Internal Pressure and External Loads in Nozzle.
Journal of Engineering Research and Applications WwwIjeraCom ISSN 2015;5.

[42] Souza FR de, Júnior PAAM. Mathematical Methods applied in Image Enhancement
using Matlab. International Journal of Advanced Engineering Research and Science
2018;5. [Link]

[43] Magalhães ALMA, Magalhães Júnior PAA. Análise de uma Célula de Combustível de
Hidrogênio com um Sistema de Energia Fotovoltaica para um Drone. Proceedings of the
XXXVIII Iberian Latin American Congress on Computational Methods in Engineering,
2017. [Link]

[44] Magalhães ALMA, Magalhães Júnior PAA. NUMERICAL DEVELOPMENT OF INTEGRATION

RULES USING MAPLE SOFTWARE. Proceedings of the XXXVIII Iberian Latin American
Congress on Computational Methods in Engineering, 2017.
[Link]

[45] Magalhaes PAA, Neto PS, Magalhães CA. A modified carre algorithm for phase shifting
interferometry. Journal of Optics (India) 2010;39. [Link]
0001-1.

[46] Reis JPC dos, Júnior PAAM. Introduction to the Method of Finite Elements by a balance
Sheet Problem: A Simplification for an Initial understanding of the Method.
International Journal of Advanced Engineering Research and Science 2018;5.
[Link]

[47] Pedro Américo AM, Magalhães CA, Neto PS. New numerical method for the
photoelastic technique. Journal of the Brazilian Society of Mechanical Sciences and
Engineering 2012;34. [Link]

[48] REIS MNE, Oliveira W, Almeida Magalhães Júnior PA. INCOMPRESSIBLE FLOW AT LOW
REYNOLDS NUMBER IN A CONICAL ENTRANCE ORIFICE PLATE, 2019.
[Link]

[49] Alexandre da Silva Scari, Pedro Américo Almeida Magalhães Junior. Influence of the
Number of Rollers on a Tapered Roller Bearing. Journal of Mechanics Engineering and
Automation 2014;4. [Link]

[50] Magalhaes PAA, Neto PS, de Barcellos CS. Phase shifting technique using generalization
of Carre algorithm with many images. Opt Rev 2009;16.
[Link]

[51] Magalhaes Junior PAA, Neto PS, Magalhães CA. New Carré equation. Metrology and
Measurement Systems 2010;17. [Link]

[52] PIGNOLATI BD, Magalhães Júnior PAA, Soares LFC. STRUCTURAL ANALYSIS OF
MONOPLANE IPUC001-CARCARÁ WING BOX, UNDER ULTIMATE LOAD AND USING THE
 FINITE ELEMENTS METHOD. Proceedings of the XXXVIII Iberian Latin American Congress
on Computational Methods in Engineering, 2017.
[Link]

[53] Andrade TL, Júnior PAAM, Paula WA de. Tensile Test: Comparison Experimental,
Analytical and Numerical Methods. International Journal of Advanced Engineering
Research and Science 2017;4. [Link]

[54] Junior PAAM, Neto PS, De Barcellos CS. Generalization of Carré equation. Optik (Stuttg)
2011;122. [Link]

[55] Magalhães PAA, Magalhães CA, Magalhães ALMA. Computational methods of phase
shifting to stress measurement with photoelasticity using plane polariscope. Optik
(Stuttg) 2017;130. [Link]

[56] Ferreira TS, Magalhães PA, Moura FL, Ferreira TS. The effect of the cavity damping on
vehicular evaluation using the finite element method. Archives of Acoustics 2016;41.
[Link]

[57] Nascimento BC do, Carneiro JRG, Júnior PAAM. AVALIAÇÃO DE UM MODELO DE STENT
PELO MÉTODO DOS ELEMENTOS FINITOS. Proceedings of the XXXVI Iberian Latin
American Congress on Computational Methods in Engineering, 2015.
[Link]

[58] Oliveira WP de, Junior PAAM, Carneiro JRG. Adhesive and abrasive wear resistance
evaluation of a ABNT 52100 steel from a Stretch Bending Roll coated with WC-Co and
NICCr applied by the HVOF Process. International Journal of Advanced Engineering
Research and Science 2018;5. [Link]

[59] Filipe Gomes Vieira, Jordana Simões Ribeiro, Pedro Américo Almeida Magalhães Júnior.
New numerical methods for photoelastic parameters evaluation with plane
polariscope. Proceedings of the 23rd ABCM International Congress of Mechanical
Engineering, 2015. [Link]

[60] Magalhães CA, Neto PS, Almeida Magalhães PA, de Barcellos CS. Separation of
isochromatics and isoclinics phasemaps for the photoelastic technique with use phase
shifting and a large number of high precision images. Metrology and Measurement
Systems 2013;20. [Link]

[61] Chamone CSF, Teixeira DM, Silva MT, Magalhães Júnior PAA. Detector of Babies and
Animals Forgotten in a Vehicle with Remote Notification System. SAE Technical Papers,
vol. Part F127082, 2016. [Link]

[62] Scari A da S, Pockszevnicki BC, Landre Junior J, Magalhaes Junior PAA. Stress-Strain
Compression of AA6082-T6 Aluminum Alloy at Room Temperature. Journal of
Structures 2014;2014. [Link]

[63] Almeida Magalhães C, Smith Neto P, Almeida Magalhães Júnior PA, Sperb de Barcellos
C. Separation of Isochromatics and Isoclinics Phasemaps for the Photoelastic Technique
with use Phase Shifting and a Large Number of High Precision Images. Metrology and
Measurement Systems 2013;20. [Link]
[64] Siqueira Mazzaro R, de Morais Hanriot S, Jorge Amorim R, Américo Almeida Magalhães
Júnior P. Numerical analysis of the air flow in internal combustion engine intake ducts
using Herschel-Quincke tubes. Applied Acoustics 2020;165.
[Link]

[65] Marino L, Jakitsch B de O, Santos R, Magalhães Jr. PAA, Landre Jr. J. ANGLE CONTROL
OF A SIMPLE STAGE INVERTED PENDULUM: COMPARISON BETWEEN PID, PI-D AND
FUZZY. Proceedings of the XXXVIII Iberian Latin American Congress on Computational
Methods in Engineering, 2017. [Link]

[66] SANTOS LH, SOUSA WPT, DE DAFÉ SSF, JUNIOR PAM. Microstructural characterization
and mechanical behavior analysis of 7075-T6 aluminum subjected to simulated
lightning strikes. Chinese Journal of Aeronautics 2021;34.
[Link]

[67] Gouvea M de P, Teixeira ACR, Costa SCS, Júnior PAAM. Internal variation temperature
analysis and thermal mapping of a central processing unit (CPU). International Journal
of Advanced Engineering Research and Science 2018;5.
[Link]

[68] Vieira FG, Scari AS, Magalhães Júnior PAA, Martins JSR, Magalhães CA. Analysis of
stresses in a tapered roller bearing using three-dimensional photoelasticity and
stereolithography. Materials 2019;12. [Link]

[69] Ulhôa FCP, da Silveira HJJ, Hanriot S de M, Magalhães Júnior PAA. DEVELOPMENT OF A
VIRTUAL CI ENGINE MODEL FOR ELECTRIC TURBO COMPOUND APPLICATIONS.
Proceedings of the XXXVIII Iberian Latin American Congress on Computational Methods
in Engineering, 2017. [Link]

[70] Machado F da S, Campos TR, Duarte TP, Arrieta FRP, Júnior PAAM. Experimental
Determination of the Convective Coefficient of Heat Transfer Using the Global
Capacitance Method. International Journal of Advanced Engineering Research and
Science 2018;5. [Link]

[71] Filho A de F, Sales V, Júnior PAAM, Filho CT de C, Vieira AL. Analysis of Elastic Recovery
in The Process of Bending Sheets of Duplex Steel SAF 2205 via Experimental Method
and Numerical Simulation. International Journal of Advanced Engineering Research and
Science 2018;5. [Link]

[72] Sousa WPT, Duarte TP, Castro BR, Rego TFU, Júnior PAAM. Analytical Calculation of the
Drives of a Flight Simulator Platform with 2 Degrees of Freedom. International Journal
of Advanced Engineering Research and Science 2018;5.
[Link]

[73] Duarte TP, Sousa WPT, Castro BR, Rego TFU, Machado FS, Junior PAAM. Simulation
Mechanism with 2 Degrees of Freedom. International Journal of Advanced Engineering
Research and Science 2018;5. [Link]

[74] Caballero GM, De Carvalho Filho OA, Hargreaves BO, Brito HHDA, Magalhães PAA,
Oliveira DD. Mandibular canine intrusion with the segmented arch technique: A finite
element method study. American Journal of Orthodontics and Dentofacial Orthopedics
2015;147. [Link]
[75] Pratis SDB, Almeida MJPA, Sérgio MP, Silva MBF da, Amaral PG, Vieira BJA.
Productiveness Evaluation of a Machine Tool Manual Setup Compared with Automated
CNC Machine. International Journal of Advanced Engineering Research and Science
2018;5. [Link]

[76] Vieira FG, Martins JSR, Ribeiro JMS, Rocha MV dos S, Magalhães Júnior PAA, da Silva TA.
Numerical methods of phase shifting to stress measurement with digital photoelasticity
using plane polariscope. Proceedings of the XXXVIII Iberian Latin American Congress on
Computational Methods in Engineering, 2017.
[Link]

[77] de Brito GM, Brito HH de A, Marra GGM, Freitas LRP, Hargreaves BO, Magalhães PAA,
et al. Pure mandibular incisor intrusion: A finite element study to evaluate the
segmented arch technique. Materials 2019;12. [Link]

[78] Júnior PAAM, Vieira FG, Magalhães CA, Ribeiro JS, Rios IG. Numerical method to digital
photoelasticity using plane polariscope. Opt Express 2016;24.
[Link]

[79] Braga RA, Magalhaes PAA. Analysis of the mechanical and thermal properties of jute
and glass fiber as reinforcement epoxy hybrid composites. Materials Science and
Engineering C 2015;56. [Link]

[80] de Lana JA, Júnior PAAM, Magalhães CA, Magalhães ALMA, de Andrade Junior AC, de
Barros Ribeiro MS. Behavior study of prestressed concrete wind-turbine tower in
circular cross-section. Eng Struct 2021;227.
[Link]

[81] Junior PAAM, Portela TM. Evaluation of the Stiffness Effect of Pipe Supports in Relation
to Static and Dynamic Loads in a Flexibility Analysis. International Journal of Advanced
Engineering Research and Science 2018;5. [Link]

[82] Magalhaes PAA, Neto PS, De Barcellos CS. Analysis of shadow moire technique with
phase shifting using generalisation of carre method. Strain 2011;47.
[Link]

[83] Ives N, Pacheco R, De Castro D, Resende R, Américo P, Magalhães A. Stability Control of

an Autonomous Quadcopter through PID Control Law. Journal of Engineering Research
and Application WwwIjeraCom ISSN 2015;5.

[84] Moreira MF de S, Brito TP, Faria LPP, da Costa LA, Gomes HC, Ribeiro MS de B, et al.
Simulação computacional de um motor a combustão interna baseado no ciclo Otto para
identificação das perdas térmicas. Proceedings of the XXXVIII Iberian Latin American
Congress on Computational Methods in Engineering, 2017.
[Link]

[85] Wang B, Bai J, Lu S, Zuo W. Structural topology optimization considering geometrical

and load nonlinearities. Comput Struct 2023;289.
[Link]

[86] Shi C, Zhao C, Zhang X, Andersson A. Analysis on dynamic performance of different

track transition forms using the discrete element/finite difference hybrid method.
Comput Struct 2020;230. [Link]
[87] Lumba R, Datta A. Scalable mesh partitioning for multibody-3D finite element based
rotary-wing structures. Comput Struct 2023;281.
[Link]

[88] Ferreira BP, Carvalho Alves AF, Andrade Pires FM. An efficient finite strain constitutive
model for amorphous thermoplastics: Fully implicit computational implementation and
optimization-based parameter calibration. Comput Struct 2023;281.
[Link]

[89] Riesselmann J, Balzani D. A simple and efficient lagrange multiplier based mixed finite
element for gradient damage. Comput Struct 2023;281.
[Link]

[90] Ferradi MK, Fliscounakis A, Arquier M, Bleyer J. Elastoplastic and limit analysis of
reinforced concrete with an equilibrium-based finite element formulation. Comput
Struct 2023;286. [Link]

[91] Cortés F, Brun M, Elejabarrieta MJ. A finite element formulation for the transient
response of free layer damping plates including fractional derivatives. Comput Struct
2023;282. [Link]

[92] He T, Ma X. An edge-based smoothed finite element method for semi-implicit coupling

of unsteady viscoelastic fluid–structure interaction. Comput Struct 2023;281.
[Link]

[93] Martins PS, Almeida Magalhães Júnior PA, Gonçalves Carneiro JR, Talibouya Ba EC,
Vieira VF. Study of Diamond-Like Carbon coating application on carbide substrate for
cutting tools used in the drilling process of an Al–Si alloy at high cutting speeds. Wear
2022;498–499:204326. [Link]