# 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.

## 3 Answers

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.

• Thanks for your answer. Generally I would like to obtain a solution for $t \geq 0$. Why would you argue that the WKB is only valid for small $t$? Also in the plot you show, it appears that the ode45 increases even stronger than the WKB, i.e. an even stronger dependence than exponential (while it should be linear for $t\rightarrow\infty$). For me it is not clear how this plot supports the correctness of ode45 results and neither how it shows that WKB is only valid on small time scales (btw text reads $x_0=1$ however plot shows $x_0=-1$). Could you please provide more details on that? Thank you. Commented Nov 1, 2018 at 14:15

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.

• Thanks for your answer. Some things remain unclear however. At first you assume that the solution can be written as a power series in $t$ with coefficients $c_n$. Expressing the time dependence as $t^n$ implies $c_n = \textrm{const}$ (especially when applying the derivative) however later on you derive $c_n \propto (-k_t)^n$ which is not constant but depends on $t$. That contradicts the initial assumption. Also from your final representation of $x(t)$ follows that $c_0 = x_0, c_1 = v_0$ and for $v_0 = 0$ the linear time dependence for $t \rightarrow \infty$ vanishes which doesn't make sense. Commented Nov 1, 2018 at 14:07

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()