Pith. sign in

REVIEW 3 major objections 5 minor 18 references

An Efficient Numerical Method for an Approximate Solution of the Beam Equation

T0 review · 3 major / 5 minor · reviewed 2026-08-07 · deepseek-v4-flash

Pith's one-line read A horizontal method of lines combining Crank-Nicolson time stepping with Hermite cubic finite elements approximates the unsteady damped beam equation and its derivatives with errors of order $O(h^4 + \tau^2)$.

desk verdict A promising HMOL/FEM scheme for the beam equation, but the printed Crank-Nicolson elimination in Eqs. (11)-(13) is algebraically wrong and must be fixed for the method to be reproducible. read the letter →

arxiv 2506.03794 v1 pith:GUT3WS2M submitted 2025-06-04 math.NA cs.NA

classification math.NAcs.NA MSC 65M6065M2074K10
keywords Euler-BernoullibeamequationhorizontalmethodoflinesCrank-NicolsonschemeHermitecubicfiniteelementsorderaccuracymeshparameteroptimizationapproximatederivatives
verification ladder T0 review T1 audit T2 compute T3 formal

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The reading

The paper proposes a horizontal method of lines (HMOL) scheme for the damped, unsteady Euler-Bernoulli beam equation and claims it produces accurate approximations of the deflection $w$ and its derivatives $w_x$, $w_t$, and $w_{xt}$. The scheme rewrites the beam equation as a first-order system in time using $z=w_t$, applies a Crank-Nicolson (trapezium) one-step update to the temporal derivative, and solves the resulting family of steady beam problems with Hermite cubic finite elements. On a test problem with a known exact solution the paper reports errors of order $O(h^4 + \tau^2)$ and recommends the mesh relation $\tau \approx 2 h^2$. This matters because velocity and slope data enter beam identification problems and are usually hard to compute stably; the method produces them directly as part of the computation.

What carries the argument

