# 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. – 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[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? – 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[1]$ 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[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? – 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[[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. – bbgodfrey Jul 5 '17 at 14:30