Using Fourier Series to acquire Nonlinear ODE Periodic Solutions

Question

For the following Cauchy problem of a non linear ODE: \begin{equation} \ddot{x}=-|x|^{1/3} \end{equation} which satisfies the initial conditions $x(0)=1, \dot{x}(0)=0$, I am aware that there are periodic solutions, since it is a conservative system.

To simplify further, it is acceptable to start without the absolute value and solve the corresponding system: \begin{equation} \ddot{x}=-x^{1/3} \end{equation} for the same ICs. I would like therefore to expand this $x(t)$ into Fourier series and acquire an approximation of these periodic solutions.

So far I am aware that Mathematica has a built in function for taking the Fourier series of a function but what I want actually to do is the following:

By considering that: \begin{equation} x(t)=\sum_{n=1}^{\infty}A_{2n-1} \cos ((2n-1)\omega t) \end{equation} and plunging it into the ODE I get the following: \begin{equation} \sum_{n=1}^{\infty}A_{2n-1}\cos \left( (2n-1)\omega t \right)=\left[ \sum_{n=1}^{\infty}\beta A_{2n-1}\cos \left( (2n-1)\omega t \right) \right]^3 \end{equation} where $\omega$ the frequency and $\beta=\left( (2n-1)\omega \right)^2 \in \mathbb{R}^+$.

Now I would like to expand and acquire the coefficients $A_n$.

I am not sure whether this piece of code exists already or if I should try to write it, but for the latter one I do not believe it would be an easy task to do since my coding skills unfortunately are limited.

How can I tell Mathematica to expand and equate the coefficients $A_n$ so that it gives me back their expression? And if there is not a general one (which is highly improbable since I have proved that these orbits exist for these particular ICs) how can I acquire at least the first say m $A_1, \cdots, A_m$, where $m<n$?

Answer 1

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.

Answer 2

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.