Vtyjkyuk

This is a follow-up question to What is the time period of an oscillator with varying spring constant? . In there a curious young student dared to ask uneasy questions and has been silenced by a more sober and older user. Since I asked myself the same kind of question again and again when I was an undergraduate student and never got any satisfactory answer to it despite all the years I spent in research let me reword the question again.

Question:

Can we find analytical solutions to a harmonic oscillator with a time dependent spring constant?
In mathematical terms we seek solutions two a linear second order ODE with coefficients depending on time only.
beginequation
ddotx(t) + omega^2(t) cdot x(t) = 0
endequation
subject to initial conditions $x(0)=x_0$ and $dotx(0)=v_0$.
We seek a class of functions $omega(t)$ for which the following equation has closed form solutions.

My approach to tackling this problem:

From my scarce mathematical knowledge (due to my poor upbringing or other factors ..) I learned the following statement:

Quotation:
Every linear second order differential equation with three singular point can be reduced to the Gaussian hypergeometric equation by a change of variables (see https://en.wikipedia.org/wiki/Hypergeometric_function).
End of quotation.

This would mean that there should be a wide class of functions $omega(t)$ (i.e. those with not more than three singularities) for which it should be possible to reduce our ODE above to the Gaussian hypergeometric equation. Unfortunately aside from the cases described in here
https://math.stackexchange.com/q/673737/ and in https://math.stackexchange.com/q/1018897/, which I devised myself after tedious work and which indeed is a solution in the form of a hypergeometric function, I am not aware of any other nontrivial functions $omega(t)$ that would yield a solution.
Can anyone enlighten me as to how to construct the relevant change of variables or if it is not possible can anyone explain to me how is that statement in the quotation above proven?

I have recently found a very promising approach [1] to tackling those questions and I thought I would present it in here. Firstly let me describe the algorithm itself and then apply it to our problem and find a whole class of closed form solutions for some functional form of $omega(t)$.

I will change the notation slightly just to be inline with that in the paper.
The idea is to seek for solutions our $bf input$ equation:
beginequation
y^''(x) = v(x)cdot y(x)
endequation
in the following form $y_1,2(x) = m(x) cdot F_1,2(xi(x))$ (ansatz) where $m(x)$ and $xi(x)$ are functions to be specified later and $F_1,2(cdot)$ are fundamental solutions to some other second order ODE to be termed the $bf target$ equation. In other words we have:
beginequation
left(fracd^2d x^2 + a_1(x) fracdd x + a_0(x)right) F_1,2(x)=0
endequation
Now by inserting the ansatz into the input equation and equating the coefficients at $F$ and at $F^'$ to zero we end up with following equations for the functions in question. We have:
begineqnarray
fracm^'m&=& frac12left(
a_1 cdot xi^' -
fracxi^''xi^'
right)\
0&=&left( a_1^2+ 2 a_1^' - 4 a_0right) cdot [xi^']^4 + 3 [xi^'']^2-2 xi^''' cdot xi^' - 4 v cdot [xi^']^2
endeqnarray
The idea is now to search for solutions to the bottom equation above in form of rational functions.
By matching the degree of that rational function (difference in degrees of the numerator and denominator) one can deduce the form of the denominator and also the degree of the numerator and then by inserting those into the equation get some non-linear equations for parameters , equations that are easily solved. I refer the reader to the paper below for details. This method seems to be quite powerful provided we properly chose the target equation. It is conceivable that the target equation should be related to the hypergeometric equation.

Now, to demonstrate the power of this method let us take $v(x) = a cdot x^2 n + b cdot x^n-1$ where $nne-1$. Then it appears that the fundamental solutions to our input equation read:
begineqnarray
left
frac1sqrtx^n M_frac-b(2n+2) sqrta,frac12n+2left( frac2n+1 sqrta x^n+1right),
frac1sqrtx^n W_frac-b(2n+2) sqrta,frac12n+2left( frac2n+1 sqrta x^n+1right)
right
endeqnarray
where $M$ and $W$ are the Whittaker functions (which in turn are related to the confluent hypergeometric equation). Below is a quick and dirty proof of this result

In[648]:= a =.; b =.; n =.; x =.;
FullSimplify[(D[#, x, 2] - (a x^(2 n) + b x^(n - 1)) #) & /@ 1/
 Sqrt[x^n] WhittakerM[-b/((2 n + 2) Sqrt[a]), 1/(2 n + 2), 
 2/(n + 1) Sqrt[a] x^(n + 1)], 
 1/Sqrt[x^n] WhittakerW[-b/((2 n + 2) Sqrt[a]), 1/(2 n + 2), 
 2/(n + 1) Sqrt[a] x^(n + 1)]]

Out[649]= 0, 0

Update: I have been playing around with this algorithm and to my great surprise I realized that the solution given above is a special case of another broad class of solutions below. Let us take $v(x) = a cdot x^2 n + b cdot x^n-1 + c cdot x^-2$ where $n ne -1$. Then the fundamental solutions to our input equation read:
begineqnarray
left
frac1sqrtx^n M_frac-b(2n+2) sqrta,fracsqrt1+4 c2n+2left( frac2n+1 sqrta x^n+1right),
frac1sqrtx^n W_frac-b(2n+2) sqrta,fracsqrt1+4 c2n+2left( frac2n+1 sqrta x^n+1right)
right
endeqnarray
Again, the line below "proves" this result. We have:

In[1249]:= a =.; b =.; c =.; x =.; n =.;
FullSimplify[(D[#, x, 
 2] - (a x^(2 n) + b x^(n - 1) + c x^-2) #) & /@ 1/Sqrt[x^n]
 WhittakerM[-(b/((2 n + 2) Sqrt[a])), Sqrt[
 1 + 4 c]/(2 n + 2) , (2 Sqrt[a])/(n + 1) x^(n + 1)], 
 1/Sqrt[x^n]
 WhittakerW[-(b/((2 n + 2) Sqrt[a])), Sqrt[
 1 + 4 c]/(2 n + 2) , (2 Sqrt[a])/(n + 1) x^(n + 1)]]


Out[1250]= 0, 0

Having said all this it would be interesting to know if this result is known and if yes by what means has it been derived.

[1] M. Bronstein & S. Lafaille, ``Solutions of linear ordinary differential equations in terms of special functions'', Proceedings of ISSAC'2002, Lille, ACM Press, 23-28. https://www-sop.inria.fr/cafe/Manuel.Bronstein/