Next Article in Journal
Newton-like Normal S-iteration under Weak Conditions
Next Article in Special Issue
Relational Contractions Involving Shifting Distance Functions with Applications to Boundary Value Problems
Previous Article in Journal
A Space-Time Legendre-Petrov-Galerkin Method for Third-Order Differential Equations
Previous Article in Special Issue
Semi-Hyers–Ulam–Rassias Stability of Some Volterra Integro-Differential Equations via Laplace Transform
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

On the Simulations of Second-Order Oscillatory Problems with Applications to Physical Systems

1
Department of Mathematics, School of Sciences, Federal College of Education, Pankshin 933105, Nigeria
2
Department of Mathematics, Faculty of Natural Sciences, University of Jos, Jos 930003, Nigeria
3
Department of Mathematics, Faculty of Sciences, University of Maragheh, Maragheh 83111-55181, Iran
4
College of Mathematics and Computer Science, Zhejiang Normal University, Jinhua 321004, China
*
Author to whom correspondence should be addressed.
Axioms 2023, 12(3), 282; https://doi.org/10.3390/axioms12030282
Submission received: 21 January 2023 / Revised: 24 February 2023 / Accepted: 28 February 2023 / Published: 8 March 2023
(This article belongs to the Special Issue Differential Equations and Related Topics)

Abstract

:
Second-order oscillatory problems have been found to be applicable in studying various phenomena in science and engineering; this is because these problems have the capabilities of replicating different aspects of the real world. In this research, a new hybrid method shall be formulated for the simulations of second-order oscillatory problems with applications to physical systems. The proposed method shall be formulated using the procedure of interpolation and collocation by adopting power series as basis function. In formulating the method, off-step points were introduced within the interval of integration in order to bypass the Dahlquist barrier, improve the accuracy of the method and also upgrade the order of consistence of the method. The paper further validated the some properties of the hybrid method derived and from the results obtained; the new method was found to be consistent, convergent and stable. The simulation results generated as a result of the application of the new method on some second-order oscillatory differential equations also showed that the new hybrid method is computationally reliable.
MSC:
65D30; 65L05; 65L07; 65L11; 65L20

1. Introduction

Over the years, second-order differential equations with oscillatory or periodic solutions have attracted the attention of researchers. In this research, a new hybrid method is formulated for the simulations of second-order oscillatory problems (with applications to physical systems) of the form,
y ( t ) = ψ t , y ( t ) , y ( t ) y ( t 0 ) = y 0 ,   y ( t 0 ) = y 0
where  ψ : × 2 n  is a sufficiently differentiable function satisfying Theorem 1. Equations of the form (1) have been known to exhibit special properties like damping, stiffness and chaos, which make them tedious to solve [1]. In fact, some of such problems do not have closed form solution(s), thus, numerical methods have to be developed in order to approximate their solutions or at least simulate their behaviors and nature.
Theorem 1.
Ref. [2]
“Let,
y ( n ) ( t ) = ψ t ,   y ( t ) ,   y ( t ) , ,   y ( n 1 ) ( t ) ,     y ( k ) ( t 0 ) = y k
where  k = 0 , 1 , 2 , , ( n 1 ) ,  y ( 0 ) = y ,  y ( t )  and  ψ  are functions defined on  n n  being a region defined by the inequalities  t 0 t t 0 + a ,     s k y k b k ,     k = 0 , 1 , , n 1 , where  y k 0  for  k > 0 . Let the function  ψ ( t , s 0 , s 1 , , s n 1 )  in (1) be continuous, non-decreasing and non-negative in  t , and non-decreasing and continuous in  s k  for  k = 0 , 1 , , n 1  in  n . Furthermore, if  ψ ( t , y 0 , y 1 , , y n 1 ) 0  in  n  for  t > t 0 , then problem (1) has at most one solution in  n ”.
According to [3], the influence upon or within an oscillating system that reduces, restricts or prevents its oscillation is called damping. On the other hand, simulation is the recreation of a scenario, often repeatedly, with the aid of a mathematical model, so that the likelihood of various outcomes can be more accurately estimated.
Second-order oscillatory problems have been found to be applicable in modeling physical systems in different areas of human endeavor. These areas include weak signal detection [4], projectile motion with aerodynamic drag and half cylinder rolling on a horizontal plane [5], vibrations [6,7], theory of stiffness [8], contractions [9], damping [10,11,12], phase-lags [13,14] and theory of oscillations [15,16]. Other areas of applications include special domains [17,18,19], quadrotor flight control [20], rings [21], dynamical systems [22] and epidemic models [23].
Some scholars proposed methods for the solutions and simulations of second order oscillatory problems. These methods among others are Laplace decomposition algorithm [24], Adomian decomposition method [25], Taylor matrix method [26], differential transform method [27,28], variational iteration method [29,30,31], Runge-Kutta method [32] and hybrid block methods [33,34,35,36,37,38,39,40]. The authors in [41] for instance, developed hybrid method for solving second order problems while [42] proposed an accuracy-preserving algorithm of order six, for integrating second order oscillatory problems. Suffice to mention that some of these methods have some disadvantages like low rate of convergence, large number of function evaluation per step, low order of accuracy, non-flexibility in terms of step size change, among others.
It is in view of the setbacks highlighted above, that a continuous hybrid method with fewer function evaluations, higher order of accuracy and high rate of convergence is proposed in this research. The proposed method is also easy to program and is efficient with respect to economy of computer time.
Definition 1.
(Hybrid Method)
“A linear multistep method is called hybrid if off-grid point(s) is/are incorporated between conventional points within integration interval.”

2. Main Method

2.1. Derivation of the Method

