Simplifying general solutions of differential equations (driven harmonic oscillator)

Question

Solving general differential equations in Mathematica usually leads to somewhat unsightly results.

As an example, consider the solution of the driven, damped harmonic oscillator:

eqn = x''[t] + β x'[t] +  ω0^2 x[t] == f0/m Exp[I ωd t];
s = DSolve[eqn, x[t], t]

Using FullSimplify helps to reduce this mess, but the result is still far away from something an engineer or physicist would recognize at first glance:

Simplify[s]

How can you transform such solutions into something nicer?

Answer 1

Referring to your own answer there is a much simpler form that you may use, recalling:

Collect[expr, var, h] applies h to the expression that forms the coefficient of each term obtained.

sol = x[t] /. s[[1]];

Collect[sol, _C, Simplify]

(E^(I t ωd) f0)/(m ω0^2 + I m β ωd - m ωd^2) + 
 E^(-(1/2) t (β + Sqrt[β^2 - 4 ω0^2])) C[1] + 
 E^(1/2 t (-β + Sqrt[β^2 - 4 ω0^2])) C[2]

---

Martin commented:

It's unfortunate that FullSimplify doesn't do this step on its own. The leaf-count is actually lower after collecting on the C[i].

One can add manipulations using TransformationFunctions therefore:

collectC[x_] := Collect[x, _C, Simplify]
SetOptions[FullSimplify, TransformationFunctions -> {Automatic, collectC}];

FullSimplify[sol]

(E^(I t ωd) f0)/(m ω0^2 + I m β ωd - m ωd^2) + 
 E^(-(1/2) t (β + Sqrt[β^2 - 4 ω0^2])) C[1] + 
 E^(1/2 t (-β + Sqrt[β^2 - 4 ω0^2])) C[2]

Answer 2

One way to tackle the problem is to recognize that the solution is one giant sum with three terms.

You can then convert this sum to a list, simplify the terms individually, and sum all elements of the list back together.

sol = x[t] /. s[[1]];
Total[ FullSimplify[ Apply[ List , sol  ] ] ]

But this "hack" is somewhat unelegant, because it might fail if the inhomogeneous solution is a sum in it's own.

A more robust solution is to Simplify first, then Collect all terms that contain a C[ ]-factor, and finally simplifying the individual terms in that sum:

sol = x[t] /. s[[1]];
Total [Simplify[Apply[List, Collect[ Simplify[ sol ], _C  ] ] ] ]

It's a bit hard to see what is going in that verbose form, so here is the same code in postfix notation (think "Unix pipes"):

sol = x[t] /. s[[1]];
sol // Simplify // Collect[#, _C]&   // Apply[List, #]&   // Simplify // Total

The output of all these commands is equivalent.

Addendum: I did another FullSimplify of the last result and noticed that it stayed the same (instead of reverting to the original). And indeed, checking the complexity of the expressions revealed 121 for the original "FullSimplified"-version, and 90 for the "Collected" version.

I suspect that FullSimplify somehow misses this simplification. Hopefully it will come in future Mathematica versions.