








A Mathematica program for the approximate analytical solution to a nonlinear undamped Duffing equation by a new approximate approach






Organization:  Shanghai University 






COMPUTER PHYSICS COMMUNICATIONS 






According to Mickens [R.E. Mickens, Comments on a Generalized Galerkin's method for nonlinear oscillators, J. Sound Vib. 118 (1987) 563], the general HB (harmonic balance) method is an approximation to the convergent Fourier series representation of the periodic solution of a nonlinear oscillator and not an approximation to an expansion in terms of a small parameter. Consequently, for a nonlinear undamped Duffing equation with a driving force B cos(omega x), to find a periodic solution when the fundamental frequency is identical to omega, the corresponding Fourier series can be written as y(x) = Sigma(m)(n=1)a(n) cos[(2n  1)omega x]. How to calculate the coefficients of the Fourier series efficiently with a computer program is still an open problem. For HB method, by substituting approximation y(x) into force equation, expanding the resulting expression into a trigonometric series, then letting the coefficients of the resulting lowestorder harmonic be zero, one can obtain approximate coefficients of approximation (x) [R.E. Mickens, Comments on a Generalized Galerkin's method for nonlinear oscillators, J. Sound Vib. 118 (1987) 563]. But for nonlinear differential equations such as Duffing equation, it is very difficult to construct higherorder analytical approximations, because the HB method requires solving a set of algebraic equations for a large number of unknowns with very complex nonlinearities. To overcome the difficulty, forty years ago, Urabe derived a computational method for Duffing equation based on Galerkin procedure [M. Urabe, A. Reiter, Numerical computation of nonlinear forced oscillations by Galerkin's procedure, J. Math, Anal. Appl. 14 (1966) 107140]. Dooren obtained an approximate solution of the Duffing oscillator with a special set of parameters by using Urabe's method [R. van Dooren, Stabilization of Cowell's classic finite difference method for numerical integration, J. Comput. Phys. 16 (1974) 186192]. In this paper, in the frame of the general HB method, we present a new iteration algorithm to calculate the coefficients of the Fourier series. By using this new method, the iteration procedure starts with a (x) cos(cox) + b(x) sin(cox), and the accuracy may be improved gradually by determining new coefficients a(1), a(2).... will be produced automatically in an onebyone manner. In all the stage of calculation, we need only to solve a cubic equation. Using this new algorithm, we develop a Mathematica program, which demonstrates following main advantages over the previous HB method: (1) it avoids solving a set of associate nonlinear equations; (2) it is easier to be implemented into a computer program, and produces a highly accurate solution with analytical expression efficiently. It is interesting to find that, generally, for a given set of parameters, a nonlinear Duffing equation can have three independent oscillation modes. For some sets of the parameters, it can have two modes with complex displacement and one with real displacement. But in some cases, it can have three modes, all of them having real displacement. Therefore, we can divide the parameters into two classes, according to the solution property: there is only one mode with real displacement and there are three modes with real displacement. This program should be useful to Study the dynamically periodic behavior of a Duffing oscillator and can provide an approximate analytical solution with highaccuracy for testing the error behavior of newly developed numerical methods with a wide range of parameters.












Duffing equation, harmonic balance method, nonlinear oscillator, Fourier series













