# The Galerkin Method for Solving Strongly Nonlinear Oscillators

Authors: Alvaro H. S. Salas
Journal: The Scientific World Journal
 (2022)

---

## Abstract

In this paper, we make use of the Galerkin method for solving nonlinear second-order ODEs that are related to some strongly nonlinear oscillators arising in physics and engineering. We derive the iterative schemes for finding the coefficients that appear in the linear Galerkin hat combination in the ansatz form solution. These coefficients may be found iteratively by solving either a quadratic or a higher degree algebraic equation. Examples are presented to illustrate the obtained results. Some exact solutions are given, and they are compared with both the Runge–Kutta numerical solution and the solution obtained using the Galerkin finite element method.


---

## Body

## 1. Introduction

The nonlinear equation describing an oscillator with a cubic nonlinearity is called the Duffing equation. This equation has a variety of applications in science and engineering, early mechanical failure signal, nonlinear circuit design [1], image processing [2], vibration of buckled beams [3], solitons [4–6], chaos [7], and many areas of physics.There are many methods for solving nonlinear differential equations. In this paper, we concentrate on the numerical solution to the Duffing equation by means of the finite element method. This method is due to Galerkin, a Russian engineer and scientist. We also derive formulas for solving a wide class of nonlinear oscillators. The Galerkin solutions are compared with the solutions obtained using the Runge–Kutta numerical method.


## 2. Finite Element Method or Galerkin Hat Method

