The Duffing Equation with Mathematica

Dr. Gustavo A. Pérez M.
Honduras National Academy of Sciences
gus@pepper.org

Perception is everything

Abstract :
We use Mathematica to obtain the analytical and numerical solutions of Duffing equation.

Introduction

We consider as Duffing problem the differential equation:
u''+u+ϵ (u^3)=0, as initial conditions, we take u(0)=(Subscript[x, 0]) and u'(0)=(Subscript[v, 0]).
The solution u of our problem is a function of the independent variable t and the parameter ϵ that is a dimensionless quantity, and its a measure of the strength of the non linearity. Then, an exact solution has a frequency ω that is a function of ϵ, the non linearity.
Mathematica solves numerically the Duffing problem.

eqdu=u''[t]+ u[t]+ϵ u[t]^3==0

eqdv=v''[t]+ v[t]+ϵ v[t]^3==0

We write the same equation two times in order to solve one for the amplitude eqdu and the other for the velocity eqdv with diferent initial conditions.

The potential energy and the solutions

The Duffing problem, as defined, is periodic and conservative, we keep the parameter ϵ in the necessary range as to satisfy these conditions.
Let define the potential energy:

F[u_,ϵ_]:=(1/2)(u^2)+ (1/4) ϵ (u^4)

Plot3D[F[u,ϵ],{u,-1.5,1.5},{ϵ,-1,1}]

Plot[F[u,-1],{u,-2,2}]

The energy equation

eqce=(1/2)v^2+F[u,ϵ]==h

The relation between amplitude and velocity

eqa=(1/2)(Subscript[v, 0]^2)==((Subscript[x, 0]^2)/2)+((ϵ (Subscript[x, 0]^4))/4)

Solve[eqa,(Subscript[x, 0])]

Where the only valid solution is (Subscript[x, 0])->(Sqrt[-(1/ϵ)+((Sqrt[1+2 ϵ (Subscript[v, 0]^2)])/ϵ)])

Whenever ϵ>0 the velocity is zero in two points and the trajectory is closed and periodic. Whenever ϵ<0 the trajectory could be closed and periodic for some values of h or could be unlimited for values of h beyond the maxima of the potential.
When the trajectory is closed and periodic the initial conditions are in some way limited by the energy conservation and a problem of initial velocity is equivalent after some time to a initial amplitude problem.

The equation of energy conservation eqce could be written as:
(1/2)(v^2)=h-(1/2)(u^2)- (1/4) ϵ (u^4)=h-F(u,ϵ)

eqe=(1/2)(v^2)==h-(1/2)(u^2)- (1/4) ϵ (u^4)

Solve[eqe,v]

If we do u'=(du/dt)=v and integrate:
t=±(∫_(Subscript[x, 0])^u)(1/(Sqrt[2 h-(u^2)-(1/2) ϵ (u^4)]))du
Introducing the amplitude in the energy conservation
h=(1/2)(Subscript[x, 0]^2)+(1/4)ϵ (Subscript[x, 0]^4)
and making the transformation
u=-(Subscript[x, 0])Cos(θ)

This takes some minutes

eq1=t1==(1/(Sqrt[1+ϵ (Subscript[x, 0]^2)]))Integrate[(1/(Sqrt[1-m Sin[θ]^2])),{θ,θ1,θ2},Assumptions->m>0&&m<1&&θ1<θ2]

eq1=t1==((-EllipticF[θ1,m]+EllipticF[θ2,m])/(Sqrt[1+a]))

This is our principal equation and gives the analytical solution to the Duffing problem other name is the first integral of movement.
The value of m is:

m=((ϵ (Subscript[x, 0]^2))/(2 (1+ ϵ (Subscript[x, 0]^2))))

Analytical solution for a initial amplitude

In our transformation a initial amplitude and zero initial velocity is characterized by θ1->0.

S1=Solve[eq1/.θ1->0,θ2]

In this case the solution is:
u[t]=-(Subscript[x, 0])JacobiCn[t1 (Sqrt[1+ϵ (Subscript[x, 0]^2)]),m] and making  a=ϵ (Subscript[x, 0]^2).

ua[t_]=-(Subscript[x, 0])JacobiCN[t (Sqrt[1+a]),m]

Analytical solution for a initial velocity

In this case we find the amplitude equivalent for the velocity

(Subscript[x, 0])=(Sqrt[-(1/ϵ)+((Sqrt[1+2 ϵ (Subscript[v, 0]^2)])/ϵ)])