The derivation of a new method for simulations of oscillatory second-order problems is carried out in this section. The method which shall be derived is of the form,
A ( 0 ) Y m ( i ) = i = 0 1 i h i e i y n ( i ) + h μ i d i f ( y n ) + b i f ( Y m )
where
Y m ( i ) = y n + i ( i )   y n + j ( i )     y n + 1 ( i ) T ,   f ( Y m ) = f n + i   f n + j     f n + 1 T y n ( i ) = y n i ( i )   y n j ( i )     y n ( i ) T ,   f ( y n ) = f n i   f n j     f n T
A ( 0 ) = ( ω 1 ) × ( ω 1 )  is an identity matrix,  μ  is the differential equation’s order, and  ( i )  is the power of derivative of the method. Note that  h  is the step-size calculated as  h = t n + 1 t n n = 0 , 1 , 2 , , N . The oscillatory problem (1) is solved in the  N  non-overlapping blocks.
At  i = 0 , the  ( ω 1 ) × ( ω 1 )  matrices  e 0 , e 1 , d 0  and  b 0  are evaluated and at  i = 1  (the first derivative), the  ( ω 1 ) × ( ω 1 )  matrices  e 1 , d 1  and  b 1  are evaluated. Let the computed solution of the oscillatory problem (1) be approximated by the power series
y ( t ) = j = 0 ω + τ 1 υ j t j
The condition  ω + τ  (where  ω  is the number of collocation points and  τ  is the number of interpolation points) is then imposed on Equation (4), which gives the polynomials of degree  m = ω + τ 1  (at most) as follows
y ( t n + τ 1 ) = y n + τ 1 = j = 0 ω + τ 1 υ j t n + τ 1 j
y ( t n + τ 2 ) = y n + τ 2 = j = 0 ω + τ 1 υ j t n + τ 2 j
y ( t n + κ ) = f n + κ = j = 0 ω + τ 1 j ( j 1 ) υ j t n + κ j 2
where  κ = 0 , p , q , r , s , u , w , 1 = 0 , 1 / 7 ,   2 / 7 ,   3 / 7 ,   4 / 7 ,   5 / 7 ,   6 / 7 ,   1  and  τ 1 , τ 2 = 5 7 , 6 7 . This leads to a system of  ω + τ  equations each of degree at most  m  which is written compactly in matrix form as
T V = U
where
V = υ 0   υ 1   υ 2     υ 9   U = y n + τ 1   y n + τ 2   f n   f n + p   f n + q   f n + r   f n + s   f n + u   f n + w   f n + 1 T
and
T = 1 t n + τ 1 t n + τ 1 2 t n + τ 1 3 t n + τ 1 4 t n + τ 1 5 t n + τ 1 6 t n + τ 1 7 t n + τ 1 8 t n + τ 1 9 1 t n + τ 2 t n + τ 2 2 t n + τ 2 3 t n + τ 2 4 t n + τ 2 5 t n + τ 2 6 t n + τ 2 7 t n + τ 2 8 t n + τ 2 9 0 0 2 6 t n 12 t n 2 20 t n 3 30 t n 4 42 t n 5 56 t n 6 72 t n 7   0 0 2 6 t n + p 12 t n + p 2 20 t n + p 3 30 t n + p 4 42 t n + p 5 56 t n + p 6 72 t n + p 7   0 0 2 6 t n + q 12 t n + q 2 20 t n + q 3 30 t n + q 4 42 t n + q 5 56 t n + q 6 72 t n + q 7   0 0 2 6 t n + r 12 t n + r 2 20 t n + r 3 30 t n + r 4 42 t n + r 5 56 t n + r 6 72 t n + r 7 0 0 2 6 t n + s 12 t n + s 2 20 t n + s 3 30 t n + s 4 42 t n + s 5 56 t n + s 6 72 t n + s 7 0 0 2 6 t n + u 12 t n + u 2 20 t n + u 3 30 t n + u 4 42 t n + u 5 56 t n + u 6 72 t n + u 7 0 0 2 6 t n + w 12 t n + w 2 20 t n + w 3 30 t n + w 4 42 t n + w 5 56 t n + w 6 72 t n + w 7 0 0 2 6 t n + 1 12 t n + 1 2 20 t n + 1 3 3   0 t n + 1 4 42 t n + 1 5 56 t n + 1 6 72 t n + 1 7
Equation (8) is then solved using the Gaussian elimination method for all the  υ j s  which are parameters to be determined. These parameters are inputted into Equation (4) to give the continuous hybrid method
y ( t ) = τ = τ 1 , τ 2 α τ ( t ) y n + τ + h 2 j = 0 1 β j ( t ) f n + j + β κ ( t ) f n + κ ,   κ = p , q , r , s , u , w
where
α τ 1 ( x ) = 6 7 x α τ 2 ( x ) = 7 x 5 β 0 ( x ) = 1 88905600 201768035 x 9 1037664180 x 8 + 2272978680 x 7 2767104480 x 6 + 2047798494 x 5 945897960 x 4 + 268939440 x 3 44452800 x 2 + 3739561 x 110490 β p ( x ) = 1 88905600 1412376245 x 9 7004233215 x 8 + 14576711100 x 7 16454389140 x 6 + 10808648928 x 5 4047797880 x 4 + 726062400 x 3 19424888 x + 1999920 β q ( x ) = 1 29635200 1412376245 x 9 6744817170 x 8 + 13341396600 x 7 14033172720 x 6 + 8320372578 x 5 2659203540 x 4 + 363031200 x 3 + 1034257 x 1073130 β r ( x ) = 1 17781120 1412376245 x 9 6485401125 x 8 + 12204907260 x 7 12046787004 x 6 + 6590226384 x 5 1913981160 x 4 + 242020800 x 3 4634042 x + 1203780 β s ( x ) = 1 17781120 1412376245 x 9 6225985080 x 8 + 11167243080 x 7 10435936896 x 6 + 5389500690 x 5 1488427920 x 4 + 181515600 x 3 + 979363 x 1352670 β u ( x ) = 1 29635200 1412376245 x 9 5966569035 x 8 + 10228404060 x 7 9141327300 x 6 + 4540310208 x 5 1216154520 x 4 + 145212480 x 3 5101628 x + 2786520 β w ( x ) = 1 88905600 1412376245 x 9 5707152990 x 8 + 9388390200 x 7 8103663120 x 6 + 3915594018 x 5 1027579980 x 4 + 121010400 x 3 + 130837 x 765330 β 1 ( x ) = 1 88905600 201768035 x 9 778248135 x 8 + 1235314500 x 7 1037664180 x 6 + 491302224 x 5 127060920 x 4 + 14817600 x 3 84434 x 20940
 and  x  is given by the expression