The central mechanism is the splitting of the beam equation into a first-order system in time followed by a spatial finite element solve at each time level. With $z = w_t$, the Crank-Nicolson update for the pair $(w,z)$ is applied and the velocity is then eliminated to produce a single steady beam equation $(r W_{n+1}'')'' + \tilde{s}(x) W_{n+1} - (j W_{n+1}')' = \tilde{g}(x)$, where $\tilde{s}(x) = 4/\tau^2 + 2\eta(x)/\tau + s(x)$. Hermite cubic basis functions carry both nodal values and first derivatives, so the slopes $w_x$ and $w_{xt}$ come from the same degrees of freedom. Integration by parts converts the fourth-order term $(r W_n'')''$ that appears in the load into a computable bilinear form, and element matrices are assembled with numerical quadrature. This separation lets the temporal integrator be swapped without rebuilding the spatial solver, which is the main design advantage the paper claims.

What would settle it

Re-derive the elimination of $Z_{n+1}$ from the Crank-Nicolson pair and compare the coefficients printed in (13), or implement the published equations exactly and check whether the reported errors of order $O(h^4 + \tau^2)$ still appear; a mismatch would mean the paper has not documented the scheme it tested.

Watch

Extended reading notes

Core claim

The central claim is that the HMOL iteration constructed from the Crank-Nicolson update and Hermite cubic spatial discretization solves the damped dynamic beam equation with errors $\|W_h - w\| = O(h^4 + \tau^2)$ and $\|Z_h - w_t\| = O(h^4 + \tau^2)$, with the velocity error deteriorating slightly on the finest meshes. The algorithm first eliminates the velocity variable $z=w_t$ algebraically from the Crank-Nicolson pair so that each time step reduces to a single steady Euler-Bernoulli boundary value problem for $W_{n+1}$. Each steady problem is solved with the same finite element machinery used for the stationary beam, and logarithmic error plots show fourth-order convergence in space and second-order convergence in time, with $\tau \approx 2 h^2$ the recommended operating point. The paper further claims that the Hermite basis supplies stable approximations of $w_x$, $w_t$, and $w_{xt}$ without any additional differentiation step.

Load-bearing premise

The whole method depends on the algebraic elimination step that turns the Crank-Nicolson update into a single steady beam equation; the printed formulas for this step are internally inconsistent, so the tests may not describe the algorithm that was actually implemented.

Editorial extensions

If this is right

  • Each time step reuses a standard steady beam finite element routine; only the right-hand side and the coefficient of the zeroth-order term change.
  • The quantities $w_x$, $w_t$, and $w_{xt}$ are by-products of the computation, so derivative data do not require separate difference schemes.
  • The observed convergence rates $O(h^4)$ in space and $O(\tau^2)$ in time hold simultaneously when $\tau \approx 2 h^2$, giving a practical rule for choosing meshes.
  • Because the temporal discretization is isolated at the first stage, a more accurate or more stable one-step scheme can be substituted without altering the finite element core.

Reading between the lines

Editorial extensions of the paper, not claims the author makes directly.

  • If the accuracy persists for non-constant damping and for loads that vary rapidly in time, the method could provide the boundary data needed by inverse problems that reconstruct loads or shear forces from measured deflection and slope.
  • The table shows the velocity error stops improving on the finest meshes, which suggests the time step, not the spatial mesh, controls $w_t$; a focused study varying $\tau$ alone for fixed small $h$ would test this directly.
  • Because $w_{xt}$ is produced directly, the scheme is a natural candidate for estimating pointwise bending moments or shear forces without post-processing, an extension the paper does not explore.
Share X Bluesky LinkedIn Reddit HN

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, and a circularity audit.

Referee Report

3 major / 5 minor

Summary. The paper proposes a horizontal method of lines (HMOL) for the damped unsteady Euler-Bernoulli beam equation (1). The PDE is rewritten as a first-order system (10), discretized in time with the Crank-Nicolson method, and each resulting steady beam equation is solved with Hermite cubic finite elements. The scheme yields simultaneous approximations of w, w_x, w_t, and w_xt. Numerical tests on a manufactured solution report O(h^4 + τ^2) accuracy with optimal mesh relation τ ≈ 2h^2. A stationary version is also analyzed and tested.

Significance. If the claims hold, the method would be an attractively simple alternative to vertical method of lines for beam problems, and the direct availability of derivative approximations could be valuable for inverse problems. The paper does not provide a convergence proof or comparisons with existing methods; its evidence is purely experimental. The MATLAB snippet for element assembly (Section 2.1) is a helpful reproducibility aid, but the full algorithm is not sufficiently specified because of an algebraic inconsistency in Eqs. (11)-(13).

major comments (3)
  1. [Section 3, Eqs. (11)-(13)] The Crank-Nicolson reduction is algebraically inconsistent. The correct first-order update is W_{n+1}=W_n+(τ/2)(Z_n+Z_{n+1}), which implies Z_{n+1}=2/τ(W_{n+1}-W_n)-Z_n; the first line of (11) drops the -Z_n term. Furthermore, eliminating Z_n+Z_{n+1} from the second component via the correct relation gives a W_n coefficient of 4/τ^2+2η/τ-s in (13), not 4/τ^2+2η/τ+s as printed. Since no code is provided, a reader cannot determine which scheme produced the reported errors, and the convergence rates quoted in Section 3.1 may not pertain to the published algorithm. This is a load-bearing reproducibility failure in the method-defining derivation.
  2. [Section 3.1, Table 1] The errors for the velocity approximation V_h-w_t (last row of Table 1) are not monotone along the claimed optimal mesh relation τ≈2h^2: the error at (h,τ)=(1/80,1/3200) is 3.07×10^-7, larger than the 1.78×10^-7 at (1/40,1/800). This contradicts the stated O(h^4+τ^2) behavior and casts doubt on the advertised stable computation of w_t. The sentence 'the error rate for Vh − wt is deteriorated with the use of finer meshes' does not explain this increase; the paper should either provide a quantitative explanation or restrict the derivative-accuracy claim.
  3. [Section 3.1, paragraph on order of accuracy] The authors state that a theoretical analysis is challenging and resort to manufactured-solution tests. In a numerical-analysis journal this is acceptable only if the experiments are unambiguous; here they are not, because of the algebraic inconsistency in Eqs. (11)-(13) and the non-monotone derivative errors in Table 1. A convergence proof, or at least a rigorous error estimate for the corrected scheme, is needed to support the conclusion that the method is 'highly accurate and reliable'.
minor comments (5)
  1. [Section 2, after Eq. (3)] The notation ℓ is used in the assumption s,j ∈ L∞(0,ℓ) while l is used elsewhere; please unify the symbol for the beam length.
  2. [Section 2.1] The partition is described by m = 1 : M+1 for |e_m|; since there are M elements, the index should be m = 1 : M (the subsequent basis definitions are consistent with M elements).
  3. [Section 3.1] The norms used in Figures 7-8 and Table 1 are not fully specified; the text defines the C(0,T;L^2) norm only after the table, and the captions only state 'Error'.
  4. [References [16],[17]] The paper refers to existence and uniqueness results for the clamped case and says the proof can be easily adapted to the hinged-hinged case, but it does not state which boundary-term modifications are involved; a brief remark would make the discussion self-contained.
  5. [Section 2.1 and Section 3] The MATLAB snippet for the element stiffness matrix is helpful, but a full pseudocode of the HMOL iteration (including the projection of initial data and the assembly with the corrected \tilde{s} coefficient) would substantially improve reproducibility.

Circularity Check

0 steps flagged · score 0.0 of 10

No circularity found: the reported convergence rates are measured against an independent manufactured exact solution, not fitted or imported from self-citations; the algebraic inconsistency in Eqs. (11)-(13) is a correctness issue, not circularity.

full rationale

The central claim is an empirical convergence study of a Crank-Nicolson/Hermite-FEM scheme. The error tables and log-log plots compare Wh, Zh, ∂xWh, ∂xZh against a manufactured exact solution w(x,t)=(1-x+sin^2 x) sin t; no constant is fitted to make errors match O(h^4+tau^2), and no parameter is calibrated to the benchmark. The temporal scheme is explicitly stated as Crank-Nicolson, and the FEM solve is standard Hermite cubics; neither step is defined in terms of the quantity being 'predicted.' The self-citations to prior beam/inverse-problem papers (refs. [7]-[10], [16]) are background existence/uniqueness and technique credits, not load-bearing for the reported convergence rates; the well-posedness remark explicitly says the proof 'can be easily adapted' from standard Galerkin arguments. The only notable defect is an internal algebraic inconsistency in the printed Crank-Nicolson elimination: Eq. (11) drops the -Z_n term and Eq. (13) carries a sign error in the s W_n coefficient. That is a correctness/reproducibility issue, not a circularity, because the reported errors are still measured against an independent exact solution rather than being forced by construction. Hence no circular step is found.

Assumptions & free parameters 1 free parameters · 3 assumptions · 0 invented entities

The scheme's convergence claims rest on standard FEM assumptions and on the correctness of the Crank-Nicolson algebra. The empirical relation tau about 2 h^2 is the only fitted choice. No new physical entities are introduced.

free parameters (1)
  • mesh parameter relation tau about 2 h^2 = tau = 2 h^2 (empirical)
    In Section 3.1, after systematic sweeps of h and tau, the authors state that the optimal mesh relation is tau about 2 h^2. This is chosen from numerical experiments, not derived, and is used to balance the observed O(h^4) and O(tau^2) errors.
assumptions (3)
  • standard math The weak problem (5) has a unique solution via the Lax-Milgram theorem under assumptions (4).
    Section 2 invokes Lax-Milgram and defines a V-norm but does not prove coercivity for the general L-infinity coefficients s and j; sign-indefinite coefficients are not excluded.
  • domain assumption Existence and uniqueness of the weak solution of (1) can be established by the Galerkin method, and the proof for the clamped case adapts to the hinged-hinged case.
    Section 3 states this without proof, relying on references [16], [17], and [18]; the adaptation to the hinged-hinged dynamic problem is not shown.
  • ad hoc to paper The Crank-Nicolson discretization algebra leading to equations (11)-(13) is correct.
    The steady equation and source term are then used as the algorithm; the displayed identities do not follow as written from the Crank-Nicolson equations above, so this is a load-bearing premise of the method.

how reviews work

0 comments
Cite this review

Pith. "Pith review of An Efficient Numerical Method for an Approximate Solution of the Beam Equation." pith.science (2026). https://pith.science/paper/GUT3WS2M

@misc{pith2026250603794,
  author       = {Pith},
  title        = {Pith review of: An Efficient Numerical Method for an Approximate Solution of the Beam Equation},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/GUT3WS2M}},
  note         = {Machine review of arXiv:2506.03794}
}
read the original abstract

