T.A.ElgoharyL.DongJ.L.Junkins and S.N.Atluri

1 Introduction

A second-order nonlinear dynamic system can be generally recast into a system of first-order ODEs as:

which can be simply rewritten as:

In explicit methods,the future unknown state is directly expressed in terms of the currently-known system state with an explicit formula.The simplest explicit method is the forward Euler-method:

which is a first-order Taylor series expansion in the time domain.

The Runge-Kutta,RK,family of methods can be considered as the most widely used explicit methods for numerical integration of general dynamical systems.The first orderRKmethod is simply the forward Euler-method given in Eq.3.The classical or the 4th-orderRKmethod,which evaluates the solution in 4 steps,is the most commonly used among variousRKmethods.In[Fehlberg(1969)],adaptive step-size 4th-orderRKmethods are developed and are now known as the Runge-Kutta-Fehlberg,RKF,methods.Several higher order adaptiveRKFmethods,[Filippi and Gräf(1986)],are widely used for very-high accuracy applications such as orbit propagation problems,see[Montenbruck(1992);Sharp(2006)].

Implicit methods put the currently-known state and the unknown future state in a set of linear or nonlinear algebraic equations,by solving which the future state can be obtained.Backward Euler-Method is an illustration of this concept:

In[Newmark(1959)],Newmark introduced the Newmark-βmethod based on the extended mean value theorem,which is among the most widely-used implicit methods for the numerically evaluating the dynamical response of engineering structures,

What is considered as a generalization of the Newmark-βmethod is introduced in the Hilber-Hughes-Taylor orHHT-αmethod,[Hilber,Hughes,and Taylor(1977)]:

For all the above mentioned explicit and implicit methods,the size of time-steps plays an important role for the accuracy of computational results.Generally speaking,numerical stability is not guaranteed for explicit methods,thus a much smaller time step is necessary for explicit methods to obtain an accurate solution.On the other hand,because an implicit method requires the solution of set of linear/nonlinear algebraic equations in each time step,the computational burden/time of implicit methods in each time step is much higher than explicit methods.The readerisreferred to[Subbarajand Dokainish(1989);Dokainish and Subbaraj(1989)]for a comprehensive review of such methods with applications in computational structural dynamics.

Besides all of the above-mentioned widely-used numerical integrators,Eq.2 as a set of first order ODEs can be numerically solved by a wide spectrum of computational methods,see[Atluri(2005)].These methods in the time domain,such as collocation,finite volume,Galerkin,MLPG,are all essentially branches from the same tree,using the concept of weighted-residual weak-forms,and with different trial and test functions,see[Dong,Alotaibi,Mohiuddine,and Atluri(2014)].Among these methods,collocation is one of the simplest and the most efficient ones.In[Dai,Schnoor,and Atluri(2012)],a collocation method with harmonic trial functions was developed for studying the periodic responses of nonlinear Duffing oscillators and aeroelastic systems.In[Elgohary,Dong,Junkins,and Atluri(2014b)],a time-domain collocation method with Radial Basis Functions as trial function was also developed for direct solution of various time-domain inverse problems such as optimal control and orbital transfer.

In this study,the previous work of[Elgohary,Dong,Junkins,and Atluri(2014b)]is further recast as a general time-domain step-wise numerical integrator,to numerically integrate the IVP in Eq.2 for arbitrary nonlinear systems.The algorithm is compared against the central differenceExplicitmethod,the classical Runge-Kutta,RK4,method,theNewmwark-βmethod,and theHHT-αmethod in terms of accuracy,step size and computational time,using various free-vibration,forced vibration,impact load problems of single-DOF as well as coupled multi-DOF Duffing oscillators.These numerical examples clearly show the advantages of the presentRBF-Collalgorithm which enables a much larger time step,and produces high solution accuracy while maintaining a relatively low computational cost.

2 Radial Basis Function&Collocation Time Integrator

wherec>0 is a shaping parameter.

The state of the dynamical system at each LGL node can therefore be expressed in terms of the undetermined coefficients of RBF basis functions:

In matrix-vector form Eq.9 can be rewritten as,

where,Dis called the derivative matrix,which is numerically evaluated from the RBF basis functions and their time-derivatives evaluated at the LGL nodes.

The algorithm as presented above is very simple and easy to implement.The set of nonlinear algebraic equations is solved with the classical Newton’s method in this study,whereas other Jacobian-inverse-free methods can also be applied as in[Liu,Yeih,Kuo,and Atluri(2009);Liu and Atluri(2012);Elgohary,Dong,Junkins,and Atluri(2014a)].With Eq.13 being solved,the unknown states at each LGL node,as well as the unknown state at the end of this time interval or step(xxxi)are obtained.In this way,the states at each time step within the entire time history,i.e.t0,t1,t2,...,tFcan be numerically evaluated by applying the above numerical algorithm in a sequential procedure.

