Application of discretization error estimators in stepped column buckling problems

<!--==Application of Discretization Error Estimators in Stepped Column Buckling Problems==

<big>Beatriz F. Souza </big><sup>a</sup>*<big>, Daniel O. Fernandes <sup>a</sup>, Carla T. M. Anflor <sup>a</sup>, Marcus V. G. de Morais <sup>b</sup></big>

''<sup>a</sup> Post Graduate Program - Integrity of Engineering Materials, Group of Experimental and Computational Mechanics, University of Brasília, Brasília, Brazil; <sup>b </sup>Department of Mechanical Engineering, University of Brasília, Brasília, Brazil.''

* [mailto:bfesouza@gmail.com bfesouza@gmail.com] (''' [[Image:Draft_Souza_246201888-image1.jpeg|6px]] ''' 0000-0002-4357-8924)

[mailto:danieldeoliveirafernandes@gmail.com danieldeoliveirafernandes@gmail.com] ( [[Image:Draft_Souza_246201888-image1.jpeg|6px]] 0000-0002-7218-1942)

[mailto:anflor@unb.br anflor@unb.br] (''' [[Image:Draft_Souza_246201888-image1.jpeg|6px]] ''' 0000-0003-3941-8335)

[mailto:mvmorais@unb.br mvmorais@unb.br] ( [[Image:Draft_Souza_246201888-image1.jpeg|6px]] 0000-0003-4800-4513)

==<br/>Application of Discretization Error Estimators in Stepped Column Buckling Problems==-->
==Abstract==

In order to reduce the discretization error, in this paper, Richardson’s Extrapolation and Convergence Error Estimator were used to investigating the buckling problem convergence. The main objective was to verify the convergence order of the stepped column problem and to define a consistent moment of inertia at the point of variation of the cross-section. The variable of interest was the critical buckling load obtained by the Finite Difference Method. The convergent solution obtained errors less than 10<sup>-8</sup>, and this work showed that the best solution is not defined by excessive mesh refinement, but by the solution convergence analysis.

'''Keywords''': Finite Difference Method, discretization error, convergence order, buckling analysis, Finite Element Method

==1. Introduction==

The use of techniques to reduce the numerical error is increasingly widespread, mainly due to the decrease in computational expenditure during the simulations. The reduction of the numerical error caused by discretization errors can be obtained in several ways: the mesh refinement and the increase in the precision order of the numerical solutions are some of these alternatives. However, there are disadvantages such as high computational cost and high level of complexity of the suggested model. In order to avoid such inconveniences, the use of extrapolation methods, such as Richardson Extrapolation (RE) [1] are very useful and applicable.

When RE is used successively, the technique is called Repeated Richardson Extrapolation (RRE). In this case, the numerical solution of the variable of interest must be defined in at least three different meshes, thus allowing two or more extrapolations. According to [2] there are many advantages when using the RRE technique that can be listed as: (a) Decreased of the discretization error; (b) Easy application and post-processing method, thus it does not interfere with the numerical solution obtained in a given mesh; (c) The computational cost (CPU time and RAM) is considered low; (d) Can be applied to existing codes or results obtained; (e) For different numerical methods and variables of interest is suitable; (f) It does not depend on priori analysis or the analytical solution of the problem; (g) Can be used in low-order numerical solutions to obtain higher-order numerical solutions.

The most interesting fact is that the technique allows defining a convergent numerical solution and its error estimators. In this sense, the true error and the estimated error of a numerical solution can be reliably reduced. In addition, the true discretization error is limited by the Richardson Error estimator [2–5].

Although the method is widely used in isotropic meshes, RE has also shown good results in improving the convergence rate (from fourth-order to sixth-order) in models of non-isotropic meshes using a finite difference compact method [6]. Furthermore, one worth noting that the convergence rates are strongly dependent on the order of the finite discrete scheme: finite schemes of a higher order lead to enriched constitutive laws of a higher order, with a higher convergence rate [7,8].

The variable of interest of this paper is the critical load of buckling in stepped columns simply supported. Columns of uniform cross section subjected to compressive loading may not be economical, so stability can be increased by removing a portion of material from the ends and increasing the cross section in the central region of the column [9]. Although the problem seems simple, several hypotheses of the stiffness in the coincident node between different cross sections can be approached, also, the convergence of the problem and the resolution method is relevant to be evaluated. According to the literature [10–13], the resolution of the problem by the Finite Differences Method (FDM) proved to be quite appropriate. Multiple-stepped beams are frequently discussed in literature to describe free vibration [14–16] or buckling behavior [17–20].

Despite the successful use of the RE method, [21] warns that the application of FDM in the problem of stepped column presents some pitfalls: (a) the convergence pattern can change when the column is divided into an even or odd number of segments and (b) formulations using a fourth-order differential can lead to incorrect solutions. Although there is a need for caution in the use of FDM, for the proposed case, precise solutions for buckling load can be obtained using a small number of discretization divisions, especially if extrapolation techniques are used.

The purpose of this work is to verify the estimate of the discretization error when the apparent error order is monotonic converging to the asymptotic order of the error. Particularly, in this work assumed: (a) the FDM in solving a one-dimensional problem of stepped columns simply supported with isotropic grids, (b) the posteriori error estimation, (c) the RRE is used to decrease and estimate the discretization errors.

The outline of this paper is as follows. The exact solution for the stepped column buckling problem is introduced in section 2, subsequently, section 3 explains the FDM for discontinued section columns. The approach of the convergence analysis and error estimators is introduced in section 4. The problem description and numerical procedure are showed in section 5. Section 6 presents the numerical results obtained for stepped columns with different stiffness. Finally, section 7 presents the final considerations.

==2. Timoshenko’s solution buckling column with section discontinuity==

