What is the time period of an oscillator with varying spring constant? It is well known that the time period of a harmonic oscillator when mass $m$ and spring constant $k$ are constant is $T=2\pi\sqrt{m/k}$. 
However, I would be interested to know what the time period is if $k$ is not constant. I have searched hours after hours for right answers from Google and came up with nothing. I am looking for an analytical solution. 
 A: Here is a solution for a spring force that varies directly with displacement. It thus varies with time implicitly, but has no explicit dependence on time or any other variable.
Givens and Assumptions


*

*oscillator with mass $m$

*amplitude of oscillation $A$

*oscillator displacement, $x$, varies with time, but $x(t)$ is unknown

*spring applies force varying with displacement, $F(x)$

*The function $F(x)$ is an odd function, that is $F(-x) = -F(x)$ (otherwise the amplitude could be different in the positive and negative directions - see below for what to do in this case)

*equilibrium position is $x=0$, that is $F(0) = 0$ (for convenience only)


Objective
Find the period of oscillation, $T$
Solution
Starting from conservation of energy, the sum of the kinetic and potential energy of the mass must be equal to the total energy, which is constant.
$$KE(x)+PE(x)=E$$
$$KE(x)=\frac{1}{2}mv^2(x)$$
$$PE(x)=\intop_0^xdx'\,F(x')$$
So $PE(x)$ is the potential energy stored in the spring, with $x'$ as just an integration variable.
We can think of $PE(x)$ as another way of defining the force-displacement relationship of the spring. We can define the potential energy versus displacement or the force versus displacement, and getting the other one is fairly easy.  
Now, at $x=A$, $KE(x=A)=0$, so $PE(A)=E$ is known.  
And so we have
$$\frac{1}{2}mv^2(x)=PE(A)-PE(x)$$
Solving for $v(x)$,
$$v(x)=\sqrt{\frac{2}{m}\left(PE(A)-PE(x)\right)}$$
Because $v=\frac{dx}{dt}$, we can also write
$$dt=\frac{dx}{v(x)}$$
Integrating both sides, the time to go from a position $x_0$ to $x_1$ is
$$\Delta t = \intop_{x_0}^{x_1}\frac{dx}{v(x)}$$
In particular, we know the time required to go from $x=0$ to $x=A$ is $T/4$, so
$$T=4\intop_0^A\frac{dx}{v(x)}$$
$$T=4\intop_0^A\frac{dx}{\sqrt{\frac{2}{m}}\sqrt{PE(A)-PE(x)}}$$
which further simplifies to...
Final Result
$$T=\sqrt{8m}\intop_0^A\frac{dx}{\sqrt{PE(A)-PE(x)}}$$
Check of Result
For the linear case, $F(x)=kx$, so $PE(x)=\frac{1}{2}kx^2$ and $PE(A)=\frac{1}{2}kA^2$, which gives
$$T=\sqrt{8m}\intop_0^A\frac{dx}{\sqrt{\frac{k}{2}}\sqrt{A^2-x^2}}
=4\sqrt{\frac{m}{k}}\intop_0^A\frac{dx}{\sqrt{A^2-x^2}}$$
This integral can be looked up in a table, to obtain
$$T=4\sqrt{\frac{m}{k}}\left(\sin^{-1}(1)-sin^{-1}(0)\right)=4\sqrt{\frac{m}{k}}\frac{\pi}{2}$$
$$T=2\pi\sqrt{\frac{m}{k}}$$
as expected. (QED)