3 Numerical Results

In this section numerical experiments are presented for three cases of nonlinear dynamical systems,namely:the highly nonlinear unforced single Duffing oscillator,a single Duffing oscillator subjected to an impact load,and a 3-Degreeof-Freedom(3-DOF)coupled nonlinear Duffing oscillator system subjected to a harmonic load.For each case,the present RBF-Collocation,RBF-Coll,algorithm is compared against the 4th order Runge-Kutta(RK4),theExplicitmethod,[Belytschko(1976);Belytschko,Lin,and Chen-Shyh(1984)],theNewmark-βmethod,[Newmark(1959)]and theHHT-αmethod,[Hilber,Hughes,and Taylor(1977)].A table is presented to show the time step∆t,the Root-Mean-Square(RMS)error of the position state and the computation time for each case.Plots of the error time history normalized by the response amplitude,Eq.14,are also presented for each case to demonstrate the accuracy of each algorithm.Furthermore,the analytical solution and the solutions obtained by various numerical methods are plotted to highlight the magnitudes of the errors for each algorithm.

For each method,two numerical cases are presented.The first,denoted by the superscript 1,uses the same∆tas theRBF-Collalgorithm,for all the other algorithms.This establishes a baseline for comparing solution accuracy and computational time of all algorithms.The second,denoted by the superscript 2,explores the ability of each algorithm to achieve/maintain higher accuracy by decreasing∆t.As expected,and also shown,all methods will require reducing the step size to achieve a higher solution accuracy except for theRBF-Collmethod,which can maintain the solution accuracy by increasing the number of collocation points without decreasing the step size.

3.1　Free Vibrations of A Highly Nonlinear Duffing Oscillator

The unforced Duffing oscillator is described by,

For a high hardening nonlinearityη≥1,Jacobi-elliptic functions provide an analytical solution for Eq.15,[Cveti´canin(2013)],as,

where,cn is the the Jacobi-elliptic function,ωthe frequency,kthe modulus,Xthe amplitude,andθthe phase angle of the function.The reader is referred to[Cveti´canin(2013)]for the detailed expression of each parameter in Eq.16.

For the parameters values shown in Tab.1,the analytical solution is compared against MATLABODE45,Explicitmethod,RK4,Newmark-β,HHT-αand the presentRBF-Collalgorithm.ForODE451andODE452the tolerances are set to 10−6and 10−9,respectively.ForRBF-Coll1,the number of collocation pointsN=7 and shaping parameterc=0.77.ForRBF-Coll2,N=35 andc=(N−1)/4∆t.Tab.2 shows the step size,the positions state,x1(t),RMS error and the computation time for all the tested methods.Fig.1 through Fig.5 show the normalized state error time history for each method,Fig.6 shows the normalized state error time history for the presentRBF-Collalgorithm and finally Fig.7 show the comparison of computations results by various methods in the last period of the freely vibrating duffing oscillator.

Table 1:Parameters for the Freely-Vibrating Duffing System

Given the nonlinear oscillator period shown in Tab.2,the numerical integration is performed for slightly above a total of 100 cycles.Clearly,all the other presented numerical methods require a much smaller step size to achieve an acceptable accuracy,whereas theRBF-Collalgorithm can maintain high accuracy with a very large step size as demonstrated in Tab.2 and Fig.6.It should also be noted thatExplicit,Newmark-β,andHHT-αmethods perform reasonably well for a short term,but accuracy is gradually lost in the long term even if a very tiny time step is used.

Figure 1:Normalized State Error Time History,ODE451,2

Figure 2:Normalized Position Error Time History,Explicit Method 1,2

Figure 3:Normalized State Error Time History,RK41,2

Figure 4:Normalized State Error Time History,Newmark-β1,2

Figure 5:Normalized State Error Time History,HHT-α1,2

Figure 6:Normalized State Error Time History,RBF-Coll1,2

Table 2:Comparison of Numerical Methods,Freely Vibrating Highly Nonlinear Duffing Oscillator

Figure 7:Solution Comparison Last Period of Integration

3.2　Forced Duffing Oscillator with Impact Loading

An impact triangular forcing function is applied to the Duffing oscillator in Eq.15 as,

wheret1defines the time interval of the applied impact force.

For the parameters shown in Tab.3,MATLABODE45is used to obtain the numerical solution which,in the absence of an analytical solution,is treated as the reference solution for Eq.17.

Table 3:Forced Duffing System Paramaters Values