Considering a clamped-free column ([[#_Ref45029125|Figure 1]]a) or the simply supported column ([[#_Ref45029125|Figure 1]]b) and using the procedures developed in [9], it can be possible to define the exact solution of the stepped column problem.

<span id="_Ref45029125"></span>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;max-width: auto;"
|-
|style="padding:10px;"| [[Image:Draft_Souza_246201888-image2-c.png|center|318px]]
|style="padding:10px;"| [[File:Figura 1b.png|centre|301x301px]]
|-
|  style="text-align: center;font-size: 75%;"|(a) 
|  style="text-align: center;font-size: 75%;"|(b) 
|- style="text-align: center; font-size: 75%;"
| colspan="2" style="padding:10px;"| '''Figure 1'''. Stepped Column. (a) clamped-free column. (b) Pinned-pinned column
|}


In the stepped columns illustrated in <span id='cite-_Ref45029125'></span>[[#_Ref45029125|Figure 1]]a and <span id='cite-_Ref45029125'></span>[[#_Ref45029125|Figure 1]]b, the bending rigidity <math display="inline">(EI)</math> changes suddenly; then, it is necessary to analyze each section separately. The governing differential equation in a clamped-free stepped column (<span id='cite-_Ref45029125'></span>[[#_Ref45029125|Figure 1]]a), for section 1 of length <math display="inline">{L}_{1}</math>, can be written as

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="margin:auto;width: 100%;"
|-
| <math>\frac{{d}^{2}{w}_{1}}{d{x}^{2}}=\frac{P}{E{I}_{1}}\left( \delta -{w}_{1}\right) .</math>
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|<span id='_Ref45031424'></span>(1)
|}


Similarly, for the section 2, the differential equation can be expressed by

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="margin:auto;width: 100%;"
|-
| <math>\frac{{d}^{2}{w}_{2}}{d{x}^{2}}=\frac{P}{E{I}_{2}}\left( \delta -{w}_{2}\right) ,</math>
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|<span id='_Ref45031426'></span>(2)
|}

where <math display="inline">E{I}_{i}</math> is the bending rigidity in section <math display="inline">i</math>, <math display="inline">P</math> is the axial compression, <math display="inline">\delta</math>  is the deflection in the top of the column and <math display="inline">{w}_{i}</math> is the lateral displacement throughout section <math display="inline">i</math>.

The general solution of the differential equations is known and given by

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center;margin:auto;width: 100%;"
|-
| <math display="inline">{w}_{1}=A\cos\left( {\psi }_{1}x\right) +B\sin\left( {\psi }_{1}x\right) +</math><math>\delta</math>,
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|<span id='_Ref45032138'></span>(3)
|}

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center;margin:auto;width: 100%;"
|-
| <math display="inline">{w}_{2}=C\cos\left( {\psi }_{2}x\right) +D\sin\left( {\psi }_{2}x\right) +</math><math>\delta</math>,
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|<span id='_Ref45113322'></span>(4)
|}

being <math display="inline">A</math>, <math display="inline">B</math>, <math display="inline">C</math> and <math display="inline">D\,</math> constants, and

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="margin:auto;width: 100%;"
|-
| <math>{\psi }_{i}=\sqrt{\frac{P}{E{I}_{i}}}.</math>
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|(5)
|}


Applying the boundary conditions,

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="margin:auto;width: 100%;"
|-
| <math>{\left( {w}_{2}\right) }_{x=0}=0,\quad \quad \quad \quad \quad {\left( \frac{{dw}_{2}}{dx}\right) }_{x=0}=</math><math>0.</math>
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|(6)
|}


From <span id='cite-_Ref45113322'></span>[[#_Ref45113322|(4)]], the constants <math display="inline">C</math> and <math display="inline">\, D\,</math> are defined by

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="margin:auto;width: 100%;"
|-
| <math>C=-\delta ,\quad \quad \quad \quad \quad \quad \quad \quad \, \, D=0.</math>
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|(7)
|}

Then, <span id='cite-_Ref45113322'></span>[[#_Ref45113322|(4)]] results as

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center;margin:auto;width: 100%;"
|-
| <math display="inline">{w}_{2}=\delta [1-\cos\,\left( {\psi }_{2}x\right) ]</math>.
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|<span id='_Ref45032004'></span>(8)
|}

One known that deflection is <math display="inline">\delta</math> for <math display="inline">x=L</math> and the same for <math display="inline">x={L}_{2}</math>, one can be written

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="margin:auto;width: 100%;"
|-
| <math>\delta =A\cos\left( {\psi }_{1}L\right) +B\sin\left( {\psi }_{1}L\right) +\delta</math> 
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|<span id='_Ref45032502'></span>(9)
|}

and

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center;margin:auto;width: 100%;"
|-
| <math display="inline">\delta \left[ 1-\cos\,\left( {\psi }_{2}{L}_{2}\right) \right] =\, A\cos\left( {\psi }_{1}{L}_{2}\right) +B\sin\left( {\psi }_{1}{L}_{2}\right) +\delta</math>.
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|<span id='_Ref45032506'></span>(10)
|}

From <span id='cite-_Ref45032502'></span>[[#_Ref45032502|(9)]] and <span id='cite-_Ref45032506'></span>[[#_Ref45032506|(10)]] can determine the constants,

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="margin:auto;width: 100%;"
|-
| <math>A=-B\tan\left( {\psi }_{1}L\right) ,\quad \quad \quad \quad \quad \, B=\frac{\delta \cos\,\left( {\psi }_{2}{L}_{2}\right) \mathrm{cos}\,\left( {\psi }_{1}L\right) }{\sin\,{\psi }_{1}{L}_{1}}.</math>
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|<span id='_Ref45033005'></span>(11)
|}

Since the tangent is the same to <math display="inline">x={L}_{2}</math> then,

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center;margin:auto;width: 100%;"
|-
| <math display="inline">\delta {\psi }_{2}\sin\,\left( {\psi }_{2}{L}_{2}\right) =-A{\psi }_{1}\mathrm{sin}\,\left( {\psi }_{1}{L}_{2}\right) +B{\psi }_{1}\cos\,({\psi }_{1}{L}_{2})</math>.
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|<span id='_Ref45033007'></span>(12)
|}