x = t t n h
Equation (9) is solved for the independent solution at  f n + j  and  f n + κ  to obtain the continuous hybrid method
Y n + j = i = 0 1 ( j h ) i i   ! y n ( i ) + h 2 j = 0 1 σ j ( t ) f n + j + σ κ ( t ) f n + κ ,   κ = p , q , r , s , u , w
where  f n + j  and  f n + κ  have the following coefficients,
σ 0 ( x ) = 1 17280 352947 x 8 1613472 x 7 + 3092488 x 6 3226944 x 5 + 1990086 x 4 735392 x 3 + 156816 x 2 17280 x σ p ( x ) = 49 17280 50421 x 8 222264 x 7 + 404740 x 6 391608 x 5 + 214368 x 4 64224 x 3 + 8640 x 2 σ q ( x ) = 49 1920 16807 x 8 71344 x 7 + 123480 x 6 111328 x 5 + 55006 x 4 14064 x 3 + 1440 x 2 σ r ( x ) = 49 17280 252105 x 8 1029000 x 7 + 1694420 x 6 1433544 x 5 + 653520 x 4 151840 x 3 + 14400 x 2 σ s ( x ) = 49 17280 252105 x 8 987840 x 7 + 1550360 x 6 1241856 x 5 + 534450 x 4 118080 x 3 + 10800 x 2 σ u ( x ) = 49 1920 16807 x 8 63112 x 7 + 94668 x 6 72520 x 5 + 30016 x 4 6432 x 3 + 576 x 2 σ w ( x ) = 49 17280 50421 x 8 181104 x 7 + 260680 x 6 192864 x 5 + 77658 x 4 16304 x 3 + 1440 x 2 σ 1 ( x ) = 1 17280 352947 x 8 1210104 x 7 + 1680700 x 6 1210104 x 5 + 477456 x 4 98784 x 3 + 8640 x 2
Evaluating (12) at  t = p , q , r , s , u , w , 1  produces discrete hybrid scheme of the form (3) where
Y m ( i ) = y n + 1 7 ( i )   y n + 2 7 ( i )   y n + 3 7 ( i )   y n + 4 7 ( i )   y n + 5 7 ( i )   y n + 6 7 ( i )   y 1 ( i ) T y n ( i ) = y n 1 7 ( i )   y n 2 7 ( i )   y n 3 7 ( i )   y n 4 7 ( i )   y n 5 7 ( i )   y n 6 7 ( i )   y n ( i ) T F ( Y m ) = f n + 1 7   f n + 2 7   f n + 3 7   f n + 4 7   f n + 5 7   f n + 6 7   f 1 T f ( y n ) = f n 1 7     f n 2 7   f n 3 7   f n 4 7   f n 5 7   f n 6 7   f n T
A ( 0 )  is an identity matrix of dimension  7 × 7  given by
A ( 0 ) = 1   0   0   0   0   0   0 0   1   0   0   0   0   0 0   0   1   0   0   0   0 0   0   0   1   0   0   0 0   0   0   0   1   0   0 0   0   0   0   0   1   0 0   0   0   0   0   0   1
At  i = 0 , we obtain
e 0 = 0   0   0   0   0   0   1 0   0   0   0   0   0   1 0   0   0   0   0   0   1 0   0   0   0   0   0   1 0   0   0   0   0   0   1 0   0   0   0   0   0   1 0   0   0   0   0   0   1 ,   e 1 = 0   0   0   0   0   0   1 7 0   0   0   0   0   0   2 7 0   0   0   0   0   0   3 7 0   0   0   0   0   0   4 7 0   0   0   0   0   0   5 7 0   0   0   0   0   0   6 7 0   0   0   0   0   0     1 d 0 = 0   0   0   0   0   0   416173 88905600 0   0   0   0   0   0   14939 1389150 0   0   0   0   0   0   18399 1097600 0   0   0   0   0   0   15824 694575 0   0   0   0   0   0   102425 3556224 0   0   0   0   0   0   597 17150 0   0   0   0   0   0   10597 259200
b 0 = 33953 3175200 341699 29635200 105943 8890560 153761 17781120 943 231525 99359 88905600 6031 44452800 27821 694575 17 675 799 27783 5881 277830 2321 231525 1916 694575 233 694575 39141 548800 24111 1097600 369 7840 7299 219520 8613 548800 4737 1097600 9 17150 71152 694575 3832 231525 11344 138915 856 19845 4912 231525 4072 694575 496 694575 59375 444528 13375 1185408 210625 1778112 130625 3556224 25 864 26875 3556224 1625 1778112 1413 8575 54 8575 267 1715 99 3430 459 8575 9 1225 9 8575 25333 129600 49 86400 245 1296 833 51840 3283 43200 2989 259200 167 64800
At  i = 1 , we obtain
e 1 = 0   0   0   0   0   0   1 0   0   0   0   0   0   1 0   0   0   0   0   0   1 0   0   0   0   0   0   1 0   0   0   0   0   0   1 0   0   0   0   0   0   1 0   0   0   0   0   0   1 ,   d 1 = 0   0   0   0   0   0     751 17280 0   0   0   0   0   0       41 980   0   0   0   0   0   0       265 6272   0   0   0   0   0   0       278 6615     0   0   0   0   0   0       265 6272 0   0   0   0   0   0       41 980 0   0   0   0   0   0       751 17280
b 1 = 139849 846720 4511 31360 123133 846720 88547 846720 1537 31360 11351 846720 275 169344 1466 6615 71 2940 68 735 1927 26460 26 735 29 2940 8 6615 1359 6272 1377 31360 5927 31360 3033 31360 1377 31360 373 31360 9 6272 1448 6615 8 245 1784 6615 106 6615 8 245 64 6615 8 6615 36725 169344 775 18816 4625 18816 13625 169344 1895 18816 275 18816 275 169344 54 245 27 980 68 245 27 980 54 245 41 980 0 3577 17280 49 640 2989 17280 2989 17280 49 640 3577 17280 751 17280
The proposed hybrid method is thus given by,
y n + 1 7 = y n + 1 7 h y n + h 2 416173 88905600 f n + 33953 3175200 f n + 1 7 341699 29635200 f n + 2 7 + 105943 8890560 f n + 3 7 153761 17781120 f n + 4 7 + 943 231525 f n + 5 7 99359 88905600 f n + 6 7 + 6031 44452800 f n + 1 y n + 2 7 = y n + 2 7 h y n + h 2 14939 1389150 f n + 27821 694575 f n + 1 7 17 675 f n + 2 7 + 799 27783 f n + 3 7 5881 277830 f n + 4 7 + 2321 231525 f n + 5 7 1916 694575 f n + 6 7 + 233 694575 f n + 1 y n + 3 7 = y n + 3 7 h y n + h 2 18399 1097600 f n + 39141 54880 f n + 1 7 24111 1097600 f n + 2 7 + 369 7840 f n + 3 7 7299 219520 f n + 4 7 + 8613 548800 f n + 5 7 4737 1097600 f n + 6 7 + 9 17150 f n + 1 y n + 4 7 = y n + 4 7 h y n + h 2 15824 694575 f n + 71152 694575 f n + 1 7 3832 231525 f n + 2 7 + 11344 138915 f n + 3 7 856 19845 f n + 4 7 + 4912 231525 f n + 5 7 4072 694575 f n + 6 7 + 496 694575 f n + 1 y n + 5 7 = y n + 5 7 h y n + h 2 102425 3556224 f n + 59375 444528 f n + 1 7 13375 1185408 f n + 2 7 + 210625 1778112 f n + 3 7 130625 3556224 f n + 4 7 + 25 864 f n + 5 7 26875 3556224 f n + 6 7 + 1625 1778112 f n + 1 y n + 6 7 = y n + 6 7 h y n + h 2 597 17150 f n + 1413 8575 f n + 1 7 54 8575 f n + 2 7 + 267 1715 f n + 3 7 99 3430 f n + 4 7 + 459 8575 f n + 5 7 9 1225 f n + 6 7 + 9 8575 f n + 1 y n + 1 = y n + h y n + h 2 10597 259200 f n + 25333 129600 f n + 1 7 + 49 86400 f n + 2 7 + 245 1296 f n + 3 7 833 51840 f n + 4 7 + 3283 43200 f n + 5 7 + 2989 259200 f n + 6 7 + 167 64800 f n + 1
and
y n + 1 7 = y n + h 751 17280 f n + 139849 8467720 f n + 1 7 4511 31360 f n + 2 7 + 123133 846720 f n + 3 7 88547 846720 f n + 4 7 + 1537 31360 f n + 5 7 11351 846720 f n + 6 7 + 275 169344 f n + 1 y n + 2 7 = y n + h 41 980 f n + 1466 6615 f n + 1 7 71 2940 f n + 2 7 + 68 735 f n + 3 7 1927 26460 f n + 4 7 + 26 735 f n + 5 7 29 2940 f n + 6 7 + 8 6615 f n + 1 y n + 3 7 = y n + h 265 6272 f n + 1359 6272 f n + 1 7 + 1377 31360 f n + 2 7 + 5927 31360 f n + 3 7 3033 31360 f n + 4 7 + 1377 31360 f n + 5 7 373 31360 f n + 6 7 + 9 6272 f n + 1 y n + 4 7 = y n + h 278 6615 f n + 1448 6615 f n + 1 7 + 8 245 f n + 2 7 + 1784 6615 f n + 3 7 106 6615 f n + 4 7 + 8 245 f n + 5 7 64 6615 f n + 6 7 + 8 6615 f n + 1 y n + 5 7 = y n + h 265 6272 f n + 36725 169344 f n + 1 7 + 775 18816 f n + 2 7 + 4625 18816 f n + 3 7 + 13625 169344 f n + 4 7 + 1895 18816 f n + 5 7 275 18816 f n + 6 7 + 275 169344 f n + 1 y n + 6 7 = y n + h 41 980 f n + 54 245 f n + 1 7 + 27 980 f n + 2 7 + 68 245 f n + 3 7 + 27 980 f n + 4 7 + 54 245 f n + 5 7 + 41 980 f n + 6 7 + ( 0 ) f n + 1 y n + 1 = y n + h 751 17280 f n + 3577 17280 f n + 1 7 + 49 640 f n + 2 7 + 2989 17280 f n + 3 7 + 2989 17280 f n + 4 7 + 49 640 f n + 5 7 + 3577 17280 f n + 6 7 + 751 17280 f n + 1
Equations (14) and (15) are simultaneously applied as the proposed method for simulating second order oscillatory problems with applications to physical systems.

2.2. Implementation of the Method

The newly derived hybrid method is self-starting, unlike the conventional predictor-corrector methods. It is implemented in block mode over a non-overlapping one-step interval  t n ,     t n + 1 , with  n = 0 , 1 , 2 , 3 , , N 1 . Set data inputs like initial conditions, step size and differential equation to be solved. Thus at  n = 0 ,  the hybrid method is employed concurrently on the subinterval  t 0 , t 1  to get  y κ i , y κ i T i = 1 , 2 , 5 . Applying the values of the known former block, we get the values of the next subinterval  t 1 , t 2  and the process goes on in this manner until the final subinterval  t N 1 , t N  is determined.

3. Theoretical Analysis of the Method

Special properties of the proposed hybrid method will be analyzed theoretically in this section. These properties include zero-stability, consistency, convergence order, error constant and region of absolute stability. These are properties that a numerical method is expected to satisfy for such a method to be computationally reliable.

