# Finding periodic ODE solution via small parameter method

I want to figure out how to find periodic solutions of ODE via the small parameter method. I will provide couple of examples of what I mean.

1. Consider equation $\ddot x + 3x = 2 \sin t + \mu \dot x^2$. The solution I want to get is $x(t, \mu) = \sin t + \mu (\frac16 - \frac12 \cos 2t) + \mu^2(-\frac16 \sin 3t + \frac12 \sin t)) + ...$.

2. Consider equation $\ddot x + 3x + x^3= 2 \mu \cos t$. The solution I want to get is $x(t, \mu) = \mu \cos t + \frac{\mu^3}{8} (\frac13 \cos t - 3 \cos t) + ...$.

I am interested in small numbers of first coefficients, namely 3-4 first terms.

This must be implemented in Mathematica, but I was not able to find anything relevant.

Update:

bbgodfrey, thanks a lot for your answer! I will post my answer to your answer here since comments are limited. Your answer is wrong, unfortunately, but it is almost correct. I can solve this kind of problems manually and I see the errors just in first two lines, but I can not fix them, since I am Mathematica newbie. I can also explain how to solve this kind of problems and what you should get. I hope that as long as first two lines are correct, everything will work correct. Anyway I can find any error in the code, but I can not fix it. Let us consider the first example. What your are doing is you trying to find the solution in the form $$x(t, \mu)=x_0(t)+\mu x_1(t)+\mu^2x_2(t)+...$$ You just plug first terms of this series into the equation and make the coefficients of $\mu$ with the same degree equal. What you should get in first equation is $$\ddot x_0 + 3x_0 = 2 \sin t, \ddot x_1 + 3x_1 = \dot x_0^2,\ddot x_2 + 3x_2 = 2\dot x_0\dot x_1.$$ But what your code gives in last condition is $$\ddot x_2 + 3x_2 = 4\dot x_0\dot x_1.$$ This $4$ is incorrect, it must be 2. You can check it manually within $3$ mins.

• For future reference, Mathematica input is far more valuable in a question than formatted HTML or TeX. – Daniel Lichtblau May 28 '17 at 15:23

Edit: Corrected error identified by Hedgehog

I am not aware of a Mathematica function to solve the two equations by expansion in μ. Nonetheless, solutions can be obtained as follows.

ODE 1

eq1 = D[x[t, μ], {t, 2}] + 3 x[t, μ] == 2 Sin[t] + μ D[x[t, μ], t]^2;


First, expand eq1 to second order in μ to obtain three ODEs

eq1μ = Take[CoefficientList[Subtract @@ eq1 /.
x -> Function[t, Sum[μ^i x[i][t], {i, 0, 2}]], μ], 3]
(* {-2*Sin[t] + 3*x[t] + x''[t],
3*x[t] - x'[t]^2 + x''[t],
(3*x[t] - 2*x'[t]*x'[t] + x''[t]} *)


Solving the first equation yields,

Flatten@DSolve[eq1μ[] == 0, x[t], t] // FullSimplify
(* {x[t] -> C Cos[Sqrt t] + Sin[t] + C Sin[Sqrt t]} *)


Evidently, the OP wishes the inhomogeneous solutions only. So, set the constants of integration to zero.

s1 = x[t] /. (% /. _C -> 0)
(* Sin[t] *)


Proceed now to the second equation, replacing x by the expression just obtained.

Flatten@DSolve[(eq1μ[] == 0) /. x'[t] -> D[s1, t],
x[t], t] // FullSimplify
s1 = x[t] /. (% /. _C -> 0)
(* {x[t] -> 1/6 - 1/2 Cos[2 t] + C Cos[Sqrt t] + C Sin[Sqrt t]} *)
(* 1/6 - 1/2 Cos[2 t] *)


And then the third equation.

Flatten@DSolve[(eq1μ[] == 0) /. x'[t] -> D[s1, t] /.
x'[t] -> D[s1, t], x[t], t] // FullSimplify
s1 = TrigReduce[x[t] /. (% /. _C -> 0)]
(* {x[t] -> C Cos[Sqrt t] + (2 Sin[t]^3)/3 + C Sin[Sqrt t]} *)
(* 1/6 (3 Sin[t] - Sin[3 t]) *)


Combining these expressions yield the desired second-order expansion.

sol1 = Sum[μ^i s1[i], {i, 0, 2}]
(* μ (1/6 - 1/2 Cos[2 t]) + Sin[t] + 1/6 μ^2 (3 Sin[t] - Sin[3 t]) *)


ODE 2

eq2 = D[x[t, μ], {t, 2}] + 3 x[t, μ] + x[t, μ]^3 == 2 μ Cos[t]


The second ode is solved much like the first.

Take[CoefficientList[Subtract @@ eq2 /.
x -> Function[t, Sum[μ^i x[i][t], {i, 0, 3}]], μ], 4]
(* {3*x[t] + x[t]^3 + x''[t],
-2*Cos[t] + 3*x[t] + 3*x[t]^2*x[t] + x''[t],
3*x[t]*x[t]^2 + 3*x[t] + 3*x[t]^2*x[t] + x''[t],
x[t]^3 + 6*x[t]*x[t]*x[t] + 3*x[t] + 3*x[t]^2*x[t] + x''[t]} *)


Although DSolve certainly can solve the first of the four equations here, it is easier to observe that the desired inhomogeneous solution is

s2 = 0;


because the first equation itself is homogeneous. Proceeding to the other three equations,

Flatten@DSolve[(eq2μ[] == 0) /. x[t] -> s2, x[t], t] // FullSimplify;
s2 = x[t] /. (% /. _C -> 0)
(* Cos[t] *)

Flatten@DSolve[(eq2μ[] == 0) /. x[t] -> s2 /.
x[t] -> s2, x[t], t] // FullSimplify;
s2 = TrigReduce[x[t] /. (% /. _C -> 0)]
(* 0 *)

Flatten@DSolve[(eq2μ[] == 0) /. x[t] -> s2 /.
x[t] -> s2 /. x[t] -> s2, x[t], t] // FullSimplify;
s2 = TrigReduce[x[t] /. (% /. _C -> 0)]
(* 1/24 (-9 Cos[t] + Cos[3 t]) *)


Thus, the desired solution is

sol2 = Sum[μ^i s2[i], {i, 0, 3}]
(* μ Cos[t] + 1/24 μ^3 (-9 Cos[t] + Cos[3 t]) *)


The corresponding expression in the question contains a typo. (Cos[t]/3 should be Cos[3 t]/3.

The approach used here probably could be packaged without too much difficulty.