Substituting <span id='cite-_Ref45033005'></span>[[#_Ref45033005|(11)]] in <span id='cite-_Ref45033007'></span>[[#_Ref45033007|(12)]] get the transcendental equation:

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="margin:auto;width: 100%;"
|-
| <math>\tan\left( {\psi }_{1}{L}_{1}\right) \tan\,\left( {\psi }_{2}{L}_{2}\right) =\frac{{\psi }_{1}}{{\psi }_{2}}.</math>
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|(13)
|}


Then, can write the equation in a more explicit form:

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="margin:auto;width: 100%;"
|-
| <math>\tan\left( \sqrt{\frac{P}{E{I}_{1}}}{L}_{1}\right) \tan\,\left( \sqrt{\frac{P}{E{I}_{2}}}{L}_{2}\right) =\sqrt{\frac{{I}_{2}}{{I}_{1}}}.</math>
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(14)
|}

Substituting <math display="inline">{L}_{1}</math> for <math display="inline">(L-</math><math>a)/2</math> and <math display="inline">{L}_{2}</math> for <math display="inline">a/2</math>, the solution for the case of the <span id='cite-_Ref45029125'></span>[[#_Ref45029125|Figure 1]]b is obtained as

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="margin:auto;width: 100%;"
|-
| <math>tan\left( \sqrt{\frac{P}{E{I}_{1}}}\frac{(L-a)}{2}\, \, \right) \tan\,\left( \sqrt{\frac{P}{E{I}_{2}\, }}\frac{a}{2}\, \, \right) =\sqrt{\frac{{I}_{2}}{{I}_{1}}}.</math>
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|<span id='eq15'></span><span id='_Ref45033979'></span>(15)
|}

The load <math display="inline">P</math> that satisfies the transcendental equation is called the buckling critical load. The solution of this equation can be obtained by numerical methods, for example, the Newton-Raphson Method. Further details for several cases can be achieved in classic literature [9].

==3. Finite difference applied to buckling column with section discontinuity==

The simply supported column (<span id='cite-_Ref47016448'></span>[[#_Ref47016448|Figure 2]]a) consists of three segments proportional to bending rigidities <math display="inline">EI</math>. Focusing only on the first buckling mode, it is sufficient to model a half column, as described in <span id='cite-_Ref47016448'></span>[[#_Ref47016448|Figure 2]]b.

<span id='_Ref47016448'></span>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 58%;"
|-
|style="padding:10px;"|  [[File:Figura 2a.png|centre|479x479px]]
|-
|  style="text-align: center;font-size: 75%;"|(a) 
|-
|style="padding:10px;"|  [[File:Figura 2b.png|centre|455x455px]]
|-
|  style="text-align: center;font-size: 75%;"|(b) 
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding:10px;"| '''Figure 2'''. Simplified finite difference model representation of stepped column. (a) Geometric description. (b) Finite difference discretization (<math display="inline">m=n/2)</math>
|}


The differential equation for lateral displacement <math display="inline">w(x)</math> is stated as

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="margin:auto;width: 100%;"
|-
| <math>\frac{d}{dx}\left( EI\left( x\right) \frac{dw\left( x\right) }{dx}\right) +P\, w\left( x\right) =0,</math>
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|<span id='_Ref32582270'></span>(16)
|}

where <math display="inline">EI(x)</math> is the bending rigidity as a function of coordinate <math display="inline">x</math>, <math display="inline">P</math> is the axial compressive load. The boundary conditions are described by

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="margin:auto;width: 100%;"
|-
| <math>w\left( 0\right) ={w}^{'}\left( 0\right) =0,</math>
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|(17)
|}

where <math display="inline">{w}^{'}\left( x\right) =dw\left( x\right) /dx</math>.

Applying the FDM, we replace the derivatives of <math display="inline">w\left( x\right)</math>  by central finite difference of error order <math display="inline">O\left( \Delta {x}^{2}\right)</math>  at the interior nodes <math display="inline">i</math> (nodes <math display="inline">1</math> to <math display="inline">n</math>):

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center;margin:auto;width: 100%;"
|-
| <math>EI\left( x\right) \left( {w}_{i-1}-2{w}_{i}+{w}_{i+1}\right) +P\Delta {x}^{2}\, {w}_{i}=</math><math>0,</math>
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|<span id='_Ref47264100'></span>(18)
|}

where <math display="inline">{w}_{i}=w\left( i\, \Delta x\right)</math>, for <math display="inline">i\in 1</math> to <math display="inline">n</math>.

At point <math display="inline">m</math> ([[#_Ref47016448|Figure 2]]b), we have difficulty determining flexural inertia <math display="inline">EI\left( x\right)</math>  for finite difference operator. As related by [12] and [21], the numerical free vibration analysis or linear buckling load determination of stepped beams using finite difference formulation need special attention at step discontinuity. If flexural rigidity average is chosen, we observe an increase of numerical error of the finite difference model. This approach does not respect the continuity conditions of bending moment for each segment. In literature, a consistent approach to discretizing the stepped beam joint by an equivalent moment of inertia is proposed as follows.

For example, to applying the differential operator <span id='cite-_Ref47264100'></span>[[#_Ref47264100|(18)]] at step discontinuity (point <math display="inline">m</math>), one proposed creation of two ghost points: (a) one ghost point at the right of the thin beam (with flexural rigidity <math display="inline">E{I}_{L}</math>), and (b) one ghost at the left of the thick beam (with flexural rigidity <math display="inline">E{I}_{R}</math>). Both differential operators are described, to thin beam:

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center;margin:auto;width: 100%;"
|-
| <math>E{I}_{L}\left( {w}_{i-1}-2{w}_{i}+{w}_{a}\right) +P\Delta {x}^{2}\, {w}_{i}=0</math>
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|<span id='_Ref24031761'></span>(19)
|}

and to thick beam:

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center;margin:auto;width: 100%;"
|-
| <math display="inline">E{I}_{R}\left( {w}_{b}-2{w}_{{i}^{'}}+{w}_{{i}^{'}+1}\right) +P\Delta {x}^{2}\, {w}_{{i}^{'}}=0</math>.
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|<span id='_Ref24031763'></span>(20)
|}

To assure the displacement and bending moment continuity, respectively, we applying the following boundary conditions:

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center;margin:auto;width: 100%;"
|-
| <math>{w}_{i}={w}_{{i}^{'}}\, ,</math>
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|<span id='_Ref24031769'></span>(21)
|-
|}

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
| 
{| style="text-align: center;margin:auto;width: 100%;"
|-
| <math>\frac{\partial {w}_{i}}{\partial x}=\frac{\partial {w}_{{i}^{'}}}{\partial x}\, \Rightarrow {w}_{i-1}+{w}_{i+1}={w}_{b}+{w}_{a}.</math>
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|<span id='_Ref24031770'></span>(22)
|}

Adding the expressions <span id='cite-_Ref24031761'></span>[[#_Ref24031761|(19)]] and <span id='cite-_Ref24031763'></span>[[#_Ref24031763|(20)]] applying the boundary conditions <span id='cite-_Ref24031769'></span>[[#_Ref24031769|(21)]] and <span id='cite-_Ref24031770'></span>[[#_Ref24031770|(22)]], we have:

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="margin:auto;width: 100%;"
|-
| <math>\left( {w}_{i-1}-2{w}_{i}+{w}_{i+1}\right) +\frac{P\Delta {x}^{2}}{E{I}_{eq}}\, {w}_{i}=0,</math>
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|(23)
|}

where the equivalent moment of inertia <math display="inline">{I}_{eq}</math> is expressed as

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="margin:auto;width: 100%;"
|-
| <math>{I}_{eq}=\frac{2{I}_{L}{I}_{R}}{{I}_{L}+{I}_{R}}.</math>
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|(24)
|}

In literature, there are other non-consistent equivalent moments of inertia <math display="inline">{I}_{eq}</math> possibilities described as

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="margin:auto;width: 100%;"
|-
| <math>{I}_{eq}\simeq {I}_{L},</math>
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|(25)
|}

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="margin:auto;width: 100%;"
|-
| <math>{I}_{eq}\simeq {I}_{R},</math>
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|(26)
|}

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="margin:auto;width: 100%;"
|-
| <math>{I}_{eq}\simeq {I}_{avg}=\frac{1}{2}\left( {I}_{L}+{I}_{R}\right) .</math>
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|(27)
|}

This work proposes to analyse each one of these moments of inertia (24)-(27) at the point <math display="inline">m\,</math> of the intersection of the stepped column, in this way, it will be possible to verify the consistency in the numerical modelling in each case.

==4. Convergence analysis and grid convergence index==

A reliable numerical solution should have a null error or very close to zero. In order to reduce the numerical error caused by discretization problems, the Richardson Error Estimator and Convergent Error Estimator were used to verify the convergence for FDM approaches. To perform the convergence methods analysis, one necessary to define 3 different meshes: coarse <math display="inline">({\phi }_{3})</math>, medium <math display="inline">({\phi }_{2})</math> and fine <math display="inline">{(\phi }_{1})</math>, with grades <math display="inline">\Delta {h}_{3}</math>, <math display="inline">\Delta {h}_{2}</math> e <math display="inline">\Delta {h}_{1}</math>, respectively.

According to [4], [5] the error estimation could be performed following the steps:

===4.1 Richardson error estimator===

The true numerical error is determined through the equation [[#eq28|(28)]] where <math display="inline">\Phi</math>, is the exact analytical solution and <math display="inline">\phi</math>  is the proposed numerical solution:

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="margin:auto;width: 100%;"
|-
| <math>\epsilon \, \left( {\phi }_{1}\right) =\Phi -{\phi }_{1}.</math>
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|<span id='eq28'></span><span id='_Ref26791192'></span>(28)
|}

The estimative of the analytical solution <math display="inline">(\Phi )</math> could be performed through

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="margin:auto;width: 100%;"
|-
| <math>{\phi }_{\infty }\left( PL\right) ={\phi }_{1}+\left( \frac{{\phi }_{1}-{\phi }_{2}}{{r}^{PL}-1}\right) ,</math>
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|<span id='_Ref26783724'></span>(29)
|}

where <math display="inline">r</math> is the grid refinement ratio:

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="margin:auto;width: 100%;"
|-
| <math>r=\frac{\Delta {h}_{2}}{\Delta {h}_{1}}=\, \frac{\Delta {h}_{3}}{\Delta {h}_{2}}</math>
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|<span id='_Ref26791973'></span>(30)
|}

and <math display="inline">PL</math> is the asymptotic order of the discretization error.

The apparent order of the uncertainty (<math display="inline">PU</math>) is calculated through

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="margin:auto;width: 100%;"
|-
| <math>PU=\frac{\mathrm{log}\,\left( \frac{{\phi }_{2}-{\phi }_{3}}{{\phi }_{1}-{\phi }_{2}}\right) }{\mathrm{log}\,(r)}</math>
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|<span id='_Ref26785381'></span>(31)
|}

and the estimative of the numerical solution <math display="inline">{\phi }_{\infty }\left( PU\right)</math>  could be obtained by

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="margin:auto;width: 100%;"
|-
| <math>{\phi }_{\infty }\left( PU\right) ={\phi }_{1}+\left( \frac{{\phi }_{1}-{\phi }_{2}}{{r}^{PU}-1}\right) .</math>
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|<span id='_Ref26785901'></span>(32)
|}

If the apparent order of the error is monotonic convergent, so the analytical solution <math display="inline">(\Phi )</math> is between <math display="inline">{\phi }_{\infty }\left( PL\right)</math>  and <math display="inline">{\phi }_{\infty }\left( PU\right)</math>.

Equations <span id='cite-_Ref26783724'></span>[[#_Ref26783724|(29)]] and <span id='cite-_Ref26785901'></span>[[#_Ref26785901|(32)]] are generalized Richardson Extrapolation. The error estimative of numerical solution ( <math display="inline">{\phi }_{1}</math>) could be given by the Richardson Error Estimator (<math display="inline">{U}_{Ri}</math>), calculated by

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="margin:auto;width: 100%;"
|-
| <math>{U}_{Ri}\left( {\phi }_{1}\right) =\, sg\left( {\phi }_{1}-{\phi }_{2}\right) \mathrm{max}\,\lbrace \left| {U}_{Ri}\left( {\phi }_{1},PL\right) \right| ;{\vert U}_{Ri}\left( {\phi }_{1},PU\right) \vert \rbrace ,</math>
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|<span id='_Ref26790570'></span>(33)
|}

been the maximum absolute value between the numerical solution error estimator, obtained by asymptotic error order is

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="margin:auto;width: 100%;"
|-
| <math>{U}_{Ri}\left( {\phi }_{1},PL\right) =\left| \frac{{\phi }_{1}-{\phi }_{2}}{{r}^{PL}-1}\right|</math> 
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|(34)
|}

and by apparent error order is

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="margin:auto;width: 100%;"
|-
| <math>{U}_{Ri}\left( {\phi }_{1},PU\right) =\left| \frac{{\phi }_{1}-{\phi }_{2}}{{r}^{PU}-1}\right| .</math>
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|<span id='_Ref26790907'></span>(35)
|}

One can be said the estimated error <math display="inline">U</math> is reliable if the condition founded

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="margin:auto;width: 100%;"
|-
| <math>\left| \frac{{U}_{Ri}\left( {\phi }_{1}\right) }{\epsilon \, ({\phi }_{1})}\right| \geq 1</math>
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|<span id='_Ref26791531'></span>(36)
|}

is fulfilled. In other words, the ratio between the estimated error <math display="inline">U</math> and the true error <math display="inline">\epsilon \,</math>  is greater or equal to the unit.

If the aim is obtaining a reliably estimated error, the numerical solution of the variable could be calculated through

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="margin:auto;width: 100%;"
|-
| <math>\phi ={\phi }_{1}+{U}_{Ri}\left( {\phi }_{1}\right) .</math>
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|<span id='_Ref26791905'></span>(37)
|}

<span id='_Ref24034320'></span><span id='_Ref32577085'></span>

===4.2 Convergent error estimator===

Using the same numerical solutions (<math display="inline">{\phi }_{1},\, {\phi }_{2},\, {\phi }_{3}</math>), is possible to calculate an estimative of the convergent numerical solution that allows minimizing the true discretization error. The numerical convergent solution can be calculated by

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="margin:auto;width: 100%;"
|-
| <math>{\phi }_{c}=\frac{{\phi }_{\infty }\left( PL\right) +{\phi }_{\infty }\left( PU\right) }{2},</math>
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|<span id='_Ref26792239'></span>(38)
|}

where <math display="inline">{\phi }_{\infty }\left( PL\right)</math>  and <math display="inline">{\phi }_{\infty }\left( PU\right)</math>  have already been determinate in equations <span id='cite-_Ref26783724'></span>[[#_Ref26783724|(29)]] and <span id='cite-_Ref26785901'></span>[[#_Ref26785901|(32)]] by Richardson extrapolation.

The convergent error estimative of the convergent numerical solution could be determined by:

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="margin:auto;width: 100%;"
|-
| <math>{U}_{c}\left( {\phi }_{c}\right) =\left| \frac{{\phi }_{\infty }\left( PL\right) -{\phi }_{\infty }\left( PU\right) }{2}\right|</math> 
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|<span id='_Ref26792678'></span>(39)
|}

and the true discretization error is calculated through

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="margin:auto;width: 100%;"
|-
| <math>\epsilon \, \left( {\phi }_{c}\right) =\Phi -{\phi }_{c}.</math>
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|<span id='_Ref26792981'></span>(40)
|}

In this case, the estimated error <math display="inline">{U}_{c}</math> is reliable if satisfy what is proposed by

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="margin:auto;width: 100%;"
|-
| <math>\left| \frac{{U}_{c}\left( {\phi }_{c}\right) }{\epsilon \, ({\phi }_{c})}\right| \geq 1.</math>
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|<span id='_Ref26793071'></span>(41)
|}