In order to capture the effect of the impact forcing function,the dynamical system in Eq.17 is solved in two intervals defined byt0≤t≤t1andt1≤t≤tF.TheODE45(reference)solution is obtained by setting the tolerances in the numerical solver to the lowest possible values,10−13.ForRBF-Coll1the number of collocation pointsN=7 and the shaping parameterc=2.5.And forRBF-Coll2,N=35 andc=(N−1)/4∆t.Tab.4 shows the comparison of step size,accuracy and computation time for various numerical methods.Fig.8 shows the normalized state error time history forExplicit1andExplicit2.Errors associated withRK41,2,Newmark-β1,2andHHT-α1,2are presented in Fig.9,Fig.10 and Fig.11,respectively.Fig.12 shows the normalized state error time history for the presentRBF-Collalgorithm and finally,Fig.13 shows the reference solution and solutions fromMethod1for the last 5 sec.of integration.

Similar to the results obtained in the previous section,the presentRBF-Collalgorithm is shown to maintain the high solution accuracy and the low computational

cost when compared to the other existing numerical methods by varying the step size and the number of collocation points.By comparing the results fromMethod1andMethod2,it can be seen that the step size is the major contributor to the computational cost associated with each method.In that senseRBF-Collis superior to all other methods as there is no need to take smaller steps in order to achieve higher solution accuracy.Thus,RBF-Collmay be a useful method to study a periodic and chaotic responses in nonlinear dynamical systems.

Table 4:Comparison of Numerical Methods,Forced Duffing Oscillator

Figure 8:Normalized Position Error Time History,Explicit Method 1,2

Figure 11:Normalized State Error Time History,HHT-α1,2

Figure 12:Normalized State Error Time History,RBF-Coll1,2

3.3　Multi Degrees of Freedom Coupled Nonlinear Dynamical System

The 3-DOF coupled nonlinear system shown in Fig.14 is analyzed in this section.WithF(t)=FcosΩt,the set of 3 coupled ODEs for this nonlinear system are,

The parameters selected for the 3-DOF system are given in Tab.5.The reference solution for this coupled nonlinear system,as in the previous section,is obtained with MATLABODE45and shwon in Fig.15.The comparison between various

Figure 14:3-DOF Nonlinear System

Figure 15:Dynamical response of the 3-DOF System

numerical integrators is given in Tab.6.ForRBF-Coll1,the number of collocation pointsN=7 and the shaping parameterc=0.35 and forRBF-Coll2,N=35 andc=0.39.Fig.16 through Fig.22 show the normalized state error time history forExplicit1,2,RK41,2,Newmark-β1,2andHHT-α1,2.TheRBF-Collnormalized state error time history is shown in Fig.24.Finally,Fig.25 shows the reference solution andMethod1solutions for the last 10 sec.of integration.

Table 5:Paramaters for the 3-DOF Coupled System

Figure 16:Normalized Position Error Time History,Explicit Method 1,2

Table 6:Comparison of Numerical Methods,3-DOF System

Figure 17:Normalized State Error Time History,RK41

Figure 18:Normalized State Error Time History,RK42

Figure 19:Normalized State Error Time History,Newmark-β1

Figure 20:Normalized State Error Time History,Newmark-β2

Figure 21:Normalized State Error Time History,HHT-α1

Figure 22:Normalized State Error Time History,HHT-α2

Figure 23:Normalized State Error Time History,RBF-Coll1

Figure 24:Normalized State Error Time History,RBF-Coll2

Figure 25:Solution Comparison,190≤t≤200

As in the previous two sections the presentRBF-Collalgorithm has the best combination of accuracy and computational cost which allows the user to take a larger time step and maintain a very high solution accuracy.In general,the computational cost ofRBF-Collis very comparable to the fastest,and as a consequence least accurate methods,which gives the algorithm a significant competitive advantage among various numerical integrators of general dynamical systems.

4 Conclusion

The presentRBF-Collalgorithm is shown to be highly accurate,fast and very simple to implement for various types of dynamical systems.Comparing the algorithm versus several implicit and explicit numerical integration methods clearly shows the advantages of such an algorithm that enables larger time step,high solution accuracy while maintaining a relatively low computational cost.For all the numerical examples shown in this study,RBF-Coll’s combination of accuracy and computational cost is superior to all other existing methods.The algorithm is shown to accurately and simply handle short and long periods of time integration,single and multi degrees of freedom system and finally transient and periodic solutions.Thus theRBF-Collalgorithm has a significant potential in handling various types of dynamical systems governed by second or higher order differential equations.Applications of the algorithm are,but not limited to,orbit propagation in celestial mechanics,dynamic buckling problems,optimal control and two-point boundary value problems.Areas of these studies will be explored in future works.

