Using Fourier Series to acquire Nonlinear ODE Periodic Solutions

Direct solution of the last equation in the question also is feasible, because the Fourier series converges very rapidly. as will be seen below. The equation for a three term expansion can be written as

Sum[a[n] (n w)^6 Cos[n w t], {n, 1, 5, 2}]^3] - Sum[a[n] Cos[n w t], {n, 1, 5, 2}] == 0

Reducing the cubic terms in terms of Cos[m w t] shows that the equation contains terms as high as Cos[15 w t] (as would be expected).

Collect[TrigReduce[Sum[a[n] (n w)^6 Cos[n w t], {n, 1, 5, 2}]^3] - 
    Sum[a[n] Cos[n w t], {n, 1, 5, 2}], _Cos, Simplify];

Consistent with the left side of the original equation, only terms up to Cos[5 w t] should be retained. Extract their coefficients and equate them to zero to provide three simultaneous equations determining the three a.

Thread[Coefficient[%, Cos[# w t] & /@ Range[1, 5, 2]] == 0]
(* {1/4 (3 w^18 a[1]^3 + 2187 w^18 a[1]^2 a[3] + 24911296875 w^18 a[3]^2 a[5] + 
    2 a[1] (-2 + 3 w^18 (531441 a[3]^2 + 11390625 a[3] a[5] + 244140625 a[5]^2))) == 0, 
    -a[3] + 1/4 w^18 (a[1]^3 + 68343750 a[1] a[3] a[5] + a[1]^2 (4374 a[3] + 46875 a[5]) + 
    2187 (531441 a[3]^3 + 488281250 a[3] a[5]^2)) == 0, 
    -a[5] + 3/4 w^18 (531441 a[1] a[3]^2 + a[1]^2 (729 a[3] + 31250 a[5]) + 
    15625 (1062882 a[3]^2 a[5] + 244140625 a[5]^3)) == 0} *)

Solutions for the three a then are obtained (slowly) by.

sol = Simplify[Solve[%, {a[1], a[3], a[5]}, Reals], 
    w > 0 && a[1] != 0 && a[3] != 0 && a[5] != 0];

Although most of the seven solutions are truly enormous in LeafCount, insight can be gained from sample numerical evaluations.

sol /. w -> 1.
(* {{a[1] -> 0, a[3] -> 0., a[5] -> 0.}, 
    {a[1] -> 0, a[3] -> 0., a[5] -> -5.91207*10^-7}, 
    {a[1] -> 0, a[3] -> 0., a[5] -> 5.91207*10^-7}, 
    {a[1] -> 0, a[3] -> -0.0000586649, a[5] -> 0.}, 
    {a[1] -> 0, a[3] -> 0.0000586649, a[5] -> 0.}, 
    {a[1] -> -1.23839, a[3] -> 0.000315137, a[5] -> -5.77387*10^-6}, 
    {a[1] -> 1.23839, a[3] -> -0.000315137, a[5] -> 5.77387*10^-6}} *)

The first result is trivial, and the next four results correspond to solutions for, in effect, harmonics of w == 1. Only the last two solutions are of interest, and one is the negative of the other. So, we can restrict attention to Last@sol. Although the values of a vary rapidly with w (shown here for a[1]),

ListLogPlot[(a[1] /. Last[sol] /. w -> #) & /@ Range[.001, 10.002, .5], 
    DataRange -> {0, 10}, AxesLabel -> {w, a[1]}, LabelStyle -> {Bold, Medium}]

the ratios {a[[2]]/a[[1]], a[[3]]/a[[1]]} are essentially constant at {-0.000254473, 4.6624*10^-6} except for w < 0.01. Thus, a[1] Cos[w t] typically is a very good approximate solution to the ODE. This assertion can be further substantiated by comparing the single-cosine approximation with the solution from the earlier symbolic solution.

GraphicsRow[Show[
    ParametricPlot[{{qp - Sqrt[fp[[1]]], x}, {qp + Sqrt[fp[[1]]], -x}, 
        {3 qp - Sqrt[fp[[1]]], -x}, {3 qp + Sqrt[fp[[1]]], x}} /. c -> #, 
        {x, 0, 100}, AspectRatio -> 1/GoldenRatio, PlotPoints -> 500], 
    Plot[(2 c/3)^(3/4) Cos[t Pi/(2 qp)] /. c -> #, {t, 0, 4 qp /. c -> #}, 
        PlotStyle -> Dashed]] & /@ {1/100, 1, 100}, 
    ImageSize -> 800]

Agreement is excellent for c in the range of 1/100 to 100.

This problem can be solved symbolically as follows. First, note that the first integral of the first ODE is given for x[t] > 0 by

x'[t]^2 + x[t]^(4/3)/(2/3) == c

where c is a constant of integration. Proof:

Simplify[D[x'[t]^2 + x[t]^(4/3)/(2/3) - c, t]/(2 x'[t])]
(* x[t]^(1/3) + x''[t] *)

Similarly, for x[t] < 0 the first integral is

x'[t]^2 + (-x[t])^(4/3)/(2/3) == c

The resulting phase-space plot for c == Range[2, 40, 2] is

Show[ContourPlot[xp^2 + x^(4/3)/(2/3), {x, -5, 5}, {xp, -5, 5}, 
         Contours -> Range[2, 40, 2], ContourShading -> None], 
     ContourPlot[xp^2 + (-x)^(4/3)/(2/3), {x, -5, 5}, {xp, -5, 5}, 
         Contours -> Range[2, 40, 2], ContourShading -> None], 
     FrameLabel -> {x, x'}, LabelStyle -> {Bold, Medium}]

To obtain the actual symbolic solution for x[t] > 0, use

DSolve[x''[t] + x[t]^(1/3) == 0, x, t]
(* Solve[(1/(2 C[1] - 3 x[t]^(4/3))) Sqrt[6] C[1]^(3/2) 
   (EllipticE[I ArcSinh[(3/2)^(1/4) Sqrt[-(1/Sqrt[C[1]])] x[t]^(1/3)], -1] - 
    EllipticF[I ArcSinh[(3/2)^(1/4) Sqrt[-(1/Sqrt[C[1]])] x[t]^(1/3)], -1])^2 
   (4 - (6 x[t]^(4/3))/C[1]) == (t + C[2])^2, x[t]] *)

If Mathematica could Solve this expression, it would have. So, extract and Simplify the unsolved equation.

f = Simplify[%[[1]] /. {C[1] -> c, C[2] -> 0, x[t] -> x}, c > 0]
(* 2 Sqrt[6] Sqrt[c] (EllipticE[ArcSin[((3/2)^(1/4) x^(1/3))/c^(1/4)], -1] - 
                      EllipticF[ArcSin[((3/2)^(1/4) x^(1/3))/c^(1/4)], -1])^2 == t^2 *)

The corresponding solution for x[t] < 0 is obtained by replacing x by -x. Note that C[2] is the initial value of t and has be set to zero without loss of generality. The plot of these expressions for c -> 2 (corresponding to the innermost curve in the plot above) is

ParametricPlot[{{Sqrt[f[[1]] /. c -> 2], x}, {-Sqrt[f[[1]] /. c -> 2], x}, 
                {Sqrt[f[[1]] /. c -> 2], -x}, {-Sqrt[f[[1]] /. c -> 2], -x}}, 
    {x, 0, 2}, PlotStyle -> Black, AxesLabel -> {t, x}, LabelStyle -> {Bold, Medium}]

The turning point is given by

Solve[First@Cases[f[[1]], ArcSin[z_] -> z, Infinity, 1] == 1, x][[1, 1]] 
(* x -> (2/3)^(3/4) c^(3/4) *)

and the quarter period of the solution by

qp = Sqrt[f[[1]] /. %]
(* 2^(3/4) 3^(1/4) c^(1/4) (EllipticE[-1] - EllipticK[-1]) *)

For c -> 2, the turning point and full period (4 qp) are 1.24081 and 6.30736, consistent with the second plot.