The numerical solution of the analyzed variable can be defined by

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="margin:auto;width: 100%;"
|-
| <math>\phi ={\phi }_{c}\pm {U}_{c}\left( {\phi }_{c}\right) .</math>
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|<span id='_Ref26793254'></span>(42)
|}

One important to highlight, according to Marchi and Silva [4,5], in a convergent interval of <math display="inline">Pu</math>, is recommended using a numerical convergent solution <math display="inline">{\phi }_{c}</math> instead using <math display="inline">{\phi }_{1}</math>, so the true error <math display="inline">\epsilon \, ({\phi }_{c})</math> is smaller than <math display="inline">\epsilon \, ( {\phi }_{1})</math>.

==5. Problem Description and Numerical Procedure==

This work analyzed the convergence of the critical buckling load in a stepped column using FDM. The cross-section of the column was square and the length (<math display="inline">L</math>) was equal to <math display="inline">100mm</math>, being composed of different cross-sections in which the intermediate section (section 2) has a moment of inertia greater than the extremity sections ([[#_Ref47016448|Figure 2]]). The material used had Young's modulus equal to <math display="inline">200GPa\,</math> and Poisson's coefficient equal to <math display="inline">0.3</math>.

The moment of inertia (<math display="inline">{I}_{L}</math>) to the left of point <math display="inline">m</math> was kept constant and equal to  <math display="inline">0.5m{m}^{4}</math>. To validate the methodology, two cases were investigated. '''CASE 1''' considering the ratio between cross-section moments of inertia (<math display="inline">{I}_{R}</math>/ <math display="inline">{\, I}_{L}</math>) equal to <math display="inline">1.5</math> and '''CASE 2''' equal to <math display="inline">2</math>. Moreover, four different moments of inertia (<math display="inline">{I}_{L},\, {I}_{R},\, {I}_{avg},\, {I}_{eq})\,</math> - discussed in section 3 - were considered in node <math display="inline">m</math> of the intersection (<span id='cite-_Ref47016448'></span>[[#_Ref47016448|Figure 2]]b) and the consistency of the modeling was verified in each case.

The FDM implementation was developed in MATLAB® and for validation, one was also simulated in ANSYS APDL® by Finite Element Method (FEM). Initially, Newton Raphson's method was used to obtain the exact solution of the transcendental equation for the buckling problem. After the numerical simulations (FDM and FEM), the data post-processing were used for the convergence analysis (discussed in section 4). Finally, the asymptotic behavior of the discretization error and possible numerical instabilities were investigated in each case of the problem proposed.
<span id='_Ref32848410'></span>

==6. Numerical Results==

This section presents a comparison between the discretization error by FEM using commercial software (ANSYS®) and FDM, posteriorly the results by FDM are showed for each moment of inertia considered in this study. Then, the analysis of the results allows verifying the convergence of the proposed problem.

===6.1 Buckling of Stepped Column===

<span id='cite-_Ref21638069'></span>[[#_Ref21638069|Figure 3]] shows a comparison between the FEM and FDM. One can be noted that convergence by FEM is the 4<sup>th</sup> order and by FDM is the 2<sup>nd</sup> order. This occurs because the approximation for the finite difference operator is 2<sup>nd</sup> order and the deflection approximation used in ANSYS® is 4<sup>th</sup> order.

{| style="width: 100%;margin: 1em auto 0.1em auto;border-collapse: collapse;" 
|-
|  style="text-align: center;width: 100%;"|[[Image:Draft_Souza_246201888-image8-c.png|372px]] 
|-
|  style="text-align: center;"|<span id='_Ref21638069'></span>Figure 3 – Convergence analysis by FEM and FDM for stepped column (CASE 1 and CASE 2)
|}


From <span id='cite-_Ref21638069'></span>[[#_Ref21638069|Figure 3]], one can be seen that the two numerical models showed monotonic convergence until reaching the computational limit of approximately <math display="inline">\, {10}^{-6}</math>. Moreover, it is clearly visible that the FDM needs further refinement of the mesh compared to the FEM. For example, to obtain an error of <math display="inline">{10}^{-4}</math>, the FDM needed approximately 1100 elements while the FEM used only 20 elements. Nevertheless, the use of extrapolation methods can cause the FDM to reach the order of convergence of the FEM from the determination of the convergent solution. This will consequently lead the FDM to achieve much smaller errors than those presented by the FEM (see section 6.2).

In both cases, there is a numerical instability is generated by the computational limit. In this case, it is advisable to always choose to use solutions outside the numerical instability zone, in addition to checking the error considered reliable for each type of problem. Finally, the convergent monotonic order is considered consistent according to the numerical model proposed by FEM and FDM.

===6.2 Convergence Analysis of Inertia Models by Finite Difference Method===

The discretization true error <math display="inline">(\epsilon )</math> and the estimated error <math display="inline">(U)</math> for CASE 1<span style="text-align: center; font-size: 75%;"> and CASE 2 </span>are shown in <span id='cite-_Ref72141256'></span>[[#_Ref72141256|Figure 4]] and <span id='cite-_Ref72141259'></span>[[#_Ref72141259|Figure 5]], respectively. The results considered four different moments of inertia ( <math display="inline">{I}_{L},\, {I}_{R},\, {I}_{avg},\, {I}_{eq})\,</math> in node <math display="inline">m</math> of the intersection.

{| style="width: 100%;border-collapse: collapse;" 
|-
|  style="vertical-align: top;"|[[File:Draft_Souza_246201888-image9-c.png|centre|282x282px]]<span style="text-align: center; font-size: 75%;">  </span>
! style="text-align: right;vertical-align: top;" |[[File:Draft_Souza_246201888-image10-c.png|centre|282x282px]]<span style="text-align: center; font-size: 75%;">  </span>
|-
|  style="text-align: center;vertical-align: top;"|(a)
|  style="text-align: center;vertical-align: top;"|(b)
|-
|  style="vertical-align: top;"|[[File:Draft_Souza_246201888-image11-c.png|centre|282x282px]]<span style="text-align: center; font-size: 75%;">  </span>
! style="text-align: right;vertical-align: top;" |[[File:Draft_Souza_246201888-image12-c.png|centre|282x282px]]<span style="text-align: center; font-size: 75%;">  </span>
|-
|  style="text-align: center;vertical-align: top;"|(c)
|  style="text-align: center;vertical-align: top;"|(d)
|-
|  colspan='2'  style="vertical-align: top;"|<span id='_Ref72141256'></span>Figure 4 – Estimated error ( <math display="inline">U</math>) and true error ( <math display="inline">\epsilon \,</math> ) for stepped column: CASE 1. a) <math display="inline">{I}_{eq}</math>, b) <math display="inline">{I}_{R}</math>, c) <math display="inline">{I}_{avg}</math>, d) <math display="inline">{I}_{L}</math>
|}


{| style="width: 100%;border-collapse: collapse;" 
|-
|  style="vertical-align: top;"|[[File:Draft_Souza_246201888-image13-c.png|centre|282x282px]]<span style="text-align: center; font-size: 75%;">  </span>
|  style="text-align: right;vertical-align: top;"|[[File:Draft_Souza_246201888-image14-c.png|centre|282x282px]]<span style="text-align: center; font-size: 75%;">  </span>
|-
|  style="text-align: center;vertical-align: top;"|(a)
|  style="text-align: center;vertical-align: top;"|(b)
|-
|  style="vertical-align: top;"|[[File:Draft_Souza_246201888-image15-c.png|centre|282x282px]]<span style="text-align: center; font-size: 75%;">  </span>
|  style="vertical-align: top;"|[[File:Draft_Souza_246201888-image16-c.png|centre|282x282px]]<span style="text-align: center; font-size: 75%;">  </span>
|-
|  style="text-align: center;vertical-align: top;"|(c)
|  style="text-align: center;vertical-align: top;"|(d)
|-
|  colspan='2' |<span id='_Ref72141259'></span>Figure 5 - Estimated error ( <math display="inline">U</math>) and true error ( <math display="inline">\epsilon \,</math> ) for stepped column CASE 2. a) <math display="inline">{I}_{eq}</math>, b) <math display="inline">{I}_{R}</math>, c) <math display="inline">{I}_{avg}</math>, d) <math display="inline">{I}_{L}</math>
|}


According to the results presented, the convergence behavior between the stepped columns (CASE 1 and CASE 2) is quite similar. The results were consistent and convergent when the equivalent moment of inertia ( <math display="inline">{I}_{eq})</math> was used, reaching 2<sup>nd</sup> convergence order for the Richardson Error Estimator and 4<sup>th</sup> convergence order using Convergence Error Estimator. Numerical instability occurs approximately for similar meshes, in both cases, that indicated an instability for computational limit.

In addition, one could be seen that the results of the convergent solution ( <math display="inline">{\phi }_{c})</math> for CASE 1 obtained an error approximately to <math display="inline">{10}^{-10}</math> and for CASE 2 the error was about <math display="inline">{10}^{-8}</math>, which the best convergence value occurs for 1024 elements. It is valid to emphasize that the matrices generated for a much large number of elements are robust, and the convergence is constrained by the computer precision used in the simulation.

As expected, using moment of inertia ( <math display="inline">{I}_{R}</math>'', '' <math display="inline">{I}_{avg}</math> and <math display="inline">{I}_{L}</math>) have inconsistent results with errors between <math display="inline">{10}^{-3}</math> and <math display="inline">{10}^{-2}</math>. This happens because the continuity conditions on bending moment are not satisfied at the point <math display="inline">m</math>of intersection between different cross-sections [12], therefore, the apparent order of convergence does not reach the asymptotic order of error.

Lastly, the results confirmed what [4] reinforced, that is in an interval convergent <math display="inline">Pu</math> is advisable to use <math display="inline">{\phi }_{c}</math> instead <math display="inline">{\phi }_{1}</math> because the true discretization error <math display="inline">\epsilon ({\phi }_{c})</math> is less than <math display="inline">\epsilon ({\phi }_{1})</math>. Furthermore, this work presented that the use of extrapolation in FDM provides better results when compared with conventional FEM results.

==7. Conclusions==

This article verified the order of convergence of the critical buckling load of a simply supported stepped column. For this, Richardson's Extrapolation and the Convergence Error Estimator were used to estimating and reducing the numerical error of the buckling problem by the Finite Difference Method.

The study contemplated two cases with different ratios between the moment of inertia of the cross-sections and evaluated the consistency of the results at the intersection of the sections (the point at which the change in the cross-section of the column occurs). To verify the consistency of the solution, different moments of inertia were proposed at the point of intersection for both cases.

The methodology using the equivalent moment of inertia at the intersection node proved to be monotonic convergent with asymptotic order 2<sup>nd</sup> for Richardson's Error Estimator and asymptotic order 4<sup>th</sup> for Convergence Error Estimator, being considered consistent. These results showed that, for the same number of elements, the convergent solution has errors twice as small as that of the numerical solution. In this case, the error of the numerical solution was approximately 10<sup>-4</sup> and the error of the convergent solution was 10<sup>-8</sup> and 10<sup>-10</sup>, depending on the ratio between moments of inertia.

This work also showed that using inconsistent moments of inertia at the point of intersection precludes the apparent order of error from reaching the asymptotic order of the proposed problem. This confirmed that the study of convergence is extremely necessary for the analysis of possible inconsistencies or numerical modeling problems.

Finally, the use of error estimators proved to be adequate to verify the reliability of the discretization errors for the proposed problem, in addition to being possible to define the regions with asymptotic or unstable behavior of the numerical solutions. Moreover, applied extrapolations technics in FDM can be possible reaching minor errors than the model by FEM using commercial software.

==8. Acknowledgements==

This study was financed in part by the Coordination for the Improvement of Higher Education Personnel – Brazil (CAPES), Brazilian Council for Scientific and Technological Development (CNPq) and Research Support Foundation of Federal District (FAPDF). The authors are grateful to Group of Experimental and Computational Mechanics (UnB-FGA/GMEC) for provide computational resources and become possible the work development. The second author acknowledges Laboratory of Vibration and Dynamic Metrology and Group of Dynamics of Systems (UnB-FT/EnM/GDS) for its support to make possible the present experiments<span style="text-align: center; font-size: 75%;">.</span>

==9. References==

[1] L. F. Richardson, “IX. The approximate arithmetical solution by finite differences of physical problems involving differential equations, with an application to the stresses in a masonry dam,” ''Philos. Trans. R. Soc. London. Ser. A, Contain. Pap. a Math. or Phys. Character'', vol. 210, pp. 307–357, 1911, doi: 10.1098/rsta.1911.0009.

[2] C. H. Marchi, L. A. Novak, C. D. Santiago, and A. P. da S. Vargas, “Highly accurate numerical solutions with repeated Richardson extrapolation for 2D laplace equation,” ''Appl. Math. Model.'', vol. 37, no. 12–13, pp. 7386–7397, 2013, doi: 10.1016/j.apm.2013.02.043.

[3] C. H. Marchi ''et al.'', “Polynomial interpolation with repeated Richardson extrapolation to reduce discretization error in CFD,” ''Appl. Math. Model.'', vol. 40, no. 21–22, pp. 8872–8885, 2016, doi: 10.1016/j.apm.2016.05.029.

[4] C. H. Marchi and A. F. C. Silva, “Multi-dimensional discretization error estimation for convergent apparent order,” ''J. Brazilian Soc. Mech. Sci. Eng.'', vol. 27, no. 4, pp. 432–439, 2005, doi: 10.1590/S1678-58782005000400012.

[5] C. H. Marchi and A. F. C. Silva, “Unidimensional Numerical Solution Error Estimation For Convergent Apparent Order,” ''Numer. Heat Transf. Part B Fundam.'', vol. 42, no. 2, pp. 167–188, 2002, doi: 10.1080/10407790190053888.

[6] Y. M. Wang, “Error analysis of a compact finite difference method for fourth-order nonlinear elliptic boundary value problems,” ''Appl. Numer. Math.'', vol. 120, no. 13, pp. 53–67, 2017, doi: 10.1016/j.apnum.2017.04.011.

[7] N. Challamel, V. Picandet, I. Elishakoff, C. M. Wang, B. Collet, and T. Michelitsch, “On Nonlocal Computation of Eigenfrequencies of Beams Using Finite Difference and Finite Element Methods,” ''Int. J. Struct. Stab. Dyn.'', vol. 15, no. 7, pp. 1–18, 2015, doi: 10.1142/S0219455415400088.

[8] N. Challamel, V. Picandet, B. Collet, T. Michelitsch, I. Elishakoff, and C. M. Wang, “Revisiting finite difference and finite element methods applied to structural mechanics within enriched continua,” ''Eur. J. Mech. A/Solids'', vol. 53, pp. 107–120, 2015, doi: 10.1016/j.euromechsol.2015.03.003.

[9] S. P. Timoshenko and J. M. Gere, “Theory of Elastic Stability-Dover Publications.” pp. 113–117, 1963.

[10] M. Soltani and A. Sistani, “Elastic stability of columns with variable flexural rigidity under arbitrary axial load using the finite difference method,” ''Journals Numer. Method Civ. Eng.'', vol. 1, no. 4, pp. 23–31, 2017, [Online]. Available: [http://nmce.kntu.ac.ir/article-1-91-en.pdf%0Ahttp://nmce.kntu.ac.ir/browse.php?a_code=A-10-50-3&sid=1&slc_lang=en. http://nmce.kntu.ac.ir/article-1-91-en.pdf%0Ahttp://nmce.kntu.ac.ir/browse.php?a_code=A-10-50-3&sid=1&slc_lang=en.]

[11] Michael O’Rourke, “Buckling load for nonuniform columns,” vol. 7, no. 4, pp. 717–720, 1977.

[12] A.Krishnan, G. George, and P. Malathi, “Use of finite difference method in the study of stepped beams,” ''Int. J. Mech. Eng. Educ.'', vol. 26, no. 1, 1996.

[13] X. Lu, Q. Wu, and Y. Khan, “The Analysis of Dynamic Buckling of an Impacted Column using Difference Methods,” ''Int. J. Comput. Methods'', vol. 15, no. 4, 2018, doi: 10.1142/S0219876218500251.

[14] Z. R. Lu, M. Huang, J. K. Liu, W. H. Chen, and W. Y. Liao, “Vibration analysis of multiple-stepped beams with the composite element model,” ''J. Sound Vib.'', vol. 322, no. 4–5, pp. 1070–1080, May 2009, doi: 10.1016/J.JSV.2008.11.041.

[15] T. D. Dang, R. K. Kapania, and M. J. Patil, “Ritz analysis of discontinuous beams using local trigonometric functions,” ''Comput. Mech.'', vol. 47, no. 3, pp. 235–250, 2011, doi: 10.1007/s00466-010-0540-y.

[16] S. M. Baǧdatli, E. Özkaya, H. A. Özyiǧit, and A. Tekin, “Non linear vibrations of stepped beam systems using artificial neural networks,” ''Struct. Eng. Mech.'', vol. 33, no. 1, pp. 15–30, 2009, doi: 10.12989/sem.2009.33.1.015.

[17] F. Yao, W. Meng, J. Zhao, Z. She, G. Shi, and H. Liu, “Buckling theoretical analysis on all-terrain crane telescopic boom with n-stepped sections,” ''J. Mech. Sci. Technol.'', vol. 32, no. 8, pp. 3637–3644, 2018, doi: 10.1007/s12206-018-0715-6.

[18] M. I. Salama, “Elastic Buckling of Single-Stepped Columns with End and Intermediate Axial Loads,” ''Comput. Eng. Phys. Model.'', vol. 1, no. 4, pp. 25–35, 2019, doi: 10.22115/CEPM.2018.147218.1046.

[19] J. Marcinowski and V. Sakharov, “Buckling Resistance of Two-Segment Stepped Steel Columns,” 2021.

[20] W. Tian, J. Hao, and W. Zhong, “Buckling of stepped columns considering the interaction effect among columns,” ''J. Constr. Steel Res.'', vol. 177, p. 106416, 2021, doi: 10.1016/j.jcsr.2020.106416.

[21] M.J. Iremonger, “Finite difference buckling analysis of non-uniform columns,” ''Comput. Struct.'', no. 3, 1980.