The differential equation problem (1) can be solved by introducing finite difference approximations for \( u'' \) and \( u' \). In case of quadratic damping one can use a geometric mean to approximate \( u'|u'| \) and thereby linearize the equations. The result of using such numerical methods is an algorithm for computing \( u(t) \) at discrete points in time. Let \( u^n \) be the approximation to \( u \) at time \( t_n=n\Delta t \), \( n=1,2,\ldots \), where \( \Delta t \) is a (small) time interval. For example, if \( \Delta t = 0.1 \), we find approximations \( u^1 \) to \( u \) at \( t=0.1 \), \( u^2 \) at \( t=0.2 \), \( u^3 \) to \( t=0.3 \), and so forth. Any value \( u^{n+1} \) can be computed if \( u^n \) and \( u^{n-1} \) are known (i.e., previously computed). The formula for \( u^{n+1} \) is, in case of linear damping \( f(u')=bu' \), $$ \begin{equation} u^{n+1} = \left(2mu^n + (\frac{b}{2}\Delta t - m)u^{n-1} + \Delta t^2(F^n - s(u^n)) \right)(m + \frac{b}{2}\Delta t)^{-1}, \tag{2} \end{equation} $$ where \( F^n \) means \( F(t) \) evaluated for \( t=t_n \). A special formula must be applied for \( n=0 \): $$ \begin{equation} u^1 = u^0 + \Delta t\, V + \frac{\Delta t^2}{2m}(-bV - s(u^0) + F^0) \thinspace . \tag{3} \end{equation} $$ For quadratic damping we have a slightly different formula, $$ \begin{align} u^{n+1} = & \left( m + b|u^n-u^{n-1}|\right)^{-1}\times\nonumber\\ &\quad \left(2m u^n - mu^{n-1} + bu^n|u^n-u^{n-1}| + \Delta t^2 (F^n - s(u^n)) \right), \tag{4} \end{align} $$ and again a special formula for \( u^1 \): $$ \begin{equation} u^1 = u^0 + \Delta t V + \frac{\Delta t^2}{2m}\left(-bV|V| - s(u^0) + F^0\right) \thinspace . \tag{5} \end{equation} $$
The implementation of the computational algorithm can make use
of an array u to represent \( u^n \) as u[n]. The force \( F(t_n) \) is
assumed to be available as an array element F[n]. The following
Python function computes u given an array t
with time points \( t_0,t_1,\ldots \),
the initial displacement I, mass m, damping parameter b,
restoring force s(u), environmental forces F as an array
(corresponding to t).