3.1. Order and Error Constant

The linear operator related to the newly derived hybrid method is established in this subsection.
Proposition 1.
The newly derived hybrid method has local truncation error  c 10 h 10 y ( 10 ) t n + O h 11 .
Proof. 
The linear difference operators associated with the hybrid method (14) and (15) are defined as
p y ( t n ) ; h = y t n + p h α τ 1 ( t ) y t n + τ 1 h + α τ 2 ( t ) y t n + τ 2 h + h 2 j = 0 1 β j ( t ) f n + j + β κ ( t ) f n + κ , κ = p , q , r , s , u , w q y ( t n ) ; h = y t n + q h α τ 1 ( t ) y t n + τ 1 h + α τ 2 ( t ) y t n + τ 2 h + h 2 j = 0 1 β j ( t ) f n + j + β κ ( t ) f n + κ , κ = p , q , r , s , u , w r y ( t n ) ; h = y t n + r h α τ 1 ( t ) y t n + τ 1 h + α τ 2 ( t ) y t n + τ 2 h + h 2 j = 0 1 β j ( t ) f n + j + β κ ( t ) f n + κ , κ = p , q , r , s , u , w s y ( t n ) ; h = y t n + s h α τ 1 ( t ) y t n + τ 1 h + α τ 2 ( t ) y t n + τ 2 h + h 2 j = 0 1 β j ( t ) f n + j + β κ ( t ) f n + κ , κ = p , q , r , s , u , w u y ( t n ) ; h = y t n + u h α τ 1 ( t ) y t n + τ 1 h + α τ 2 ( t ) y t n + τ 2 h + h 2 j = 0 1 β j ( t ) f n + j + β κ ( t ) f n + κ , κ = p , q , r , s , u , w w y ( t n ) ; h = y t n + w h α τ 1 ( t ) y t n + τ 1 h + α τ 2 ( t ) y t n + τ 2 h + h 2 j = 0 1 β j ( t ) f n + j + β κ ( t ) f n + κ , κ = p , q , r , s , u , w 1 y ( t n ) ; h = y t n + h α τ 1 ( t ) y t n + τ 1 h + α τ 2 ( t ) y t n + τ 2 h + h 2 j = 0 1 β j ( t ) f n + j + β κ ( t ) f n + κ , κ = p , q , r , s , u , w
Assuming  y ( t )  is sufficiently differentiable; Equation (16) is expanded in power of  h  with the aid of Taylor series. It is important to state that the first non-zero term of each formula in Equation (16) is  c 10 h 10 y ( 10 ) t n + O h 11 . The same procedure applies to the derivative scheme in Equation (15). 
Definition 2.
Ref. [43]
“A linear multistep method for a second order problem is of order  p  if it satisfies   c 0  =   c 1  =   c 2  =   c 3  = … =   c p  =   c p + 1  =   0 ,   c p + 2 0 , where,
c 0 = j = 0 k α j c 1 = j = 0 k j α j β j   .   .   . c p = j = 0 k 1 p ! j p α j 1 ( p 1 ) ! j p 1 β j ,   p = 2 , 3 , , q + 1
The parameter   c p + 2 0  is referred to as the error constant with the local truncation error defined as   t n + k = c p + 2 h p + 2 y ( p + 2 ) ( t n ) + O h p + 3 ”.
The hybrid method is expanded about the point  t n  via Taylor series as follows
j = 0 1 7 h j j ! y n j y n 1 7 h y n 416173 88905600 h 2 y n j = 0 h j + 2 j ! y n j + 2 33953 3175200 1 7 j 341699 29635200 2 7 j + 105943 8890560 3 7 j 153761 17781120 4 7 j + 943 231525 5 7 j 99359 88905600 6 7 j + 6031 44452800 1 j j = 0 2 7 h j j ! y n j y n 2 7 h y n 14939 1389150 h 2 y n j = 0 h j + 2 j ! y n j + 2 27821 694575 1 7 j 17 675 2 7 j + 799 27783 3 7 j 5881 277830 4 7 j + 2321 231525 5 7 j 1916 694575 6 7 j + 233 694575 1 j j = 0 3 7 h j j ! y n j y n 3 7 h y n 18399 1097600 h 2 y n j = 0 h j + 2 j ! y n j + 2 39141 54880 1 7 j 24111 1097600 2 7 j + 369 7840 3 7 j 7299 219520 4 7 j + 8613 548800 5 7 j 4737 1097600 6 7 j + 9 17150 1 j j = 0 4 7 h j j ! y n j y n 4 7 h y n 15824 694575 h 2 y n j = 0 h j + 2 j ! y n j + 2 71152 694575 1 7 j 3832 231525 2 7 j + 11344 138915 3 7 j 856 19845 4 7 j + 4912 231525 5 7 j 4072 694575 6 7 j + 496 694575 1 j j = 0 5 7 h j j ! y n j y n 5 7 h y n 102425 3556224 h 2 y n j = 0 h j + 2 j ! y n j + 2 59375 444528 1 7 j 13375 1185408 2 7 j + 210625 1778112 3 7 j 130625 3556224 4 7 j + 25 864 5 7 j 26875 3556224 6 7 j + 1625 1778112 1 j j = 0 6 7 h j j ! y n j y n 6 7 h y n 597 17150 h 2 y n j = 0 h j + 2 j ! y n j + 2 1413 8575 1 7 j 54 8575 2 7 j + 267 1715 3 7 j 99 3430 4 7 j + 459 8575 5 7 j 9 1225 6 7 j + 9 8575 1 j j = 0 h j j ! y n j y n h y n 10597 259200 h 2 y n j = 0 h j + 2 j ! y n j + 2 25333 129600 1 7 j + 49 86400 2 7 j + 245 1296 3 7 j 833 51840 4 7 j + 3283 43200 5 7 j + 2989 259200 6 7 j + 167 64800 1 j = 0 0 0 0 0 0 0
Corollary 1.
The local truncation error of the hybrid method is
Local   truncation   error = 1.9603 × 10 11 h 10 y ( 10 ) t n + O h 11 4.8794 × 10 11 h 10 y ( 10 ) t n + O h 11 7.6453 × 10 11 h 10 y ( 10 ) t n + O h 11   1.0439 × 10 10 h 10 y ( 10 ) t n + O h 11 1.3262 × 10 10 h 10 y ( 10 ) t n + O h 11 1.5931 × 10 10 h 10 y ( 10 ) t n + O h 11 1.9558 × 10 10 h 10 y ( 10 ) t n + O h 11    
Table 1 summarizes the order as well as error constant of the hybrid method.
Therefore, the hybrid method is of order eight.

3.2. Consistency

A continuous linear multistep method is consistent if it has order  p 1  [44].
Thus, the new hybrid method is consistent because its order  p = 8 .

3.3. Zero-Stability

Definition 3.
Ref. [45]
“A method is zero-stable, if the first characteristic polynomial   ρ ( z )  with roots   z s ,   s = 1 , 2 , , k  given by   ρ ( z ) = det ( z A ( 0 ) e 0 )  satisfies   z s 1  and every root satisfying   z s = 1  has multiplicity not more than the order of the differential equation”.
The first characteristic polynomial of the new hybrid method is given as
ρ ( z ) = z 1       0       0       0       0       0       0 0       1       0       0       0       0       0 0       0       1       0       0       0       0 0       0       0       1       0       0       0 0       0       0       0       1       0       0 0       0       0       0       0       1       0 0       0       0       0       0       0       1 0         0       0       0       0       0       1 0         0       0       0       0       0       1 0         0       0       0       0       0       1 0         0       0       0       0       0       1 0         0       0       0       0       0       1 0         0       0       0       0       0       1 0         0       0       0       0       0       1
Evaluating Equation (20), we obtain
z 6 ( z 1 ) = 0
Solving Equation (21) for  z  gives  z 1 = z 2 = z 3 = z 4 = z 5 = z 6 = 0  and  z 7 = 1 . The hybrid method is therefore zero-stable.