Let us consider a polynomial second-order damped and forced ode:(1)x¨+2εx˙+Px=ft,x0=0,x′0=x˙0,where x≡xt and P≡Px is a polynomial whose coefficients depend on t, say(2)Px=∑s=0Mastxs,x=xt.Given the i.v.p.(3)y¨+2εy˙+Qy=ft,y0=y0≠0,andy′0=y˙0,whereQy=∑r=0Mbrtyr,y=yt.we define(4)yt=y0+xt,x0=0,x′0=y0.Then,(5)Qy=∑r=0Mbry0+xtr=∑r=0Mbr∑s=0rrsy0r−sxst=∑s=0M∑r=sMrsy0r−sbrxst,=∑s=0Masxs≕Px,x=xt,whereas=∑r=sMrsy0r−sbr.Thus, the problem reduces to (1). Some particular cases to i.v.p. (1) are(6)x¨+2εx˙+nt+ptx=ft,x¨+2εx˙+nt+ptx+qtx2=ft,x¨+2εx˙+nt+ptx+qtx2+rtx3=ft.We will use the same idea as for the linear case(7)x¨+2εx˙+ptx=ft,x0=0,x′0=x˙0,that is, we will assume an approximate analytical solution in the ansatz form(8)x=xt=∑k=1nckφkt,where the functions φkt are the so-called linear Galerkin hats. Let 0≤t≤T. Choose some positive integer n≥2 and define the step h=T/n and let ξj=jh=jT/n for j=0,1,2,… The functions φkt are defined as follows.(9)φkt=t−ξk−1hifξk−1≤t≤ξk,−t−ξk+1hifξk≤t≤ξk+1,fork=1,2,…,n.For an illustration, see Figure1.Figure 1
Galerkin hats forn=7 and T=5.Some properties of these functions can be illustrated as follows.(10)φjtφktdt=0,for j−k≥2 and t∈0,T, and(11)∫0Tφjpt=2Tp+1n,for j≥1 and p=1,2,3,….In general,(12)∫0Tφjrtφkst=Tr!s!nr+s+1!,forj−k=1,r,s=0,1,2,3,….Using the formula(13)∫0TφjtxtNdt=hN+1N+2∑k=0N−1k+1cj−1N−k+cj+1N−kcjk+2N+1cjNfor anyN≥0,c0=cn+1=0,and assuming that ajt≡aj= const, we may evaluate easily the following integral:(14)∫0TPxφjtdt=∑N=0maN∫0TφjtxtNdtfor anyj.Moreover,(15)∫0Tx′tφjtdt=−∫0Tx′tφj′tdt=cj−1−2cj+cj+1hfor anyj,∫0Tx′tφjtdt=12cj+1−cj−1.Thus, for example, ifa0, a1, …are independent of t,(16)∫0Tx″t+δx′t+a0+a1xt+a2x2t+a3x3t+a4x4t+a5x5tφjtdt=ncj−1−2cj+cj+1T+δ2cj+1−cj−1+Tna0+a1A1+a2A2+a3A3+a4A4+a5A5,where(17)A1=T6ncj−1+4cj+cj+1,A2=T12ncj−12+2cjcj−1+6cj2+cj+12+2cjcj+1,A3=T20ncj−13+2cjcj−12+3cj2cj−1+8cj3+cj+13+2cjcj+12+3cj2cj+1,A4=T30ncj−14+2cjcj−13+3cj2cj−12+4cj3cj−1+10cj4+cj+14+2cjcj+13+3cj2cj+12+4cj3cj+1,A5=T42ncj−15+2cjcj−14+3cj2cj−13+4cj3cj−12+5cj4cj−1+12cj5+cj+15+2cjcj+14+3cj2cj+13+4cj3cj+12+5cj4cj+1.In general, ifar does not depend on t,(18)∫0Tx′′t+δx′t+∑s=0Masxstφjtdt=ncj−1−2cj+cj+1T+12δcj+1−cj−1+Ta0n+∑s=0MAsas,where(19)As=Tns+1s+2∑k=0s−1k+1cj−1s−k+cj+1s−kcjk+2s+1cjs.Let us consider the forced and damped oscillator (1). Following are useful expressions for different forces:(20)∫0Tftφjtdt=4nTΩ2sin2TΩ2nF0Tn+F1cosjTΩn+F2sinjTΩn,ft=f0+f1cosΩt+f2sinΩt,Ω≠0,∫0Tftφjtdt=ne−j+1λT/nTλ2+Ω22−2λΩsinj+1TΩn+λ−Ωλ+ΩeλT/ncosj−1TΩn+λ−Ωλ+Ωcosj+1TΩn−2eλT/nλΩeλT/nsinj−1TΩn−2sinjTΩn+λ−Ωλ+ΩcosjTΩn,ft=e−λtcosΩt,Ω2+λ2≠0,∫0Tftφjtdt=ne−j+1λT/nTλ2+Ω22λ−Ωλ+ΩeλT/neλT/nsinj−1TΩn−2sinjTΩn+λ−Ωλ+Ωsinj+1TΩn+2λΩe2λT/ncosj−1TΩn−4λΩeλT/ncosjTΩn+2λΩcosj+1TΩn,ft=e−λtsinΩt,Ω2+λ2≠0,∫0Tftφjtdt=λ−ρ−2TλT1−jΓρ+1,j−1Tλn+2jΓρ+1,jTλn−j+1Γρ+1,j+1Tλn+nΓρ+2,j−1Tλn−2Γρ+2,jTλn+Γρ+2,j+1Tλn,ft=e−λttρ,∫0Tftφjtdt≈n9TwA−1cos3j−1Twn−2cos3jTwn+cos3j+1Twn−18AcosTwn−1cosjTwn,ft=cnΩ,m≈Acoswt+1−Acos3wt.w=27m2−92m+64Ω23m2−60m+64.A=215m2−48m+3227m2−92m+64,m≤0.5,∫0Tftφjtdt≈n9TwA−1cos3j−1Twn−2cos3jTwn+cos3j+1Twn−18AcosTwn−1cosjTwn,ft=cnΩt,m≈Acoswt+1−Acos3wt.w=27m2−92m+64Ω23m2−60m+64.A=215m2−48m+3227m2−92m+64,m≤0.5,∫0Tftφjtdt≈3−4A2n18AcosTΩ/3n−4An−1sinjTΩ/3n−4An−A−1sin3j−1TΩ/4A−3n+sin3j+1TΩ/4A−3n+2sin3jTΩ/3n−4An9TΩ2,ft=snΩt,m≈AsintΩ4A−3−1−Asin3tΩ4A−3,A=−−25m2+144m−12829m2−152m+128,m≤0.25.Other useful formulas when theaj depend on time are(21)∫0Ttxtφjtdt=112hcj−12a+h2j−1+8cja+hj+cj+12a+2hj+h,∫0Tt2xtφjtdt=160h340j2+4cj+10j−1j+3cj−1+10jj+1+3cj+1,∫0Tt3xtφjtdt=160h44j10j2+3cj+2j−15j−1j+2cj−1+2j+15jj+1+2cj+1,∫0Ttxt2φjtdt=160h25j−3cj−12+25j−2cjcj−1+30jcj2+5j+3cj+12+25j+2cjcj+1,∫0Tt2xt2φjtdt=160h330j2+2cj2+j5j−6+2cj−12+2j5j−4+1cjcj−1+j5j+6+2cj+12+2j5j+4+1cjcj+1,∫0Tt3xt2φjtdt=1420h4425j3+jcj2+7jj5j−9+6−10cj−12+27jj5j−6+3−4cjcj−1+7jj5j+9+6+10cj+12+27jj5j+6+3+4cjcj+1,∫0Ttxt3φjtdt=160h23j−2cj−13+32j−1cjcj−12+33j−1cj2cj−1+24jcj3+3j+2cj+13+32j+1cjcj+12+33j+1j2cj+1,∫0Tt2xt3φjtdt=1420h3821j2+1cj3+7j3j−4+10cj−13+67j−1j+2cjcj−12+37j3j−2+3cj2cj−1+7j3j+4+10cj+13+67jj+1+2cjcj+12+37j3j+2+3cj2cj+1,∫0Tt3xt3φjtdt=1280h4167j3+jcj3+2j7j−2j+10−5j−13+2j−114j−1j+5cjcj−12+32j7j−1j+3−1cj2cj−1+2j7jj+2+10+5cj+13+2j+114jj+1+5cjcj+12+32j7jj+1+3+1cj2cj+1.Other formulas for calculating are given in the Appendix.


