Я не вижу ничего плохого в вашем решении, и действительно, быстрая программа Matlab подтверждает это (я использовал ode45 с$k=1$,$\sigma=1$,$x_{0}=1$и$v_{0}=0$). Ваша проблема в том, что вы слишком многого ожидаете от своего приближения. Это нормально, но только в короткие сроки.

Эту проблему можно решить, просто используя степенные ряды. Чтобы упростить запись, скажите$k_t=k*e^{-\frac {t^2}{2\sigma^2}}$. Тогда предположим, что решение может быть записано как:

$$x(t)=\sum_{n=0}^\infty c_n t^n$$

Затем, вводя это в данное дифференциальное уравнение:

$$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}$$ Оба равны, вы получаете рекурсивное отношение: $$c_n=\frac {-k_t \cdot c_{n-2}}{(n)(n-1)}$$

Итак, как и ожидалось, вы получаете две бесплатные константы,$c_0$и$c_1$, и вы можете разделить ряд на четные и нечетные члены. Делая это, вы получаете:

$$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)!}$$

Выстраивая как четные, так и нечетные ряды, вы получаете:

$$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)!}$$

Ряды представляют собой хорошо известные разложения по степеням синуса и косинуса, но теперь с другим параметром, поэтому результат:

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

Это точное решение, приближение не требуется, и в этом случае, если взять предел$t \to \infty$соответственно вы должны получить линейную зависимость от времени.

Один из способов решения проблемы — построить ряд решений, каждое из которых справедливо в диапазоне, для которого зависящую от времени часть выражения можно считать постоянной.

Это означает, что мы строим решения$x_n(t)$которые действительны на$[t_n, t_{n+1}]$для$\exp(-t_n^2/2\sigma^2) \approx \textrm{const}$. Решения имеют вид

$$\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}$$

Чтобы соединить решения на границах$t_n$мы используем непрерывность$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)$$

Отсюда следует рекуррентное соотношение для коэффициентов$a_n, b_n$для$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}$$

Эта система уравнений имеет определитель$\omega_n \neq 0$и, следовательно, имеет единственное решение. Для$n = 0$у нас есть$a_0 = x_0$и$b_0 = v_0 / \omega_0$.

Теперь осталось выбрать значения для$t_n$(с$t_0 = 0$зафиксированный). Поскольку обязательным условием была временная зависимость$Q(t) \equiv \exp(-t^2/2\sigma^2)$быть приблизительно постоянной на соответствующих интервалах, мы можем сделать оценку через производную$\dot{Q}(t)$:

$$\left\lvert \frac{\dot{Q}(t_n)\cdot (t_{n+1} - t_n)}{Q(t_n)} \right\rvert \leq \epsilon$$

с выбранным лимитом$\epsilon$. Отсюда получаем:

$$t_{n+1} = \left\lvert\frac{Q(t_n)}{\dot{Q}(t_n)}\right\rvert \cdot \epsilon + t_n$$

Однако это работает только для$\dot{Q}(t) \neq 0$. В случае, если это не выполнено, мы можем просто выбрать постоянное значение обновления вместо этого.

$$t_{n+1} = t_n + \delta$$

за небольшую стоимость$\delta$; например:

$$\delta = \left\lvert\frac{Q(t_{\textrm{max}})}{\dot{Q}(t_{\textrm{max}})}\right\rvert\cdot\epsilon$$

где$t_{\textrm{max}} = arg\,max \lvert\dot{Q}(t)\rvert$.

Пример

С использованием$x_0 = 1\,\textrm{mm}, v_0 = 0, k = 1\,\textrm{s}^{-2}, \sigma = 2\pi\,\textrm{s}$и схема постоянного обновления с$\delta = 0.001\sigma$мы получаем:

Код:

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