3.4. Convergence

The convergence of a method is analyzed by employing Theorem 2.
Theorem 2.
Ref. [45]
Consistency and zero-stability are necessary and sufficient conditions a linear multistep method must satisfy for it to be convergent”.
Therefore, the hybrid method is convergent because it is consistent and zero-stable.

3.5. Region of Absolute Stability

Definition 4.
Ref. [45]
“The region of absolute stability of a numerical method is the set of complex values   λ h  for which all solutions of the test problem  y = λ 2 y  will remain bounded as   n ”.
The boundary locus method is adopted in generating the stability polynomial of the hybrid method. The polynomial is
h ¯ ( w ) = h 14 1 195328244980512 w 7 + 383 39065648996102404266493378560 w 6                               h 12 64313 39065648996102400 w 7 + 14662561 117196946988307200 w 6                               h 10 21601 256261545892800 w 7 + 54664709 1025046183571200 w 6                               + h 8 2021 1394620657920 w 7 36041177 3486551644800 w 6 h 6 757 1778852880 w 7 + 430601 444713220 w 6                             h 4 5077 123480 w 6 23 246960 w 7 h 2 2 147 w 7 + 94 147 w 6 + w 7 2 w 6
Figure 1 shows the plot of stability region of the proposed hybrid method.
From Figure 1, it is clear the hybrid method is A–stable. The stability region is the exterior of the curve.

4. Numerical Results and Discussion

We shall simulate some second order oscillatory problems with applications to physical systems using the newly derived hybrid method.

4.1. Problem 1 (Highly Stiff Oscillatory System)

Consider the highly stiff system of oscillatory problem
y 1 ( t ) y 2 ( t ) +       13     12 12               13 y 1 ( t ) y 2 ( t ) = 12 ε 2 5       3   2 2             3 y 1 ( t ) y 2 ( t ) + ε 2 f 1 ( t ) f 2 ( t ) ,   y 1 ( 0 ) y 2 ( 0 ) = ε 2 ε 2 ,   y 1 ( 0 ) y 2 ( 0 ) = 4       6
where  ε = 10 3 f 1 ( t ) = 36 / 5 sin t + 24 sin 5 t  and  f 2 ( t ) = 24 / 5 sin t 36 sin 5 t . The analytical solution is given as
y 1 ( t ) y 2 ( t ) = sin t sin 5 t + ε cos t sin t + sin 5 t + ε cos t
The results obtained in Table 2 clearly showed that the proposed hybrid method (which is of order eight) outperformed the block hybrid method of order eleven developed by [41]. The solution curves in Figure 2 also showed that the approximate solution obtained by the hybrid method converges to the analytical solution on the interval [0, 100]. Notice that the global error at  y 1  and  y 2  are same. This is because the exact solution at  y 1  and  y 2  are almost the same.

4.2. Problem 2 (Projectile Motion with Aerodynamic Drag)

The projectile motion (the X and Y displacements) of an object is given by the differential equations
m x + c d x x 2 + y 2 ( 1 / 2 ) = 0
m y + c d y x 2 + y 2 ( 1 / 2 ) = w
The differential equations were solved using the following coefficients of aerodynamic drag  c d = 0.1 ,   0.3 ,   0.5   N / m . Take  m = 10   kg x ( 0 ) = 100   ms 1 y ( 0 ) = 10   ms 1 x ( 0 ) = 0 y ( 0 ) = 0 g = 9.81   ms 2  and  w = m g , see [42].
The authors [6] simulated the problem above by converting the equations into their equivalent first order systems and then employed ode45 solver to generate their solution curves.
It is obvious from the plots in Figure 3 that as  c d  increases, the displacement’s maximum value decreases in the  y  direction. In the  x  direction, the value of the displacement for which the displacement in the  y  direction is maximum, also decreases. The solution curves generated using the new hybrid method clearly converges to those of the curves generated by the ode45 solver at different values of the coefficients of aerodynamics  c d . This shows that the hybrid method is computationally reliable.

4.3. Problem 3 (Half Cylinder Rolling on a Horizontal Plane)

The governing equation for motion of a half cylinder (or semi cylinder) with mass m and radius  R  rolling on a horizontal plane is given by the oscillatory problem,
3 m R 2 2 2 m b R sin θ θ m b R cos θ θ 2 = w b cos θ
The differential Equation (27) is solved by generating a solution curve of angle of oscillation  θ  against time  t . The mass center is at  G  located at distance  b = 4 R 3 π m  from the cylindrical axis at  O . Take  m = 5 kg ,   R = 0.1   m g = 9.81   ms 2  and  w = m g . Figure 4 shows a typical half-cylinder on a horizontal plane with an angle of inclination  θ .
The authors [6] also simulated the problem above by converting the equation into its equivalent first order systems and then employed ode45 solver to generate the solution curve.
The solution curves obtained by using the new hybrid method on Problem 3 converge to those of the curves generated by the ode45 solver showing that the new hybrid method is computationally reliable, see Figure 5.

4.4. Problem 4 (The Stiefel and Bettis Oscillatory Problem)

Consider the Stiefel and Bettis oscillatory system of equations,
y 1 ( t ) + y 1 ( t ) = 0.001 cos ( t ) ,   y 1 ( 0 ) = 1 ,   y 1 ( 0 ) = 0 y 2 ( t ) + y 2 ( t ) = 0.001 sin ( t ) ,   y 2 ( 0 ) = 1 ,   y 2 ( 0 ) = 0.9995
The exact solution is given by,
y 1 ( t ) = 0.0005 t sin   ( t ) + cos ( t ) y 2 ( t ) = 0.0005 t cos   ( t ) + sin ( t )

4.5. Problem 5 (Damped Duffing Oscillatory Problem)

Consider the damped Duffing oscillatory problem,
y ( t ) + y ( t ) + y ( t ) + y 2 ( t ) y ( t ) = 2 cos ( t ) cos 3 ( t ) ,     y ( 0 ) = 0 ,     y   ( 0 ) = 1
defined on the  0 t 100 . The analytical solution is given by,
y ( t ) = sin   ( t )

4.6. Problem 6 (Undamped Duffing Oscillatory Problem)

Consider the following undamped Duffing oscillatory problem,
y ( t ) + y ( t ) + y 3 ( t ) = ξ cos ( λ t ) ,   y ( 0 ) = ϕ 0 ,     y   ( 0 ) = 0
where,
α = 2.00426728067 × 10 1 ,   ξ = 2.0 × 10 3 ,     λ = 1.01
The analytical solution to the problem is
y ( t ) = i = 0 3 ϕ 2 i + 1 C o s 2 i + 1 λ t
where,
ϕ 1 ,   ϕ 3 ,   ϕ 5 ,   ϕ 7 = 2.00179477536 × 10 1 ,   2.4946143 × 10 4 ,   3.04014 × 10 7 ,       3.74 × 10 10
Table 3, Table 4 and Table 5 display the performance of the new hybrid method on Problems 4–6 in comparison to other existing methods. The numerical results clearly showed that the new hybrid method performed better than the methods we juxtaposed our results with in terms of the number of steps taken, absolute errors, maximum errors, average errors and computation time. In Figure 5 for instance, it is obvious that the approximate solution obtained using the hybrid method agree with the solution obtained using ode45 solver for Problem 3. For Problem 4, it is clear that the approximate solution using the hybrid method converges to the exact solution, see Figure 6. Thus, the new hybrid method outperformed the methods we juxtaposed our results with. The execution time (in seconds) further showed that the new hybrid method is more efficient in terms of the economy of computer time (and also more accurate than the OSHBM developed by [38]) for Problems 4 and 6, see Table 3 and Table 5 respectively. Table 4 also shows that the hybrid method is more accurate than the VOVSM developed by [38] in terms of maximum and average errors. Also, the hybrid method requires fewer numbers of steps than the VOVSM. Figure 7 and Figure 8 further buttress the fact that the hybrid method is more accurate for Problems 5 and 6.