## 3. Applications

### 3.1. Linear Oscillator

This is the ode(22)x¨+α+βx=0,β≠0,x0=0,x′0=x˙0,0≤t≤T.The exact solution is given by(23)xt=1βαcosβt−1+βx˙0sinβt.Assume the ansatz(24)x=xt=∑k=1nckφkt.Definec−1=c0=0,c1=T/nx˙0 and(25)Rj=∫0Tx¨+α+βxφjtdt,j=0,1,2,…,n−1.Then,(26)Rj=6αT26nT+6n2+βT26nTcj−1+4cβT2−3n26nTcj+βT6n+nTcj+1=0for anyj.This is a linear recurrence that may be solved in closed form:(27)cj=16β−6α+3αT3ββT2−12n2−2βT+6n26n2+βT2j+3α6n2−T3ββT2−12n2+2βT6n2+βT2j+1nT2/3−4n2/βT3ββT2−12n2−2βT+6n2j−6n2−T3ββT2−12n2+2βTj6n2+βT2−jx˙06n2+βT2+3αnT,j=1,2,3,…,n.Example 1.
Letα=0, β=1, x˙0=1, and T=2π. The exact solution to (14) is xt=sint. From (19), we obtain(28)cj=17440418π2/n2+26160627295378481881−26160627π2/n2sinjtan−15906π78481881−26160627π2/n2n26160627−34880836π2/n2.
In Figure2, we compare the exact solution with the approximate solution (pairs 2π/nj,cj,j=1,2,3,…,n) for n=30.Figure 2
Comparison between the Galerkin solution and the Runge–Kutta numerical solution.


### 3.2. Undamped and Unforced Helmholtz Oscillator

This is the ode(29)y¨+α+βy+γy2=0,y0=y0,y′0=y˙0,0≤t≤T.The exact solution to (22) may be expressed in any of the following equivalent forms:(30)yt=−β2γ−6℘t+t0;1/12β2−4αγ,1/216β+2γy0β2−6αγ−2γy0β+γy0−6γ2y˙02γ,wheret0=±℘−11/12−β−2γy0;1/12β2−4αγ,1/216β+2γy0−6αγ+β2−2γy0β+γy0−6γ2y˙02,yt=A−6α+AAγ+β12℘t+t0;1/12β2−4αγ,1/216β+2Aγβ2−2Aγβ−2γγA2+3α+2Aγ+β,where−2γA3+3βA2+6αA+6αy0−3βy02−2γy03−3y˙02=0,t0=℘−16α+4A2γ+5Aβ+2Aγy0+βy012A−y0;112β2−4αγ,12162Aγ+β−2γ3α+A2γ−2Aβγ+β2,yt=A+y0−Acnωt,m+b1snωt,mdnωt,m1+b2sn2ωt,m2,where(31)b1=y˙02ωy0−A,b2=y˙02−2A−y0α+2Aω+y0β+γy0−2ω8ωA−y02,ω=−−2A2β−4A2γy0−6αA−2Aβy0+2Aγy02+6αy0+4βy02+2γy03+3y˙0242m−1A−y02,m=121−32Aγ+β22Aγ+β2β−2Aγ2Aγ+3β−16αγ,2γA3+3βA2+6αA−6αy0−3βy02−2γy03−3y˙02=0.The period equalsT=2Km/wIn order to obtain approximate solution by means of the Galerkin method, let(32)yt=y0+xt,x0=0,x′0=y˙0.Then,(33)x¨+α+βy0+γy02+β+2γy0xt+γxt2=0,x0=0,x′0=y˙0.Assume the ansatz(34)x=xt=∑k=1nckφkt.Definec−1=c0=0,c1=T/nx˙0 and(35)Rj=∫0Tx¨+α+βy0+γy02+β+2γy0xt+γxt2φjtdt.Then,(36)Rj=8cj−3n2+βT2+2γT2y0+2cj−1γT2cj+2y0+6n2+βT2+γT2cj−12+6γT2cj2+12T2α+y0β+γy012nT+γT2cj+2y0+6n2+βT26nTcj+1+γT12ncj+12=0.Observe that the equationRj=0 is a quadratic equation in z=cj+1. It is clear that the system R0=R1=…=Rn−1=0 may be solved recursively. We first find c2 letting j=1 in (29):(37)c2=−6n2−T2β+2γy0±6n2+βT22+γT2−2nTy˙06n2+βT2+2γT2y0−8n2y0−3n2+βT2+γT2y0−γT4y˙02/n2−12αγT4γT2.We choose the value ofc2 that is closest to c1=T/nx˙0. Next, we set j=2 in (37) and then we will find z=c3. We choose the closest to c2 solution to the quadratic equation in (29). We continue this procedure and then we will find all values of cjj=1,2,3,…n. Since xξj=xT/nj=cj, all pairs ξj,cj will lie on the graph of the solution for sufficiently large n. Plotting these points, we obtain the graph of the solution.On the other hand, ifn is large enough, the values of cj and cj+1 will be close to each other. We may use the following approximate expression for cj+1 in terms of cj−1 and cj:(38)cj+1≈−8cj−3n2+βT2+2γT2y0+2cj−1γT2cj+2y0+6n2+βT2+γT2cj−12+5γT2cj2+12T2α+y0β+γy022γT2cj+y0+6n2+βT2.Example 2.
Let(39)y¨+2yt−y2t=0∧y0=0∧y′0=1.
The exact solution is obtained from (25) and (26):(40)yt=−0.641784+0.801114cn0.574208t|−0.744849+1.08694dn0.574208t|−0.744849sn0.574208t|−0.74484921+0.420436sn0.574208t|−0.74484922.
The solution is periodic and its period equalsT=4.737423705838371 (see Figure 3).Figure 3
Comparison between the exact solution and the finite element method forn=30.


