Simple Harmonic Motion—Quadratic Perturbation

Author

Eric W. Weisstein
May 3, 2004

This notebook downloaded from http://mathworld.wolfram.com/notebooks/DifferentialEquations/SimpleHarmonicMotionQuadraticPerturbation.nb.

For more information, see Eric's MathWorld entry http://mathworld.wolfram.com/SimpleHarmonicMotionQuadraticPerturbation.html.

©2005 Wolfram Research, Inc. except for portions noted otherwise

Initialization

<<Graphics`Colors`

Solution of the Perturbed ODE

x[t_]:=x0[t]+ϵ  x1[t];

x''[t]+ω0^2 x[t]-α ϵ x[t]^2

c=CoefficientList[x''[t]+ω0^2 x[t]-α ϵ x[t]^2,ϵ]

HoldForm[Evaluate[c Table[ϵ^i,{i,0,3}]]]/.List->Plus

(ϵ^0)part

DSolve[{x0''[t]+((ω0)^2) x0[t]==0,x0[0]==A,x0'[0]==0},x0[t],t]

(ϵ^1)part

-α ((x0[t])^2)+((ω0)^2) x1[t]+x1''[t]==0/.x0->Function[t,A  Cos[t ω0]]

Simplify[TrigReduce//@FullSimplify[DSolve[
-(A^2) α ((Cos[t ω0])^2)+((ω0)^2) x1[t]+x1''[t]==0,
x1[t],t]]]

Simplify[-(A^2) α ((Cos[t ω0])^2)+((ω0)^2) x1[t]+x1''[t]/.x1->Function[t,(((A^2) α)/(2 ((ω0)^2)))+C[1] Cos[t ω0]-(((A^2) α Cos[2 t ω0])/(6 ((ω0)^2)))+C[2] Sin[t ω0]]]

Plots

tmax=100;
f=NDSolve[{Evaluate[x''[t]+ω0^2 x[t]-α ϵ x[t]^2==0 /.{ω0->1,α->1,ϵ->0}],x[0]==0,x'[0]==1},x[t],{t,0,tmax}][[-1,-1,-1]];
g=NDSolve[{Evaluate[x''[t]+ω0^2 x[t]-α ϵ x[t]^2==0 /.{ω0->1,α->1,ϵ->.2}],x[0]==0,x'[0]==1},x[t],{t,0,tmax}][[-1,-1,-1]];

h=NDSolve[{Evaluate[x''[t]+ω0^2 x[t]-α ϵ x[t]^2==0 /.{ω0->1,α->1,ϵ->.0001}],x[0]==0,x'[0]==1},x[t],{t,0,tmax}][[-1,-1,-1]];

dmax=40;
Show[GraphicsArray[{
{Plot[{f,h},{t,0,dmax},PlotStyle->{Black,Red},DisplayFunction->Identity,AxesLabel->TraditionalForm/@{t,x[t]}],
Plot[f-h,{t,0,dmax},DisplayFunction->Identity,AxesLabel->TraditionalForm/@{t,Δ x(t)}]},
{
Plot[{f,g},{t,0,dmax},PlotStyle->{Black,Red},DisplayFunction->Identity,AxesLabel->TraditionalForm/@{t,x[t]}],
Plot[f-g,{t,0,dmax},DisplayFunction->Identity,AxesLabel->TraditionalForm/@{t,Δ x(t)}]
}}]]