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 →
The pith
A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.
The reading
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.
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
- 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.
Editorial analysis
A structured set of objections, weighed in public.
Referee Report
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)
- [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.
- [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.
- [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)
- [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.
- [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).
- [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'.
- [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.
- [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
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
free parameters (1)
- mesh parameter relation tau about 2 h^2 =
tau = 2 h^2 (empirical)
assumptions (3)
- standard math The weak problem (5) has a unique solution via the Lax-Milgram theorem under assumptions (4).
- 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.
- ad hoc to paper The Crank-Nicolson discretization algebra leading to equations (11)-(13) is correct.
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 from the paper (5 more)
Reference graph
Works this paper leans on
-
[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
work page 2007
-
[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
work page 2000
-
[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
work page 2015
-
[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
work page 2013
-
[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
work page 2012
-
[6]
Kwon Y W and Bang H, The Finite Element Method Using Matlab. CRC Press, 2000
work page 2000
-
[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
work page 2024
-
[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
work page 2016
Show all 18 references
-
[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
2019
-
[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
2024
-
[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
2017
-
[12]
Springer, 2011
Brezis H, Functional Analysis, Sobolev Spaces and Partial Differential Equations. Springer, 2011
2011
-
[13]
Cambridge University Press, 2012
S¨ uli E and Mayers D F,An Introduction to Numerical Analysis. Cambridge University Press, 2012
2012
-
[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
2012
-
[15]
Wiley, 2020
Asadzadeh M, An Introduction to the Finite Element Method (FEM) for Differential Equations. Wiley, 2020
2020
-
[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
2019
-
[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
2021
-
[18]
Rhode Island: American Mathematical Society, 2002
Evans L C, Partial Differential Equations. Rhode Island: American Mathematical Society, 2002. 10
2002
Reviewed August 7, 2026 · model on record in the stance chip above.
Discussion (0). Sign in to comment.