### 3.3. Duffing–Helmholtz Oscillator

Let(41)q¨+α+βq+γq2+δq3=0,q0=q0,q′0=q˙0.The exact solution to i.v.p. (34) is given by(42)qt=A+B1+C℘t+t0;g2,g3,where(43)B=−6A3δ+A2γ+Aβ+α3A2δ+2Aγ+β,C=123A2δ+2Aγ+β,t0=±℘−1q0−A−BCA−q0;g2,g3,g2=112−3A4δ2−4A3γδ−6A2βδ−12Aαδ−4αγ+β2,g3=1216−3A4δγ2−3βδ−4A3γγ2−3βδ+6A2β3βδ−γ2−12Aαγ2−3βδ+27α2δ−6αβγ+β3,3δA4+4γA3+6βA2+12Aα−12q0α+6q02β+4q03γ+3q04δ+6q˙02=0.This solution is valid even ifα=γ=0. The solution is periodic and its period equals(44)T=2∫a∞14x3−g2x−g3,where a is the greatest real root to the cubic 4x3−g2x−g3=0. Assume that 4x3−g2x−g3=4x−ax−bx−c. Then, the period may be evaluated using the formulas(45)∫audxx−ax−bx−c=2c−bFsin−1b−cu−ab−au−c,b−ab−c,∫a∞dxx−ax−bx−c=limu⟶∞2c−bFsin−1b−cu−ab−au−c,b−ab−c=2c−bKb−ab−c,∫audxx−ax−bx−b¯=Fa−α2+β24cos−1a−u+a−α2+β2−a+u+a−α2+β2,α−a+a−α2+β22a−α2+β2,∫a∞dxx−ax−bx−b¯=limu⟶∞∫audxx−ax−bx−b¯=2a−α2+β24K−a+α+a−α2+β22a−α2+β2,α=Reb,β=Imb.Another expression for the exact solution is given by(46)qt=A+B1+v0cna+bv02t,bv02/2a+bv02,where(47)v0=q0−A−BA−q0,3δA4+4γA3+6βA2+12αA−12q0α+6q02β+4q03γ+3q04δ=0,a=A4−γ2−3βδ−2A3βγ−9αδ−3A212q0αδ+6q02βδ+4q03γδ+3q04δ2−2αγ+β2−2A12q0αγ+6q02βγ+4q03γ2+3q04γδ+3αβ−12q0αβ−6q02β2−4q03βγ−3q04βδ−9α2A−q03A3δ+A23q0δ+4γ+A4q0γ+3q02δ+6β+6q0β+4q02γ+3q03δ+12α,b=−3A3δ+A2γ+Aβ+α2A−q03A3δ+A23q0δ+4γ+A4q0γ+3q02δ+6β+6q0β+4q02γ+3q03δ+12α,B=−A−q03A3δ+A23q0δ+4γ+A4q0γ+3q02δ+6β+6q0β+4q02γ+3q03δ+12α3A3δ+A2γ+Aβ+α.Letq0=0. Assume the ansatz(48)q=qt=∑k=1nckφkt.Definec0=cn+1=0. We have(49)∫0Tq¨t+α+βqt+γq2t+δq3tφjtdt=cj−1−2cj+cj+1h+αh+16βhcj−1+4cj+cj+1+112γhcj−12+2cjcj−1+6cj2+cj+12+2cjcj+1+120δhcj−13+2cjcj−12+3cj2cj−1+8cj3+cj+13+2cjcj+12+3cj2cj+1forj=1,2,3,….n.This gives an algebraic system of nonlinear equations. The initial data are(50)c0=0,c1=c1=Tq˙0n.Observe that algebraic system (39) may be solved recursively using the Tartaglia formula for the cubic. Indeed, we may write(51)120δhz3+160h6δcj+5γz2+h210β+10γcj+9δcj2+6060hz+160hcj−1h2cj9δcj+10γ+10βh2+6+40βh2cj+h2cj−126δcj+5γ+6h2cj24δcj+5γ+3δh2cj−13−120cj+60αh2=0,wherez=cj+1.Example 3.
Letα=1, β=−1, γ=1, δ=1, q0=0, q˙0=0.3, and T=20 (see Figure 4 for a comparison between the exact solution and the Galerkin method for different values of n).Figure 4
Galerkin hat method for different values of n.


