Oscillator with decaying restoring force Suppose a system that is described by the equation of motion:
$$
\ddot{x} = -k\cdot x\cdot \exp\left(-\frac{t^2}{2\sigma^2}\right).
$$
(For example a spring with decaying stiffness.)
I'd like to obtain the corresponding motion $x(t)$ however I'm not sure how to tackle that problem. I tried WKB approach to first order which gives me:
$$
x(t) = \left[x_0\cos\left(\sqrt{\pi k}\,\sigma\cdot\textrm{erf}(t/2\sigma)\right) + \frac{v_0}{\sqrt{k}}\sin\left(\sqrt{\pi k}\,\sigma\cdot\textrm{erf}(t/2\sigma)\right)\right]\cdot\exp\left(\frac{t^2}{8\sigma^2}\right)
$$
This seems incorrect since for $t\rightarrow\infty$ the restoring force approaches zero and hence I'd expect a linear $t$-dependence in that limit.
Are there any other techniques which are appropriate and can be used to solve this equation of motion?
Edit
As the question has been put on hold, I will clarify my intents in the following. They include an implicit and an explicit question:


*

*(implicit) - Why does the WKB approach yield an incorrect result for that case? Including higher orders doesn't seem to help since the first order already contains an exponential dependence. Are there any criteria for the explicit time dependence which need to be fulfilled in order to make WKB an appropriate approach?

*(explicit) - If WKB doesn't work, what other approach can be used to solve that kind of differential equation? Here I provided an example of a specific explicit time dependence, however, I am also interested in general solutions for arbitrary explicit time dependencies, if possible.


Note This is no homework, but pure private curiosity.
 A: I don't see anything wrong with your solution, and indeed a quick Matlab program confirms this (I used ode45 with $k=1$, $\sigma=1$, $x_{0}=1$ and $v_{0}=0$). Your problem is that you expect too much from your approximation. It is okay, but only on short time scales.

A: It is possible to solve this problem simply by using power series. To simplify notation say $k_t=k*e^{-\frac {t^2}{2\sigma^2}}$. Then assume the solution can be written as:
$$x(t)=\sum_{n=0}^\infty c_n t^n$$
Then, inputing this into the given differential equation:
$$r.h.s.:~~\ddot x(t)=\sum_{n=2}^\infty (n)(n-1)c_n t^n\\l.h.s.:~~-k_t\cdot x(t)=-\sum_{n=0}^\infty k_t \cdot c_n t^n=-\sum_{n=2}^\infty k_t \cdot c_{n-2} t^{n-2}$$
Both being equal you get the recursive relation: $$c_n=\frac {-k_t \cdot c_{n-2}}{(n)(n-1)}$$
So as expected you get two free constants, $c_0$ and $c_1$, and you can divide the series in even and odd terms. Doing that, you get:
$$c_2=\frac {-k_t}{(2)(1)}c_0,~~c_4=\frac {(-k_t)}{(4)(3)}\frac {(-k_t)}{(2)(1)}c_0,... \implies c_{2n}=c_0 \frac {(-k_t)^n}{(2n)!}$$
$$c_3=\frac {(-k_t)}{(3)(2)}c_1,~~c_5=\frac {(-k_t)}{(5)(4)}\frac {(-k_t)}{(3)(2)}c_1,... \implies c_{2n+1}=c_1 \frac {(-k_t)^n}{(2n+1)!}$$
Arranging both even and odd series you get:
$$x(t)=c_0 \sum_{n=0}^\infty\frac{(-1)^n [k_t^{1/2} \cdot t]^{2n}}{(2n)!}+\frac {c_1}{k_t^{1/2}} \sum_{n=0}^\infty\frac{(-1)^n [k_t^{1/2} \cdot t]^{2n+1}}{(2n+1)!}$$
The series are the well known sine and cosine power series expansions, but now with a different parameter, so the result is:
$$x(t)=c_0 \cdot cos(\sqrt k *e^{-\frac {t^2}{4\sigma^2}}\cdot t)+\frac {c_1}{\sqrt k *e^{-\frac {t^2}{4\sigma^2}}} \cdot sin(\sqrt k *e^{-\frac {t^2}{4\sigma^2}}\cdot t)$$
This is an exact solution, no approximation needed, and in this case if you take the limit $t \to \infty$ appropriately you should get a linear time dependence.
A: One way to go about the problem is to construct a series of solutions each of which is valid on a range for which the time-dependent part of the expression can be considered constant.
That means we construct solutions $x_n(t)$ which are valid on $[t_n, t_{n+1}]$ for $\exp(-t_n^2/2\sigma^2) \approx \textrm{const}$. The solutions are of the form
$$
\begin{aligned}
x_n(t) &= a_n\cos(\omega_n t) + b_n\sin(\omega_n t) \\
\dot{x}_n(t) &= -a_n \omega_n\sin(\omega_n t) + b_n\omega_n\cos(\omega_n t) \\
\textrm{with}\;\, \omega_n &= \sqrt{k}\cdot\exp(-t_n^2/2\sigma^2)
\end{aligned}
$$
To connect the solutions at the boundaries $t_n$ we use continuity of $x(t), \dot{x}(t)$:
$$
x_n(t_n) \stackrel{!}{=} x_{n-1}(t_n) \;\;\; \textrm{and} \;\;\; \dot{x}_n(t_n) \stackrel{!}{=} \dot{x}_{n-1}(t_n)
$$
From that follows a recurrence relation for the coefficients $a_n, b_n$ for $n \geq 1$:
$$
\begin{pmatrix} \cos(\omega_n t_n)     & \sin(\omega_n t_n) \\
                -\omega_n\sin(\omega_n t_n) & \omega_n\cos(\omega_n t_n) \end{pmatrix}
