# Using Fourier Series to acquire Nonlinear ODE Periodic Solutions

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$?

• Note that \beta depends on n. – mmeent May 11 '17 at 9:39
• True. I will correct it, thanks! – Mitscaype May 11 '17 at 9:40

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^3 + 2187 w^18 a^2 a + 24911296875 w^18 a^2 a +
2 a (-2 + 3 w^18 (531441 a^2 + 11390625 a a + 244140625 a^2))) == 0,
-a + 1/4 w^18 (a^3 + 68343750 a a a + a^2 (4374 a + 46875 a) +
2187 (531441 a^3 + 488281250 a a^2)) == 0,
-a + 3/4 w^18 (531441 a a^2 + a^2 (729 a + 31250 a) +
15625 (1062882 a^2 a + 244140625 a^3)) == 0} *)


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

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


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

sol /. w -> 1.
(* {{a -> 0, a -> 0., a -> 0.},
{a -> 0, a -> 0., a -> -5.91207*10^-7},
{a -> 0, a -> 0., a -> 5.91207*10^-7},
{a -> 0, a -> -0.0000586649, a -> 0.},
{a -> 0, a -> 0.0000586649, a -> 0.},
{a -> -1.23839, a -> 0.000315137, a -> -5.77387*10^-6},
{a -> 1.23839, a -> -0.000315137, a -> 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),

ListLogPlot[(a /. Last[sol] /. w -> #) & /@ Range[.001, 10.002, .5],
DataRange -> {0, 10}, AxesLabel -> {w, a}, LabelStyle -> {Bold, Medium}] the ratios {a[]/a[], a[]/a[]} are essentially constant at {-0.000254473, 4.6624*10^-6} except for w < 0.01. Thus, a 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[]], x}, {qp + Sqrt[fp[]], -x},
{3 qp - Sqrt[fp[]], -x}, {3 qp + Sqrt[fp[]], 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.

• Amazing job. What's more you have been explaining the whole way what you are doing, thank you. I have a question, why the plot of the solution from your previous answer regarding the symbolic solution does not coincide with the ones at the end of your Fourier analysis? – Mitscaype May 14 '17 at 17:34
• Also, I get an error about the "qp" parameter in your last Parametric plot:/ – Mitscaype May 14 '17 at 18:25
• @Mitscaype Actually, the plots are the same, if account is taken of the periodicity of the solution. You encountered the qp error, because I neglected to define qp. Now, I have done so. See the end of the earlier answer for its definition. – bbgodfrey May 14 '17 at 22:47
• @Ahh, I see. Therefore, if changing the period to $\omega=5.30304$ which corresponds to the solution of $x(0)=1, \dot{x}(0)=0$ will give me the correct results, since you showed and I verified that the ratio $A_3/A_1, A_5/A_1$ is a constant. Moreover, on your previous reply, $C$ is $x(0)$, right? – Mitscaype May 15 '17 at 4:31
• @Mitscaype The amplitude of the oscillatory solution is given by (2/3)^(3/4) c^(3/4), not by c. See the second to the last equation in my earlier answer. – bbgodfrey May 15 '17 at 13:40

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 - 3 x[t]^(4/3))) Sqrt C^(3/2)
(EllipticE[I ArcSinh[(3/2)^(1/4) Sqrt[-(1/Sqrt[C])] x[t]^(1/3)], -1] -
EllipticF[I ArcSinh[(3/2)^(1/4) Sqrt[-(1/Sqrt[C])] x[t]^(1/3)], -1])^2
(4 - (6 x[t]^(4/3))/C) == (t + C)^2, x[t]] *)


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

f = Simplify[%[] /. {C -> c, C -> 0, x[t] -> x}, c > 0]
(* 2 Sqrt 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 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[] /. c -> 2], x}, {-Sqrt[f[] /. c -> 2], x},
{Sqrt[f[] /. c -> 2], -x}, {-Sqrt[f[] /. c -> 2], -x}},
{x, 0, 2}, PlotStyle -> Black, AxesLabel -> {t, x}, LabelStyle -> {Bold, Medium}] The turning point is given by

Solve[First@Cases[f[], 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[] /. %]
(* 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.

• Thank you for taking the time to provide an answer. Nevertheless, I have some questions. At the part where you use {DSolve} you have some Elliptic functions involved. How did this come up? Also, the way I understand it, is that you did not took at all the Fourier approximation of $x(t)$, but you found period in a different way, is that correct? – Mitscaype May 12 '17 at 5:04
• @Mitscaype The Elliptic functions are part of the solution obtained by DSolve. They are not uncommon in ODEs that can be reduced to first integrals of the motion. I found the period by computing the point at which x'[t] vanished. Given the results in my answer, it would not be difficult to obtain the corresponding Fourier coefficients, but I felt that you would not need them, given the complete answer. I can provide answers to other questions, if any, tomorrow. It is late here now. – bbgodfrey May 12 '17 at 5:10
• So if I understand it correctly, the fact that you could compute a period means that there are indeed periodic solutions to the problem, right? I really would like the procedure to acquire the Fourier $A_n$ along with $\omega$ and then plot this solution to the one acquired by DSolve. Again, many thanks for your time and help. Anytime tomorrow is fine :) – Mitscaype May 12 '17 at 5:20
• @Mitscaype The solutions provided in my answer all are periodic, and all periods (or frequencies) are allowed. Because the equation is nonlinear, the waveform depends on the frequency. I shall look at the Fourier amplitudes question later today. – bbgodfrey May 12 '17 at 12:45
• @Mitscaype Cases[f[], ArcSin[z_] -> z, Infinity, 1] extracts the argument of ArcSin from the left side of f. That argument must lie between -1 and 1 for the ArcSin to have a real value. Thus turning points occur when the argument of ArcSin has one of those two values. I chose 1 for convenience. – bbgodfrey Jul 5 '17 at 14:30