### 3.4. Duffing Equation with Damping and Forcing

Let(52)x¨+2εx˙+αx+βx3=Fcosωt,x0=0,x′0=x˙0,0≤t≤T.Assume the ansatz(53)x=xt=∑k=1nckφkt.Definec0=cn+1=0. We have(54)∫0Tx¨+2εx˙+αx+βx3−Fcosωtφjtdt=120βhcj+13+110βhcjcj+12+320βhcj2+ε+αh6+1hcj+1+160h85cjαh2−3+3βh2cj3−30Fω2sin2hω2coshjω+cj−19βh2cj2+10αh2−6εh+6+3βh2cj−13+6βh2cjcj−12=0,j=1,2,3,….We may solve this system recursively using the Tartaglia formula for the cubic. The initial data are(55)c0=0,c1=Tnx˙0.


### 3.5. Forced Van der Pol–Duffing Equation

Let(56)x¨−ε1−x2x˙+αx+βx3=Fcosωt,x0=0,x′0=x˙0,0≤t≤T.Assume the ansatz(57)x=xt=∑k=1nckφkt.Definec0=cn+1=0. We have(58)∫0Tx¨−ε1−x2x˙+αx+βx3−Fcosωtφjtdt=1605ε+3βhcj+13+160cj5ε+6βhcj+12+160cj25ε+9βh−ε2+αh6+1hcj+1+160h85cjαh2−3+3βh2cj3−30Fω2sin2hω2coshjω+cj−1hcj29βh−5ε+10αh2+3εh+6+hcj−133βh−5ε+hcjcj−126βh−5ε=0,forj=1,2,3,….We may solve this system recursively using the Tartaglia formula for the cubic. The initial data are(59)c0=0,c1=nx˙0T.


### 3.6. Conservative Nonlinear Oscillators

Conservative single-degree-of-freedom nonlinear oscillators are modelled by second-order autonomous ordinary differential equations of the form(60)x¨+Fx=0,x0=x0,x′0=0.Using Chebyshev polynomials or minimization techniques, we may approximate the functionFx by means of a cubic polynomial, say Fx=α+βx+γx2+δx3, and this allows us to study i.v.p. (44) using the solution to the i.v.p.(61)x¨+α+βx+γx2+δx3=0,x0=x0,x′0=0.For example, let us consider a mass attached to two stretched elastic springs [8]. For this problem, the function Fx has the form(62)Fx=x−λx1+x2,0<λ≤1.Suppose thatx≤A. We have Fx≈px+qx3, where(63)p=1+15λ20A2+21sinh−1A−3AA2+12A2+732A5,q=35λAA2+12A2+15−34A2+5sinh−1A32A7.Then, the problem reduces to that of solving a Duffing equation.


## 3.1. Linear Oscillator