We can dispense with the assumption that $F(x)$ is odd if we define two amplitude values: $A_+ > 0$ for the amplitude in the positive direction and $A_- < 0$ for the amplitude in the negative direction.
The total oscillator energy, $E = PE(A_+) = PE(A_-)$, so we can still call it $PE(A)$ as long as we remember what that means now.
Then, the period is twice the time required to go from $A_-$ to $A_+$, and so
$$T=2\intop_{A_-}^{A_+}\frac{dx}{\sqrt{\frac{2}{m}}\sqrt{PE(A)-PE(x)}}$$
$$T=\sqrt{2m}\intop_{A_-}^{A_+}\frac{dx}{\sqrt{PE(A)-PE(x)}}$$
A: From Newton's second law we have (whether $k$ is constant or not) that:
\begin{equation}
m\ddot{x}+kx=0
\end{equation}
The only difference is whether or not $k$ is a function of $t$ or not. If it is a function of $t$, the only general way to solve this differential equation is by using Taylor expansions. Let us take:
\begin{equation}
x\left(t\right)=\sum_{n=0}^\infty a_nt^n
\end{equation}
and:
\begin{equation}
k\left(t\right)=\sum_{n=0}^\infty b_nt^n
\end{equation}
Our differential equation then becomes:
\begin{equation}
\begin{aligned}
m\ddot{x}+kx&=0\\
\implies\sum_{n=2}^\infty mn\left(n-1\right)a_nt^{n-2}+\left(\sum_{n=0}^\infty b_nt^n\right)\left(\sum_{n=0}^\infty a_nt^n\right)&=0\\
\implies\sum_{n=0}^\infty\left[m\left(n+2\right)\left(n+1\right)a_{n+2}+\sum_{i=0}^na_ib_{n-i}\right]t^n&=0\\
\implies m\left(n+2\right)\left(n+1\right)a_{n+2}+\sum_{i=0}^na_ib_{n-i}&=0\forall n\\
\implies a_{n+2}&=-\frac{\sum_{i=0}^na_ib_{n-i}}{m\left(n+2\right)\left(n+1\right)}\forall n
\end{aligned}
\end{equation}
As the $k\left(t\right)$ is known all of the $b_n$ are known, and if we know two of our initial conditions two of the $a_n$ are known (let us say $a_0$ and $a_1$). Using this recurrence relation, one can read off all of the $a_n$--that is, one knows all of the coefficients of the Taylor series for $x$. You can't really see too much more analytically in this super general case (to find a period, one would have to find a $k\left(t\right)$ that generated $a_n$ such that $x\left(t\right)$ was periodic, and read off the period from that function), but a good sanity check is to check if we recover our same answer when $k$ is a constant $k_c$; that is, when $b_0=k_c$ and $b_n=0$ for all $n>0$. In this case we find that:
\begin{equation}
\begin{aligned}
a_2&=-\frac{a_0k_c}{2m}\\
a_3&=-\frac{a_1k_c}{6m}\\
a_4&=\frac{a_0k_c^2}{24m^2}\\
&\vdots
\end{aligned}
\end{equation}
Following the pattern, we notice that the $a_n$ for even $n$ give the Taylor series for $a_0\cos\left(\sqrt{\frac{k_c}{m}}t\right)$ and the $a_n$ for odd $n$ give the Taylor series for $a_1\sin\left(\sqrt{\frac{k_c}{m}}t\right)$, yielding an angular frequency of $\sqrt{\frac{k_c}{m}}$ and therefore a period of $2\pi\sqrt{\frac{m}{k_c}}$.
A: To obtain some kind of practical answer, you have to determine how k varies.  For example, if k varies with temperature, I would determine its value at -50, 0, and 50 degrees, then use those values and calculate T (which varies inversely as the square root of k).  I would Use more points if a higher accuracy is required.  If a formula is required, I would use the "best fit curve" through the points, to generate it. 
A: Alright, according to my knowledge there are some cases with a time dependent spring constant where a closed form solution is known. One of my favorites is the following where the spring constant is a power function. Assume that $k/m = \omega^2/t^\beta$ where $\omega \in {\mathbb R}$ and $\beta \ge 0$. Then the ODE in question take a form:
\begin{equation}
\ddot{x}_t + \frac{\omega^2}{t^\beta} x_t = 0
\end{equation}
and it has the following linearly independent solutions:
\begin{equation}
\sqrt{t} J_{\pm \frac{1}{2 \beta-2}}\left[\frac{\omega}{\beta-1} t^{1-\beta} \right]
\end{equation}
where $J_\cdot[]$ is a Bessel function.
The result is proven in several different ways in https://math.stackexchange.com/questions/673737/solution-to-a-second-order-ordinary-differential-equation .
A quick and dirty "proof" is given by the following Mathematica code:
In[166]:= w =.; b =.; t =.;
Table[FullSimplify[(D[#, {t, 2}] + w^2/t^(2 b) #) & /@ {Sqrt[
      t] BesselJ[eps/(2 b - 2), w/(b - 1) t^(1 - b)]}], {eps, -1, 1, 
  2}]

Out[167]= {{0}, {0}}

Now, having said all this it would be natural to generalize our ODE by by replacing the power function by a linear combination of two power functions. In other words we seek solutions to the following ODE:
\begin{equation}
\ddot{x}_t + \left(\frac{\omega_1^2}{t^\beta_1} + \frac{\omega_2^2}{t^\beta_2}\right)x_t = 0
\end{equation}
I have long been struggling  to solve this problem without any considerable success (however see this https://math.stackexchange.com/questions/1018897/an-ordinary-differential-equation-with-time-varying-coefficients). The usual way "plug the result into Mathematica" or "ask a mathematician" has led to nothing. As it seems mathematics as usual lags several decades behind physics and it will never catch up...;-).