In this paper, we propose a horizontal type method of lines numerical scheme for the unsteady Euler-Bernoulli beam equation. The problem is initially reformulated as a first order system of initial value problems and a suitable one-step difference scheme is used for the highest order temporal derivative which leads to a system of steady beam equations. Then resulted family of steady problems is solved iteratively by the finite element method with Hermite cubic basis functions. This iterative procedure leads to approximations for both the solution of the unsteady problem and its derivatives. All these approximations are compared with the exact ones to illustrate the performance of the proposed method. Moreover, the optimization of the mesh parameters is discussed for both steady and unsteady problems by logarithmic scale plot.

Figures

Figures reproduced from arXiv: 2506.03794 by the authors.

Figure 1
Figure 1. Geometry of the problem: Hinged-hinged beam model. The non-homogenous dynamic Euler-Bernoulli beam with supported (hinged) boundary conditions is illustrated in [PITH_FULL_IMAGE:figures/full_fig_p002_1.png] view at source ↗
Figure 2
Figure 2. Comparison of the exact solution and approximate solution with derivatives (left), error plot in log scale (right). [PITH_FULL_IMAGE:figures/full_fig_p005_2.png] view at source ↗
Figure 3
Figure 3. Approximate solution Wh(x, t) (left) and corresponding error (right) [PITH_FULL_IMAGE:figures/full_fig_p007_3.png] view at source ↗
Figures from the paper (5 more)
Figure 4
Figure 4. Figure 4: Approximate solution ∂xWh(x, t) (left) and corresponding error (right) [PITH_FULL_IMAGE:figures/full_fig_p007_4.png]
Figure 5
Figure 5. Figure 5: Approximate solution Zh(x, t) (left) and corresponding error (right). increases when capturing higher order derivatives, often making it nearly impossible to obtain stable results. However, the proposed algorithm approximates accurately not only wt(x, t) but also wxt(x…
Figure 6
Figure 6. Figure 6: Approximate solution ∂xZh(x, t) (left) and corresponding error (right). shown in Figures (5) and (6). This potential can be regarded as one of the most significant advantages of the algorithm. For hybrids methods, generally, the order of accuracy cannot be determined d…
Figure 7
Figure 7. Figure 7: Error in log scale for Wh (left) and Zh (right) throughout temporal step size τ . Following table confirms that the errors Wh −w and Vh −wt in C(0, T; L 2 (0, l)) norm are O(h 4 +τ 2 ). As noted in error plots, the error rate for Vh − wt is deteriorated with the use of…
Figure 8
Figure 8. Figure 8: Error in log scale for Wh (left) and Zh (right) throughout spatial step size h. (h, τ ) (1/10, 1/50) (1/20, 1/200) (1/40, 1/800) (1/80, 1/3200) ||Wh − w|| 1.78 10−6 1.12 10−7 6.99 10−9 4.34 10−10 ||Vh − wt || 4.48 10−5 2.79 10−6 1.78 10−7 3.07 10−7 4. Conclusion In con…

Discussion (0). Sign in to comment.

Reference graph

Works this paper leans on

18 extracted references · 18 canonical work pages

  1. [1]

    Vinod K G, Gopalakrishnan S and Ganguli R, Free vibration and wave propagation analysis of uniform and tapered rotating beams using spectrally formulated finite elements.International Journal of Solids and Structures, 44 5875-5893, 2007. 9

  2. [2]

    Thesis (UK: University of Bristol), 2000

    Antognozzi M, Investigation of the shear force contrast mechanism in transverse dynamic force microscopy Ph.D. Thesis (UK: University of Bristol), 2000

  3. [3]

    Nguyen T and et. al. Estimation of the shear force in transverse dynamic force microscopy using a sliding mode observer. AIP Advances, 5 097157, 2015

  4. [4]

    Inverse Problems, 29 015006, 2013

    Bao G and Xu X, An inverse random source problem in quantifying the elastic modulus of nanoma- terials. Inverse Problems, 29 015006, 2013

  5. [5]

    and Sankar A, A Finite Element Solution of the Beam Equation via MATLAB

    Gunakala S R, Comissiong D M G, Jordan K. and Sankar A, A Finite Element Solution of the Beam Equation via MATLAB. International Journal of Applied Science and Technology, Vol. 2 No. 8, 2012

  6. [6]

    CRC Press, 2000

    Kwon Y W and Bang H, The Finite Element Method Using Matlab. CRC Press, 2000

  7. [7]

    Hasanov A, Kawano A and Baysal O,Exponential stability of damped Euler-Bernoulli beam controlled by boundary springs and dampers , Journal of Mathematical Analysis and Applications, Volume 533, Issue 2, 128031, 2024

  8. [8]

    Automatica 71 (2016) 106- 117, 2016

    Hasanov A and Baysal O, Identification of unknown temporal and spatial load distributions in a vibrating Euler-Bernoulli beam from Dirichlet boundary measured data. Automatica 71 (2016) 106- 117, 2016

Show all 18 references
  1. [9]

    Inverse Ill-posed Probl

    Hasanov A, Baysal O and Itou H, Identification of an unknown shear force in a cantilever Euler- Bernoulli beam from measured boundary bending moment J. Inverse Ill-posed Probl. 27(6), 859-876, 2019

  2. [10]

    Mathematics in Engineering, 6(1), 137–154, 2024

    Hasanov A, Kawano A and Baysal O, Reconstruction of shear force in Atomic Force Microscopy from measured displacement of the cone-shaped cantilever tip. Mathematics in Engineering, 6(1), 137–154, 2024

  3. [11]

    CRC Press, 2017

    Surana K S and Reddy J N The Finite Element Method for Initial Value Problems Mathematics and Computations. CRC Press, 2017

  4. [12]

    Springer, 2011

    Brezis H, Functional Analysis, Sobolev Spaces and Partial Differential Equations. Springer, 2011

  5. [13]

    Cambridge University Press, 2012

    S¨ uli E and Mayers D F,An Introduction to Numerical Analysis. Cambridge University Press, 2012

  6. [14]

    Mayer A, A Simplified Calculation of Reduced HCT–Basis Functions in a Finite Element Context , Computational Methods in Applied Mathematics 12-4, 486-499, 2012

  7. [15]

    Wiley, 2020

    Asadzadeh M, An Introduction to the Finite Element Method (FEM) for Differential Equations. Wiley, 2020

  8. [16]

    Applied Math- ematics Letters, vol

    Baysal O and Hasanov A, Solvability of the clamped Euler-Bernoulli beam equation. Applied Math- ematics Letters, vol. 93, 85-90, 2019

  9. [17]

    2nd ed, Springer, New York, 2021

    Hasanov A and Romanov A G, Introduction to Inverse Problems for Differential Equations. 2nd ed, Springer, New York, 2021

  10. [18]

    Rhode Island: American Mathematical Society, 2002

    Evans L C, Partial Differential Equations. Rhode Island: American Mathematical Society, 2002. 10

Pith tools

Reviewed August 7, 2026 · model on record in the stance chip above.