This is the ode(22)x¨+α+βx=0,β≠0,x0=0,x′0=x˙0,0≤t≤T.The exact solution is given by(23)xt=1βαcosβt−1+βx˙0sinβt.Assume the ansatz(24)x=xt=∑k=1nckφkt.Definec−1=c0=0,c1=T/nx˙0 and(25)Rj=∫0Tx¨+α+βxφjtdt,j=0,1,2,…,n−1.Then,(26)Rj=6αT26nT+6n2+βT26nTcj−1+4cβT2−3n26nTcj+βT6n+nTcj+1=0for anyj.This is a linear recurrence that may be solved in closed form:(27)cj=16β−6α+3αT3ββT2−12n2−2βT+6n26n2+βT2j+3α6n2−T3ββT2−12n2+2βT6n2+βT2j+1nT2/3−4n2/βT3ββT2−12n2−2βT+6n2j−6n2−T3ββT2−12n2+2βTj6n2+βT2−jx˙06n2+βT2+3αnT,j=1,2,3,…,n.Example 1.
Letα=0, β=1, x˙0=1, and T=2π. The exact solution to (14) is xt=sint. From (19), we obtain(28)cj=17440418π2/n2+26160627295378481881−26160627π2/n2sinjtan−15906π78481881−26160627π2/n2n26160627−34880836π2/n2.
In Figure2, we compare the exact solution with the approximate solution (pairs 2π/nj,cj,j=1,2,3,…,n) for n=30.Figure 2
Comparison between the Galerkin solution and the Runge–Kutta numerical solution.


## 3.2. Undamped and Unforced Helmholtz Oscillator

This is the ode(29)y¨+α+βy+γy2=0,y0=y0,y′0=y˙0,0≤t≤T.The exact solution to (22) may be expressed in any of the following equivalent forms:(30)yt=−β2γ−6℘t+t0;1/12β2−4αγ,1/216β+2γy0β2−6αγ−2γy0β+γy0−6γ2y˙02γ,wheret0=±℘−11/12−β−2γy0;1/12β2−4αγ,1/216β+2γy0−6αγ+β2−2γy0β+γy0−6γ2y˙02,yt=A−6α+AAγ+β12℘t+t0;1/12β2−4αγ,1/216β+2Aγβ2−2Aγβ−2γγA2+3α+2Aγ+β,where−2γA3+3βA2+6αA+6αy0−3βy02−2γy03−3y˙02=0,t0=℘−16α+4A2γ+5Aβ+2Aγy0+βy012A−y0;112β2−4αγ,12162Aγ+β−2γ3α+A2γ−2Aβγ+β2,yt=A+y0−Acnωt,m+b1snωt,mdnωt,m1+b2sn2ωt,m2,where(31)b1=y˙02ωy0−A,b2=y˙02−2A−y0α+2Aω+y0β+γy0−2ω8ωA−y02,ω=−−2A2β−4A2γy0−6αA−2Aβy0+2Aγy02+6αy0+4βy02+2γy03+3y˙0242m−1A−y02,m=121−32Aγ+β22Aγ+β2β−2Aγ2Aγ+3β−16αγ,2γA3+3βA2+6αA−6αy0−3βy02−2γy03−3y˙02=0.The period equalsT=2Km/wIn order to obtain approximate solution by means of the Galerkin method, let(32)yt=y0+xt,x0=0,x′0=y˙0.Then,(33)x¨+α+βy0+γy02+β+2γy0xt+γxt2=0,x0=0,x′0=y˙0.Assume the ansatz(34)x=xt=∑k=1nckφkt.Definec−1=c0=0,c1=T/nx˙0 and(35)Rj=∫0Tx¨+α+βy0+γy02+β+2γy0xt+γxt2φjtdt.Then,(36)Rj=8cj−3n2+βT2+2γT2y0+2cj−1γT2cj+2y0+6n2+βT2+γT2cj−12+6γT2cj2+12T2α+y0β+γy012nT+γT2cj+2y0+6n2+βT26nTcj+1+γT12ncj+12=0.Observe that the equationRj=0 is a quadratic equation in z=cj+1. It is clear that the system R0=R1=…=Rn−1=0 may be solved recursively. We first find c2 letting j=1 in (29):(37)c2=−6n2−T2β+2γy0±6n2+βT22+γT2−2nTy˙06n2+βT2+2γT2y0−8n2y0−3n2+βT2+γT2y0−γT4y˙02/n2−12αγT4γT2.We choose the value ofc2 that is closest to c1=T/nx˙0. Next, we set j=2 in (37) and then we will find z=c3. We choose the closest to c2 solution to the quadratic equation in (29). We continue this procedure and then we will find all values of cjj=1,2,3,…n. Since xξj=xT/nj=cj, all pairs ξj,cj will lie on the graph of the solution for sufficiently large n. Plotting these points, we obtain the graph of the solution.On the other hand, ifn is large enough, the values of cj and cj+1 will be close to each other. We may use the following approximate expression for cj+1 in terms of cj−1 and cj:(38)cj+1≈−8cj−3n2+βT2+2γT2y0+2cj−1γT2cj+2y0+6n2+βT2+γT2cj−12+5γT2cj2+12T2α+y0β+γy022γT2cj+y0+6n2+βT2.Example 2.
Let(39)y¨+2yt−y2t=0∧y0=0∧y′0=1.
The exact solution is obtained from (25) and (26):(40)yt=−0.641784+0.801114cn0.574208t|−0.744849+1.08694dn0.574208t|−0.744849sn0.574208t|−0.74484921+0.420436sn0.574208t|−0.74484922.
The solution is periodic and its period equalsT=4.737423705838371 (see Figure 3).Figure 3
Comparison between the exact solution and the finite element method forn=30.