4.7. Problem 7 (The Van der Pol Model)

The Van der Pol equation, a nonlinear stiff equation is given by,
y ( t ) + μ 1 + y 2 ( t ) y ( t ) + y ( t ) = 0
This nonlinear model does not have a closed form solution, thus, we need to obtain an approximate solution using the new hybrid method. Firstly, Equation (36) is transformed to an equivalent system of equations of the form
y 1 ( t ) = y 2 ( t ) ,     y 1 ( 0 ) = 2 y 2 ( t ) = y 1 ( t ) + ( 1 y 1 2 ( t ) ) y 2 ( t ) / ,       y 2 ( 0 ) = 0
defined on the interval [0, 70]. This transformation is achieved by substituting  y 1 ( t ) = φ ( t ) ,   y 2 ( t ) = μ φ ( t )  and  t = t / μ , where  = 1 / μ 2  is a parameter that controls stiffness. In this paper, we assume  = 500 .
Since the Van der Pol model does not have a closed form solution, the approximate solution using the new hybrid method was compared with those of FPHBI developed by [36] and ode 15s at the end points  x = 1 x = 5 x = 10  and  x = 20 . The numerical solution and accuracy curves generated clearly showed that the new hybrid method is accurate and computationally reliable, see Table 6 and Figure 9 respectively.
It is therefore obvious that the proposed hybrid method is accurate and efficient. The research works earlier carried out by [46,47,48,49,50,51,52] on the derivations, implementations and applications of numerical methods on physical systems may also be consulted for more details.

5. Conclusions

In this research, a hybrid method was formulated for the simulations of second order oscillatory problems with applications to physical systems. The results clearly showed that the new method is efficient, accurate and reliable. The properties of the method were also analyzed and it is obvious that the method is zero-stable, consistent and convergent. The region of absolute stability of the method also showed that the method is A-stable, thus, making it suitable for the solutions and simulations of problems of the form (1). In view of the foregoing, the proposed hybrid method is therefore recommended for the simulations and solutions of physical systems of second order. Further research may also extend this work, by considering the possibility of proposing new methods for simulating larger systems of equations.

Author Contributions

Conceptualization, L.J.K., J.S., J.N.N., A.S. and Y.W.; methodology, L.J.K., J.S. and J.N.N.; software, J.S., J.N.N. and L.J.K.; validation, A.S. and Y.W.; formal analysis, A.S., J.N.N. and Y.W.; writing—original draft preparation, L.J.K., J.S. and J.N.N.; writing—review and editing, A.S., Y.W. and J.N.N.; supervision, J.S. and J.N.N.; project administration, A.S., J.N.N. and J.S. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the National Natural Science Foundation of China (Grant No.: 12171435).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Acknowledgments

All authors would like to thank the Reviewers and Editors for all the constructive comments, observations and suggestions, which helped to improve the quality of the manuscript.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

The following abbreviations are used in Table 2, Table 3, Table 4, Table 5 and Table 6.
  t evaluation point
TOLtolerance level
ABSERRabsolute error
AVERRaverage error
MAXERRmaximum error
APPSOLapproximate solution
VOVSMvariable order variable step size method proposed by [35]
FPHBIfour-point hybrid block integrator by [36]
OSHBMone-step hybrid method developed by [38]
BHMblock hybrid method of order 11 developed by [41]
HMthe proposed hybrid method in Equations (14) and (15)
Timeexecution time in seconds