Acknowledgement:This work is supported by the US Army Research Labs VT Division,under an AFRL/UCI collaborative research agreement.The encouragement of Messrs.Dy Le and Jarret Riddick is thankfully acknowledged.This work is also supported by the Texas A&M Institute for Advanced Study(TIAS).It was initiated while S.N.Atluri visited TIAS briefly in January,2014.

Atluri,S.N.(2005):Methods of computer modeling in engineering&the sciences,volume 1.Tech Science Press Palmdale.

Belytschko,T.(1976): A survey of numerical methods and computer programs for dynamic structural analysis.Nuclear Engineering and Design,vol.37,no.1,pp.23–34.

Belytschko,T.;Lin,J.I.;Chen-Shyh,T.(1984): Explicit algorithms for the nonlinear dynamics of shells.Computer methods in applied mechanics and engineering,vol.42,no.2,pp.225–251.

Buhmann,M.D.(2003):Radial basis functions:theory and implementations,volume 5.Cambridge university press Cambridge.

Cveti´canin,L.(2013):Ninety years of duffing’s equation.Theoretical and Applied Mechanics,vol.40,no.1,pp.49–63.

Dai,H.-H.;Schnoor,M.;Atluri,S.N.(2012):A simple collocation scheme for obtaining the periodic solutions of the duffing equation,and its equivalence to the high dimensional harmonic balance method:subharmonic oscillations.Computer Modeling in Engineering and Sciences,vol.84,no.5,pp.459–497.

Dokainish,M.;Subbaraj,K.(1989):Asurvey ofdirecttime-integration methods in computational structural dynamics-i.explicit methods.Computers&Structures,vol.32,no.6,pp.1371–1386.

Dong,L.;Alotaibi,A.;Mohiuddine,S.;Atluri,S.(2014): Computational methods in engineering:a variety of primal&mixed methods,with global&local interpolations,for well-posed or ill-posed bcs.CMES:Computer Modeling in Engineering&Sciences,vol.99,no.1,pp.1–85.

Elgohary,T.A.;Dong,L.;Junkins,J.L.;Atluri,S.N.(2014): Solution of post-buckling&limit load problems,without inverting the tangent stiffness matrix&without using arc-length methods.CMES:Computer Modeling in Engineering&Sciences,vol.98,no.6,pp.543–563.

Elgohary,T.A.;Dong,L.;Junkins,J.L.;Atluri,S.N.(2014):Time domain inverse problems in nonlinear systems using collocation&radial basis functions.CMES:Computer Modeling in Engineering&Sciences,vol.100,no.1,pp.59–84.

Fehlberg,E.(1969): Low-order classical runge-kutta formulas with stepsize control and their application to some heat transfer problems.Technical report,NASA,1969.

Filippi,S.;Gräf,J.(1986): New runge–kutta–nyström formula-pairs of order 8(7),9(8),10(9)and 11(10)for differential equations of the formy”=f(x,y).Journal of computational and applied mathematics,vol.14,no.3,pp.361–370.

Hilber,H.M.;Hughes,T.J.;Taylor,R.L.(1977):Improved numerical dissipation for time integration algorithms in structural dynamics.Earthquake Engineering&Structural Dynamics,vol.5,no.3,pp.283–292.

Liu,C.-S.;Atluri,S.N.(2012):A globally optimal iterative algorithm using the best descent vector˙x=λ?αcF+BTF?,with the critical valueαc,for solving a system of nonlinear algebraic equationsF(x)=0.Computer Modeling in Engineering and Sciences,vol.84,no.6,pp.575.

Liu,C.S.;Yeih,W.;Kuo,C.L.;Atluri,S.N.(2009): A scalar homotopy method for solving an over/under determined system of non-linear algebraic equations.Computer Modeling in Engineering and Sciences,vol.53,no.1,pp.47–71.

Montenbruck,O.(1992): Numerical integration methods for orbital motion.Celestial Mechanics and Dynamical Astronomy,vol.53,no.1,pp.59–69.

Newmark,N.M.(1959): A method of computation for structural dynamics.Journal of the Engineering Mechanics Division,vol.85,no.3,pp.67–94.

Noor,A.K.;Lambiotte Jr,J.J.(1979):Finite element dynamic analysis on cdc star-100 computer.Computers&Structures,vol.10,no.1,pp.7–19.

Sharp,P.W.(2006):N-body simulations:The performance of some integrators.ACM Transactions on Mathematical Software(TOMS),vol.32,no.3,pp.375–395.

Subbaraj,K.;Dokainish,M.(1989):A survey of direct time-integration methodsin computationalstructuraldynamics-ii.implicitmethods.Computers&Structures,vol.32,no.6,pp.1387–1401.