## 3.3. Duffing–Helmholtz Oscillator

Let(41)q¨+α+βq+γq2+δq3=0,q0=q0,q′0=q˙0.The exact solution to i.v.p. (34) is given by(42)qt=A+B1+C℘t+t0;g2,g3,where(43)B=−6A3δ+A2γ+Aβ+α3A2δ+2Aγ+β,C=123A2δ+2Aγ+β,t0=±℘−1q0−A−BCA−q0;g2,g3,g2=112−3A4δ2−4A3γδ−6A2βδ−12Aαδ−4αγ+β2,g3=1216−3A4δγ2−3βδ−4A3γγ2−3βδ+6A2β3βδ−γ2−12Aαγ2−3βδ+27α2δ−6αβγ+β3,3δA4+4γA3+6βA2+12Aα−12q0α+6q02β+4q03γ+3q04δ+6q˙02=0.This solution is valid even ifα=γ=0. The solution is periodic and its period equals(44)T=2∫a∞14x3−g2x−g3,where a is the greatest real root to the cubic 4x3−g2x−g3=0. Assume that 4x3−g2x−g3=4x−ax−bx−c. Then, the period may be evaluated using the formulas(45)∫audxx−ax−bx−c=2c−bFsin−1b−cu−ab−au−c,b−ab−c,∫a∞dxx−ax−bx−c=limu⟶∞2c−bFsin−1b−cu−ab−au−c,b−ab−c=2c−bKb−ab−c,∫audxx−ax−bx−b¯=Fa−α2+β24cos−1a−u+a−α2+β2−a+u+a−α2+β2,α−a+a−α2+β22a−α2+β2,∫a∞dxx−ax−bx−b¯=limu⟶∞∫audxx−ax−bx−b¯=2a−α2+β24K−a+α+a−α2+β22a−α2+β2,α=Reb,β=Imb.Another expression for the exact solution is given by(46)qt=A+B1+v0cna+bv02t,bv02/2a+bv02,where(47)v0=q0−A−BA−q0,3δA4+4γA3+6βA2+12αA−12q0α+6q02β+4q03γ+3q04δ=0,a=A4−γ2−3βδ−2A3βγ−9αδ−3A212q0αδ+6q02βδ+4q03γδ+3q04δ2−2αγ+β2−2A12q0αγ+6q02βγ+4q03γ2+3q04γδ+3αβ−12q0αβ−6q02β2−4q03βγ−3q04βδ−9α2A−q03A3δ+A23q0δ+4γ+A4q0γ+3q02δ+6β+6q0β+4q02γ+3q03δ+12α,b=−3A3δ+A2γ+Aβ+α2A−q03A3δ+A23q0δ+4γ+A4q0γ+3q02δ+6β+6q0β+4q02γ+3q03δ+12α,B=−A−q03A3δ+A23q0δ+4γ+A4q0γ+3q02δ+6β+6q0β+4q02γ+3q03δ+12α3A3δ+A2γ+Aβ+α.Letq0=0. Assume the ansatz(48)q=qt=∑k=1nckφkt.Definec0=cn+1=0. We have(49)∫0Tq¨t+α+βqt+γq2t+δq3tφjtdt=cj−1−2cj+cj+1h+αh+16βhcj−1+4cj+cj+1+112γhcj−12+2cjcj−1+6cj2+cj+12+2cjcj+1+120δhcj−13+2cjcj−12+3cj2cj−1+8cj3+cj+13+2cjcj+12+3cj2cj+1forj=1,2,3,….n.This gives an algebraic system of nonlinear equations. The initial data are(50)c0=0,c1=c1=Tq˙0n.Observe that algebraic system (39) may be solved recursively using the Tartaglia formula for the cubic. Indeed, we may write(51)120δhz3+160h6δcj+5γz2+h210β+10γcj+9δcj2+6060hz+160hcj−1h2cj9δcj+10γ+10βh2+6+40βh2cj+h2cj−126δcj+5γ+6h2cj24δcj+5γ+3δh2cj−13−120cj+60αh2=0,wherez=cj+1.Example 3.
Letα=1, β=−1, γ=1, δ=1, q0=0, q˙0=0.3, and T=20 (see Figure 4 for a comparison between the exact solution and the Galerkin method for different values of n).Figure 4
Galerkin hat method for different values of n.