References

  1. Musaev, H.K. The Cauchy problem for degenerate parabolic convolution equation. TWMS J. Pure Appl. Math. 2021, 12, 278–288. [Google Scholar]
  2. Wend, D.V.V. Uniqueness of solution of ordinary differential equations. Am. Math. Mon. 1967, 74, 27–33. [Google Scholar] [CrossRef]
  3. Raisinghania, M.D. Ordinary and Partial Differential Equations; S. Revised edition; Chand and Company LTD: New-Delhi, India, 2014. [Google Scholar]
  4. Abolfazl, J.; Hadi, F. The application of Duffing oscillator in weak signal detection. ECTI Trans. Electr. Eng. Electron. Commun. 2011, 9, 1–6. [Google Scholar]
  5. Soraya, H. Pendulum with aerodynamic and viscous damping. J. Appl. Inf. Commun. Technol. 2016, 3, 43–47. [Google Scholar] [CrossRef]
  6. Harihara, P.; Childs, D.N. Solving Problems in Dynamics and Vibration Using MATLAB; Department of Mechanical Engineering, Texas A and M University College Station: College Station, TX, USA, 2020; pp. 1–104. [Google Scholar]
  7. Zhang, W.B. Differential Equations, Bifurcations, and Chaos in Economics; World Scientific Publishing: Hackensack, NJ, USA, 2005. [Google Scholar]
  8. Sunday, J.; Chigozie, C.; Omole, E.O.; Gwong, J.B. A pair of three-step hybrid block methods for the solutions of linear and nonlinear first-order systems. Util. Math. 2021, 118, 1–15. [Google Scholar] [CrossRef]
  9. Adiguzel, R.S.; Aksoy, U.; Karapinar, E.; Erthan, I.M. On the solutions of fractional differential equations via Geraghty type hybrid contractions. Appl. Comput. Math. 2021, 20, 313–333. [Google Scholar]
  10. Aliev, F.A.; Aliyev, N.A.; Hajiyeva, N.S.; Mahmudov, N.I. Some mathematical problems and their solutions for the oscillating systems with liquid dampers: A review. Appl. Comput. Math. 2021, 20, 339–365. [Google Scholar]
  11. Shokri, A.; Tahmourasi, M. A new two-step Obrechkoff method with vanished phase-lag and some of its derivatives for the numerical solution of radial Schrödinger equation and related IVPs with oscillating solutions. Iran. J. Math. Chem. 2017, 8, 137–159. [Google Scholar]
  12. Marian, D.; Ciplea, S.A.; Lungu, N. Ulam-Hyers stability of Darboux-Lonescu problem. Carpathian J. Math. 2021, 37, 211–216. [Google Scholar] [CrossRef]
  13. Adeyefa, E.O.; Omole, E.O.; Shokri, A. Numerical solution of second-order nonlinear partial differential equations originating from physical phenomena using Hermite based block methods. Results Phys. 2023, 46, 106270. [Google Scholar] [CrossRef]
  14. Noor, M.A.; Noor, K.I. Some new classes of strongly generalized preinvex functions. TWMS J. Pure Appl. Math. 2021, 12, 181–192. [Google Scholar]
  15. Akbarov, S.D.; Mehdiyev, M.A.; Zeynalov, A.M. Dynamic of the moving ring-load acting in the interior of the bi-layered hollow cylinder with imperfect contact between the layers. TWMS J. Pure Appl. Math. 2021, 12, 223–242. [Google Scholar]
  16. Tongxing, L.; Yuriy, V.R.; Shuhong, T. Oscillation of second-order nonlinear differential equations with damping. Math. Slovaca 2014, 64, 1227–1236. [Google Scholar]
  17. Juraev, D.A. On the Cauchy problem for matrix factorizations of the Helmholtz equation in an unbounded domain in R2. Sib. Electron. Math. Rep. 2018, 15, 1865–1877. [Google Scholar]
  18. Marian, D.; Ciplea, S.A.; Lungu, N. Ulam-Hyers stability of Euler’s equation in the calculus of variations. Mathematics 2021, 9, 3320. [Google Scholar] [CrossRef]
  19. Juraev, D.A. The Cauchy problem for matrix factorization of the Helmholtz equation in an unbounded domain. Sib. Electron. Math. Rep. 2017, 14, 752–764. [Google Scholar]
  20. Larin, V.B.; Tunik, A.A. Two-stage procedure of H∞-parameterization of stabilizing controllers applied to quadrotor flight control. TWMS J. Pure Appl. Math. 2021, 12, 199–208. [Google Scholar]
  21. Davvaz, B.; Mohammadi, F. Different types of ideals and homomorphisms of (m, n)-semirings. TWMS J. Pure Appl. Math. 2021, 12, 209–222. [Google Scholar]
  22. Pankov, P.S.; Zheentaeva, Z.K.; Shirinov, T. Asymptotic reduction of solution space dimension for dynamical systems. TWMS J. Pure Appl. Math. 2021, 12, 243–253. [Google Scholar]
  23. Chen, W.; Wu, W.X.; Teng, Z.D. Complete dynamics in a nonlocal dispersal two-strain SIV epidemic model with vaccination and latent delays. Appl. Comput. Math. 2020, 19, 360–391. [Google Scholar]
  24. Yusufoglu, E. Numerical solution of Duffing equation by the Laplace decomposition algorithm. Appl. Math. Comput. 2006, 177, 572–580. [Google Scholar]
  25. Iskandarov, S.; Komartsova, E. On the influence of integral perturbations on the boundedness of solutions of a fourth-order linear differential equation. TWMS J. Pure Appl. Math. 2022, 13, 3–9. [Google Scholar]
  26. Hamidov, S.I. Optimal trajectories in reproduction models of economic dynamics. TWMS J. Pure Appl. Math. 2022, 13, 16–24. [Google Scholar]
  27. Nourazar, S.; Mirzabeigy, A. Approximate solution for nonlinear Duffing oscillator with damping effect using the modified differential transform method. Sci. Iran. B 2013, 20, 364–368. [Google Scholar]
  28. Kalsoom, H.U.; Ali, M.A.; Abbas, M.U.; Budak, H.Ü.; Murtaza, G.H. Generalized quantum Montgomery identity and Ostrowski type inequalities for preinvex functions. TWMS J. Pure Appl. Math. 2022, 13, 72–90. [Google Scholar]
  29. He, J.H. Variational iteration method. A kind of nonlinear analytical technique. Int. J. Nonlinear Mech. 1999, 34, 699–708. [Google Scholar] [CrossRef]
  30. He, J.H. Variational iteration method for autonomous ordinary differential systems. Appl. Math. Comput. 2000, 114, 115–123. [Google Scholar] [CrossRef]
  31. Goharee, F.; Babolian, E. Modified variational iteration method for solving Duffing equations. Indian J. Sci. Res. 2014, 6, 25–29. [Google Scholar]
  32. Agocs, F.J.; Handley, W.J.; Lasenby, A.N.; Hobson, M.P. Efficient method for solving highly oscillatory ordinary differential equations with applications to physical systems. Phys. Rev. Res. 2020, 2, 013030. [Google Scholar] [CrossRef] [Green Version]
  33. Ukpebor, L.A.; Omole, E.O. Three-step optimized block backward differentiation formulae (TOBBDF) for solving ordinary differential equations. Afr. J. Math. Comput. Sci. Res. 2020, 13, 51–57. [Google Scholar] [CrossRef] [Green Version]
  34. Obarhua, F.O.; Kayode, S.J. Continuous explicit hybrid method for solving second order ordinary differential equations. Pure Appl. Math. J. 2020, 9, 26–31. [Google Scholar] [CrossRef]
  35. Rasedee, A.F.N.; Ishak, N.; Hamzah, S.T.; Ijam, H.M.; Suleiman, M.; Ibrahim, Z.B.; Sathar, M.H.A.; Ramli, N.A.; Kamaruddin, N.S. Variable order variable step-size algorithm for solving nonlinear Duffing oscillator. J. Phys. Conf. Ser. 2017, 890, 012045. [Google Scholar] [CrossRef] [Green Version]
  36. Sunday, J.; Kumleng, G.M.; Kamoh, N.M.; Kwanamu, J.A.; Skwame, Y.; Sarjiyus, O. Implicit four-point hybrid block integrator for the simulations of stiff models. J. Niger. Soc. Phys. Sci. 2022, 4, 287–296. [Google Scholar] [CrossRef]
  37. Khashan, M.M.; Syam, M.I.; Alomari, A.K. Numerical solution of cubic free undamped Duffing oscillator equation using continuous implicit hybrid method. Ital. J. Pure Appl. Math. 2019, 42, 588–597. [Google Scholar]
  38. Kwari, L.J.; Sunday, J.; Ndam, J.N.; James, A.A. On the numerical approximation and simulation of damped and undamped Duffing oscillators. Sci. Forum J. Pure Appl. Sci. 2021, 21, 503–515. [Google Scholar] [CrossRef]
  39. Yakubu, S.D.; Yahaya, Y.A.; Lawal, K.O. Three-step block hybrid linear multistep method for solution of special second order ordinary differential equations. J. Niger. Math. Soc. 2021, 40, 149–160. [Google Scholar]
  40. Awari, Y.S. Some generalized two-step block hybrid Numerov method for solving general second order ordinary differentialequations without predictors. Sci. World J. 2015, 12, 12–18. [Google Scholar]
  41. Jator, S.N.; King, K.L. Integrating oscillatory general second-order initial value problems using a block hybrid method of order 11. Math. Probl. Eng. 2018, 2018, 3750274. [Google Scholar] [CrossRef]
  42. Sunday, J.; Ndam, J.N.; Kwari, L.J. An accuracy-preserving block hybrid algorithm for the integration of second-order physical systems with oscillatory solutions. J. Niger. Soc. Phys. Sci. 2023, 5, 1017. [Google Scholar] [CrossRef]
  43. Lambert, J.D. Numerical Methods for Ordinary Differential Systems: The Initial Value Problem; John Wiley and Sons Ltd.: Hoboken, NJ, USA, 1991. [Google Scholar]
  44. Lambert, J.D. Computational Methods in Ordinary Differential Equations; John Wiley & Sons, Inc.: New York, NY, USA, 1973. [Google Scholar]
  45. Fatunla, S.O. Numerical integrators for stiff and highly oscillatory differential equations. Math. Comput. 1980, 34, 373–390. [Google Scholar] [CrossRef]
  46. Dahlquist, G.G. Convergence and stability in the numerical integration of ordinary differential equations. Math. Scand. 1956, 4, 33–50. [Google Scholar] [CrossRef] [Green Version]
  47. Shokri, A.; Mehdizadeh Khalsaraei, M.; Molayi, M. Nonstandard dynamically consistent numerical methods for MSEIR model. J. Appl. Comput. Mech. 2022, 8, 196–205. [Google Scholar]
  48. Sunday, J.; Shokri, A.; Akinola, R.O.; Joshua, K.V.; Nonlaopon, K. A convergence-preserving non-standard finite difference scheme for the solutions of singular Lane-Emden equations. Results Phys. 2022, 42, 106031. [Google Scholar] [CrossRef]
  49. Rufai, M.A.; Shokri, A.; Omole, E.O. A one-point third-derivative hybrid technique for solving second-order oscillatory and periodic problems. J. Math. 2023, 2023, 2343215. [Google Scholar] [CrossRef]
  50. Kuboye, J.O.; Omar, Z. Derivation of a six-step block method for direct solution of second order ordinary differential equations. Math. Comput. Appl. 2015, 20, 151–159. [Google Scholar]
  51. Sunday, J.; Shokri, A.; Kwanamu, J.A.; Nonlaopon, K. Numerical integration of stiff differential systems using non-fixed step-size strategy. Symmetry 2022, 14, 1575. [Google Scholar] [CrossRef]
  52. Argyros, I.K.; Santhosh, G. Extended Kung-Traub-type method for solving equations. TWMS J. Pure Appl. Math. 2021, 12, 193–198. [Google Scholar]
