8
$\begingroup$

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

Improve this question
$\endgroup$
2
  • 1
    $\begingroup$ Note that \beta depends on n. $\endgroup$ Commented May 11, 2017 at 9:39
  • 1
    $\begingroup$ True. I will correct it, thanks! $\endgroup$ Commented May 11, 2017 at 9:40

2 Answers 2

Reset to default
10
$\begingroup$

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}]

enter image description here

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]

enter image description here

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

Improve this answer
$\endgroup$
17
  • $\begingroup$ 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? $\endgroup$ Commented May 14, 2017 at 17:34
  • $\begingroup$ Also, I get an error about the "qp" parameter in your last Parametric plot:/ $\endgroup$ Commented May 14, 2017 at 18:25
  • $\begingroup$ @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. $\endgroup$ Commented May 14, 2017 at 22:47
  • $\begingroup$ @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? $\endgroup$ Commented May 15, 2017 at 4:31
  • $\begingroup$ @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. $\endgroup$ Commented May 15, 2017 at 13:40
7
$\begingroup$

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}]

enter image description here

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}]

enter image description here

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.

Improve this answer
$\endgroup$
12
  • $\begingroup$ 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? $\endgroup$ Commented May 12, 2017 at 5:04
  • 1
    $\begingroup$ @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. $\endgroup$ Commented May 12, 2017 at 5:10
  • $\begingroup$ 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 :) $\endgroup$ Commented May 12, 2017 at 5:20
  • 1
    $\begingroup$ @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. $\endgroup$ Commented May 12, 2017 at 12:45
  • 1
    $\begingroup$ @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. $\endgroup$ Commented Jul 5, 2017 at 14:30

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.