## 3.4. Duffing Equation with Damping and Forcing

Let(52)x¨+2εx˙+αx+βx3=Fcosωt,x0=0,x′0=x˙0,0≤t≤T.Assume the ansatz(53)x=xt=∑k=1nckφkt.Definec0=cn+1=0. We have(54)∫0Tx¨+2εx˙+αx+βx3−Fcosωtφjtdt=120βhcj+13+110βhcjcj+12+320βhcj2+ε+αh6+1hcj+1+160h85cjαh2−3+3βh2cj3−30Fω2sin2hω2coshjω+cj−19βh2cj2+10αh2−6εh+6+3βh2cj−13+6βh2cjcj−12=0,j=1,2,3,….We may solve this system recursively using the Tartaglia formula for the cubic. The initial data are(55)c0=0,c1=Tnx˙0.


## 3.5. Forced Van der Pol–Duffing Equation

Let(56)x¨−ε1−x2x˙+αx+βx3=Fcosωt,x0=0,x′0=x˙0,0≤t≤T.Assume the ansatz(57)x=xt=∑k=1nckφkt.Definec0=cn+1=0. We have(58)∫0Tx¨−ε1−x2x˙+αx+βx3−Fcosωtφjtdt=1605ε+3βhcj+13+160cj5ε+6βhcj+12+160cj25ε+9βh−ε2+αh6+1hcj+1+160h85cjαh2−3+3βh2cj3−30Fω2sin2hω2coshjω+cj−1hcj29βh−5ε+10αh2+3εh+6+hcj−133βh−5ε+hcjcj−126βh−5ε=0,forj=1,2,3,….We may solve this system recursively using the Tartaglia formula for the cubic. The initial data are(59)c0=0,c1=nx˙0T.


## 3.6. Conservative Nonlinear Oscillators

Conservative single-degree-of-freedom nonlinear oscillators are modelled by second-order autonomous ordinary differential equations of the form(60)x¨+Fx=0,x0=x0,x′0=0.Using Chebyshev polynomials or minimization techniques, we may approximate the functionFx by means of a cubic polynomial, say Fx=α+βx+γx2+δx3, and this allows us to study i.v.p. (44) using the solution to the i.v.p.(61)x¨+α+βx+γx2+δx3=0,x0=x0,x′0=0.For example, let us consider a mass attached to two stretched elastic springs [8]. For this problem, the function Fx has the form(62)Fx=x−λx1+x2,0<λ≤1.Suppose thatx≤A. We have Fx≈px+qx3, where(63)p=1+15λ20A2+21sinh−1A−3AA2+12A2+732A5,q=35λAA2+12A2+15−34A2+5sinh−1A32A7.Then, the problem reduces to that of solving a Duffing equation.


## 4. Analysis and Discussion

We have described the way to solve strongly nonlinear oscillators by means of the Galerkin method. In general, any second-order ordinary differential equation may be solved using the Runge–Kutta numerical method. In general, the Runge-Kutta numerical solution offers a more accurate than the Galerkin solution. On the other hand, the Galerkin method offers the possibility to write the solution as a linear combination of hat functions. In this sense, the Galerkin method is a kind of analytical method. For a given conservative oscillator (45), we may approximate the function Fx by a cubic polynomial and then we replace the original problem with problem (45), which has exact analytical solution. However, the exact solution demands the evaluation of either a Jacobian or elliptic Weierstrass function, which has extra costs.


## 5. Conclusions

The Galerkin method offers a way to obtain semianalytical solution to a given nonlinear oscillator of special form (1). The exact solution to it is not known in general. There are other ways to solve it. Perturbative methods like the Lindstedt–Poincaré method or the Krylov–Bogoliubov–Mitropolsky method are also possible for this end. The advantage of the proposed method consists of the possibility to solve the Galerkin equations iteratively using an algebraic equation. The Galerkin hat method is usually used to solve second-order linear ODEs. In this work, we extended it to a class of nonlinear oscillators. Other methods for solving nonlinear differential equations may be found in [8–10].



---

*Source: 8141227-2022-09-07.xml*