Figure 1. Region of absolute stability of the hybrid method.
Figure 1. Region of absolute stability of the hybrid method.
Axioms 12 00282 g001
Figure 2. Solution curves for Problem 1 on the interval [0, 100].
Figure 2. Solution curves for Problem 1 on the interval [0, 100].
Axioms 12 00282 g002
Figure 3. Solution curves for Problem 2 showing the  X  and  Y  displacements on the interval [0, 100].
Figure 3. Solution curves for Problem 2 showing the  X  and  Y  displacements on the interval [0, 100].
Axioms 12 00282 g003
Figure 4. Schema of a half cylinder rolling on a horizontal plane.
Figure 4. Schema of a half cylinder rolling on a horizontal plane.
Axioms 12 00282 g004
Figure 5. Solution Curves for Problem 3 showing the angle of inclination against time.
Figure 5. Solution Curves for Problem 3 showing the angle of inclination against time.
Axioms 12 00282 g005
Figure 6. Solution curves for Problem 4 on the interval [0, 80].
Figure 6. Solution curves for Problem 4 on the interval [0, 80].
Axioms 12 00282 g006
Figure 7. Accuracy curves for Problem 5.
Figure 7. Accuracy curves for Problem 5.
Axioms 12 00282 g007
Figure 8. Accuracy curves for Problem 6.
Figure 8. Accuracy curves for Problem 6.
Axioms 12 00282 g008
Figure 9. Solution curves for Problem 7 on the interval [0, 70].
Figure 9. Solution curves for Problem 7 on the interval [0, 70].
Axioms 12 00282 g009
Table 1. Order and error constant of the hybrid method.
Table 1. Order and error constant of the hybrid method.
  y t n + j OrderError Constant
  j = 1 / 7   8   1.9603 × 10 11
  j = 2 / 7   8   4.8794 × 10 11
  j = 3 / 7   8   7.6453 × 10 11
  j = 4 / 7   8   1.0439 × 10 10
  j = 5 / 7   8   1.3262 × 10 10
  j = 6 / 7   8   1.5931 × 10 10
  j = 1   8   1.9558 × 10 10
Table 2. Comparison of global errors for Problem 1 using our method (HM) and BHM on the interval [0, 100].
Table 2. Comparison of global errors for Problem 1 using our method (HM) and BHM on the interval [0, 100].
  h   HM   y 1   BHM   y 1   HM   y 2   BHM   y 2
  1 2.04 × 10−21.57 × 10−12.04 × 10−21.57 × 10−1
  1 / 2 1.21 × 10−79.38 × 10−51.21 × 10−79.38 × 10−5
  1 / 4 1.57 × 10−95.29 × 10−81.57 × 10−95.29 × 10−8
  1 / 8 3.14 × 10−131.54 × 10−113.14 × 10−131.54 × 10−11
  1 / 16 6.89 × 10−175.07 × 10−146.89 × 10−175.07 × 10−14
Table 3. Comparison of the absolute errors and execution time for Problem 4 using our method (HM) and OSHBM.
Table 3. Comparison of the absolute errors and execution time for Problem 4 using our method (HM) and OSHBM.
t HM   y 1 HM   y 2 OSHBM   y 1 OSHBM y 2 Time in HMTime in OSHBM
0.10006.944134 × 10−142.053113 × 10−141.256535 × 10−122.826929 × 10−120.03210.4261
0.20007.080203 × 10−142.759703 × 10−142.113954 × 10−125.899430 × 10−120.04930.4362
0.30007.764101 × 10−145.796163 × 10−142.376410 × 10−126.830921 × 10−120.06670.4621
0.40007.888446 × 10−146.826273 × 10−143.424150 × 10−121.499123 × 10−120.10560.4891
0.50007.810730 × 10−146.264544 × 10−143.394396 × 10−121.839452 × 10−120.12290.5001
0.60007.678613 × 10−146.505462 × 10−143.343548 × 10−121.655884 × 10−110.14190.5201
0.70007.554268 × 10−146.540990 × 10−144.294920 × 10−121.247034 × 10−110.15940.5578
0.80007.457679 × 10−148.378897 × 10−144.257394 × 10−128.431255 × 10−110.17680.5881
0.90007.397727 × 10−148.072475 × 10−145.234413 × 10−125.323966 × 10−110.19390.6101
1.00007.377743 × 10−148.038814 × 10−146.226530 × 10−123.212587 × 10−110.28850.6312
Table 4. Comparison of the number of steps, maximum error and average errors for Problem 5 using our method (HM) and VOVSM.
Table 4. Comparison of the number of steps, maximum error and average errors for Problem 5 using our method (HM) and VOVSM.
TOLMethodStepsMAXERRAVERR
  10 2 HM2092.12651 × 10−54.81903 × 10−6
VOVSM2171.07600 × 10−13.03894 × 10−2
  10 4 HM2643.19061 × 10−71.41210 × 10−8
VOVSM2841.24649 × 10−31.92899 × 10−4
  10 6 HM3161.29139 × 10−91.09310 × 10−9
VOVSM3301.28039 × 10−53.26641 × 10−6
  10 8 HM4864.90431 × 10−107.19038 × 10−10
VOVSM4997.27324 × 10−71.46368 × 10−7
  10 10 HM7404.34750 × 10−134.12474 × 10−13
VOVSM7727.83863 × 10−92.03721 × 10−9
Table 5. Comparison of the end-point absolute errors and execution time for Problem 6 using our method (HM) and OSHBM.
Table 5. Comparison of the end-point absolute errors and execution time for Problem 6 using our method (HM) and OSHBM.
  h ABSERRABSERRTimeTime
in HMin OSHBMin HMin OSHBM
  m / 500 6.000145 × 10−172.256151 × 10−161.26762.0032
  m / 1000 4.010991 × 10−171.678127 × 10−161.62722.2001
  m / 2000 2.102903 × 10−176.061690 × 10−172.12873.1245
  m / 3000 1.211451 × 10−175.167032 × 10−172.27813.6234
  m / 4000 3.190380 × 10−184.267817 × 10−173.01824.0151
  m / 5000 2.201181 × 10−182.762717 × 10−174.01676.1930
Note that  m = 10  in Table 5.
Table 6. Comparison of approximate solutions for Problem 7 at  h = 0.1 .
Table 6. Comparison of approximate solutions for Problem 7 at  h = 0.1 .
  t   y i APPSOL (HM)APPSOL (FPHBI)APPSOL (ode 15s)
  1   y 1   1.8650951120   1.8650950931   1.8650950571
  y 2   0.7524845401       0.7524845342   0.7524845299
  5   y 1   1.8985234609   1.8985234584   1.8985234421
  y 2   0.7289532599   0.7289532571   0.7289532451
  10   y 1   1.7865365312   1.7865365214   1.7865365103
  y 2   0.8156276678   0.8156276595   0.8156276438
  20   y 1   1.5075643427   1.5075643297   1.5075643177
  y 2   1.1911230055   1.1911230041   1.1911230003
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Kwari, L.J.; Sunday, J.; Ndam, J.N.; Shokri, A.; Wang, Y. On the Simulations of Second-Order Oscillatory Problems with Applications to Physical Systems. Axioms 2023, 12, 282. https://doi.org/10.3390/axioms12030282

AMA Style

Kwari LJ, Sunday J, Ndam JN, Shokri A, Wang Y. On the Simulations of Second-Order Oscillatory Problems with Applications to Physical Systems. Axioms. 2023; 12(3):282. https://doi.org/10.3390/axioms12030282

Chicago/Turabian Style

Kwari, Lydia J., Joshua Sunday, Joel N. Ndam, Ali Shokri, and Yuanheng Wang. 2023. "On the Simulations of Second-Order Oscillatory Problems with Applications to Physical Systems" Axioms 12, no. 3: 282. https://doi.org/10.3390/axioms12030282

APA Style

Kwari, L. J., Sunday, J., Ndam, J. N., Shokri, A., & Wang, Y. (2023). On the Simulations of Second-Order Oscillatory Problems with Applications to Physical Systems. Axioms, 12(3), 282. https://doi.org/10.3390/axioms12030282

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop