Finding periodic ODE solution via small parameter method

Question

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.

Please help me to do it. Thanks a lot for you help!

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.

Can you please fix your code? Thanks a lot!

Answer 1

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[0][t] + x[0]''[t], 
    3*x[1][t] - x[0]'[t]^2 + x[1]''[t], 
    (3*x[2][t] - 2*x[0]'[t]*x[1]'[t] + x[2]''[t]} *)

Solving the first equation yields,

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

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

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

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

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

And then the third equation.

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

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

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

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

Flatten@DSolve[(eq2μ[[4]] == 0) /. x[0][t] -> s2[0] /. 
    x[1][t] -> s2[1] /. x[2][t] -> s2[2], x[3][t], t] // FullSimplify;
s2[3] = TrigReduce[x[3][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.