\begin{pmatrix} a_n \\ b_n \end{pmatrix}
=
\begin{pmatrix} a_{n-1}\cos(\omega_{n-1} t_n) + b_{n-1}\sin(\omega_{n-1} t_n) \\
                -a_{n-1}\omega_{n-1}\sin(\omega_{n-1} t_n) + b_{n-1}\omega_{n-1}\cos(\omega_{n-1} t_n) \end{pmatrix}
$$
This system of equations has determinant $\omega_n \neq 0$ and hence has a unique solution.
For $n = 0$ we have $a_0 = x_0$ and $b_0 = v_0 / \omega_0$.
Now it remains to choose the values for $t_n$ (with $t_0 = 0$ fixed). Since a prerequisite was the time dependence $Q(t) \equiv \exp(-t^2/2\sigma^2)$ to be approximately constant on the corresponding intervals we can make an estimate via the derivative $\dot{Q}(t)$:
$$
\left\lvert \frac{\dot{Q}(t_n)\cdot (t_{n+1} - t_n)}{Q(t_n)} \right\rvert \leq \epsilon
$$
with a chosen limit $\epsilon$. From this we obtain:
$$
t_{n+1} = \left\lvert\frac{Q(t_n)}{\dot{Q}(t_n)}\right\rvert \cdot \epsilon + t_n
$$
However this only works for $\dot{Q}(t) \neq 0$. In case this is not fulfilled we can just choose a constant update value instead
$$
t_{n+1} = t_n + \delta
$$
for a small value $\delta$; for example:
$$
\delta = \left\lvert\frac{Q(t_{\textrm{max}})}{\dot{Q}(t_{\textrm{max}})}\right\rvert\cdot\epsilon
$$
where $t_{\textrm{max}} = arg\,max \lvert\dot{Q}(t)\rvert$.
Example
Using $x_0 = 1\,\textrm{mm}, v_0 = 0, k = 1\,\textrm{s}^{-2}, \sigma = 2\pi\,\textrm{s}$ and a constant update scheme with $\delta = 0.001\sigma$ we obtain:

Code:
import matplotlib.pyplot as plt
from numpy import cos, sin, sqrt, exp, pi
import numpy as np


def xns(x0, v0, k, s, delta=0.001):
    a = x0
    w = sqrt(k)
    b = v0 / w

    def func(t):
        return a*cos(w*t) + b*sin(w*t)

    t = delta * s

    yield t, func

    while True:
        v = w
        w = sqrt(k) * exp(-t**2 / (2 * s**2))
        a, b = np.linalg.solve([
            [cos(w*t), sin(w*t)],
            [-w*sin(w*t), w*cos(w*t)]],

            [a*cos(v*t) + b*sin(v*t),
            -a*v*sin(v*t) + b*v*cos(v*t)]
        )
        t += delta*s
        yield t, func


x0, v0, k, s = 0.001, 0, 1, 2*pi
ts = np.linspace(0, 7*pi*sqrt(k), 1000)
xns = xns(x0, v0, k, s)

xs = []
tn, xn = next(xns)
for t in ts:
    while t >= tn:
        tn, xn = next(xns)
    xs.append(xn(t))
xs = np.asarray(xs)

plt.figure()
plt.title(r'$\rm x_0 = {:.2f}\,[mm],\;'
      r'     v_0 = {:.2f}\,[mm\,s^{{-1}}],\;'
      r'     k   = {:.2f}\,[s^{{-2}}],\;'
      r'\sigma   = {:.2f}\,[s]$'
          .format(x0*1e3, v0*1e3, k, s))
plt.xlabel('t [s]')
plt.ylabel('x [mm]')
plt.plot(ts, xs * 1e3, '-', lw=3)
plt.grid(lw=1.5)
plt.twinx()
plt.ylabel(r'$\rm\sqrt{k}\cdot\exp\left(-t^2/2\sigma^2\right)$', color='#ff7f0e')
plt.tick_params('y', colors='#ff7f0e')
plt.plot(ts, sqrt(k) * exp(-ts**2 / (2 * s**2)), '--', color='#ff7f0e', lw=2)
plt.savefig('example.png', bbox_inches='tight', pad_zeros=0)
plt.show()

