# Using Fourier Series to acquire Nonlinear ODE Periodic Solutions

For the following Cauchy problem of a non linear ODE: $$\ddot{x}=-|x|^{1/3}$$ 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: $$\ddot{x}=-x^{1/3}$$ 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: $$x(t)=\sum_{n=1}^{\infty}A_{2n-1} \cos ((2n-1)\omega t)$$ and plunging it into the ODE I get the following: $$\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$$ 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. Commented May 11, 2017 at 9:39
• True. I will correct it, thanks! Commented May 11, 2017 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[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.

• 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? Commented May 14, 2017 at 17:34
• Also, I get an error about the "qp" parameter in your last Parametric plot:/ Commented May 14, 2017 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. Commented May 14, 2017 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[1]$ is $x(0)$, right? Commented May 15, 2017 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. Commented May 15, 2017 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[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.

• 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? Commented May 12, 2017 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. Commented May 12, 2017 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 :) Commented May 12, 2017 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. Commented May 12, 2017 at 12:45
• @Mitscaype Cases[f[[1]], 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. Commented Jul 5, 2017 at 14:30