This solution is characterized by θ1->(π/2)

S2=Solve[eq1/.θ1->(π/2),θ2]

The solution is

uv[t_]=-(Subscript[x, 0])JacobiCN[EllipticK[(a/(2 (1+a)))]+t (Sqrt[1+a]),(a/(2 (1+a)))]

uv[t_]=JacobiSD[t (Sqrt[1+ϵ (-(1/ϵ)+((Sqrt[1+2 ϵ (Subscript[v, 0]^2)])/ϵ))]),((ϵ (-(1/ϵ)+((Sqrt[1+2 ϵ (Subscript[v, 0]^2)])/ϵ)))/(2 (1+ϵ (-(1/ϵ)+((Sqrt[1+2 ϵ (Subscript[v, 0]^2)])/ϵ)))))] (Sqrt[-(1/ϵ)+((Sqrt[1+2 ϵ (Subscript[v, 0]^2)])/ϵ)]) (Sqrt[1-((ϵ (-(1/ϵ)+((Sqrt[1+2 ϵ (Subscript[v, 0]^2)])/ϵ)))/(2 (1+ϵ (-(1/ϵ)+((Sqrt[1+2 ϵ (Subscript[v, 0]^2)])/ϵ)))))])

The Period

Since the time needed to move from u=-(Subscript[x, 0]) to u=(Subscript[x, 0]) is one half the period then

T=(2/(Sqrt[1+ϵ (Subscript[x, 0]^2)]))Integrate[(1/(Sqrt[1-m Sin[θ]^2])),{θ,0,π},Assumptions->m<1]

Ta[m_,a_]=T/.ϵ (Subscript[x, 0]^2)->a

Numerical Solutions

Mathematica solves the numerical Duffing problem, here we obtain the numerical solutions and its comparison with the analytical results.

(Subscript[x, 0])=1;ϵ=1;(Subscript[v, 0])=1;

a=ϵ (Subscript[x, 0]^2);

This is the numerical solution for (Subscript[x, 0])=1

m=((ϵ (Subscript[x, 0]^2))/(2 (1+ ϵ (Subscript[x, 0]^2))))

sol1=NDSolve[{eqdu/.ϵ->1,u[0]==-1,u'[0]==0},u,{t,0,5}]

The form of the solution for the amplitude (Subscript[x, 0]) is:

Plot[ua[t],{t,0,5}]

Note here that even if the solution is a Cosine the frequency is a function of the non lineality ϵ.

Comparison with the analytical solution

Plot[ua[t]-Evaluate[u[t]/.sol1],{t,0,5}]

The numerical solution for (Subscript[v, 0])

(Subscript[x, 0])=(Sqrt[-(1/ϵ)+((Sqrt[1+2 ϵ (Subscript[v, 0]^2)])/ϵ)])

m=((ϵ (Subscript[x, 0]^2))/(2 (1+ ϵ (Subscript[x, 0]^2))))

sol2=NDSolve[{eqdv/.ϵ->1,v[0]==0,v'[0]==1},v,{t,0,5}]

The form of the solution for (Subscript[v, 0]) is:

Plot[uv[t],{t,0,5}]

Plot[uv[t]-Evaluate[v[t]/.sol2],{t,0,5}]

Values of the Exact Period

The values of the exact period could be tabulated as a function of m and a the non linearity.

mt={0,0.02,0.04,0.06,0.08,0.1}

at={0,0.042,0.087,0.136,0.190,0.25}

Ta[mt,at]

Grid[{mt,at,Ta[mt,at]},Frame->All]

This proves that the period and frequency depend on the non linearity.

Conclusions

We use Mathematica to find the analytical solutions of the Duffing problem without friction. With the numerical capabilities of Mathematica,  we check our solutions. 

The Jacobi Elliptic functions  result of the variable transformation and the solution of the first integral for θ2,  that is a Jacobi Amplitude.

Mathematica could be used to bring the students of the level of the Nayfeh book to use the Jacobi Elliptic functions without effort.

References

Nayfeh, Ali Hasan. Introduction to Perturbation Techniques. John Wiley & Sons 1993. Page 113 to 118.
Shivamoggi, Bhimsen K. Perturbation Methods for Differential Equations. Birkhausser Boston 2003. Page 25 to 28.
Introduction to the Jacobi elliptic functions. http://functions.wolfram.com/EllipticFunctions/JacobiNC/introductions/JacobiPQs/ShowAll.html