# Example on interactive live documents versus traditional static documents **hpl** Date: **May 25, 2015** ## Physics

We consider a vibrating mechanical system as shown in [Figure](#ipynbex:physics:fig).

Oscillating system with spring and damper.

Mathematics

A differential equation for the system in Figure reads

$$ mu'' + f(u') + s(u) = F(t),\quad u(0)=I,\ u'(0)=V\thinspace . $$

In the linear damping case, where $f(u')=bu'$ for some constant $b\geq 0$, we can solve the problem numerically by a scheme

$$ \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}, \label{ipynbex:u:scheme:lin} \tag{1} \end{equation} $$

for $n=0,1,\ldots$. A special formula is required for $n=0$:

$$ \begin{equation} u^1 = u^0 + \Delta t\, V + \frac{\Delta t^2}{2m}(-bV - s(u^0) + F^0) \thinspace . \label{ipynbex:u:scheme0:lin} \tag{2} \end{equation} $$

Implementation

The formulas (1) and (2) can be implemented as follows in a Python function:



In [1]:

    
from numpy import zeros



In [2]:

    
def solver_linear_damping(I, V, m, b, s, F, t):
    N = t.size - 1              # No of time intervals
    dt = t[1] - t[0]            # Time step
    u = zeros(N+1)              # Result array
    u[0] = I
    u[1] = u[0] + dt*V + dt**2/(2*m)*(-b*V - s(u[0]) + F[0])

    for n in range(1,N):
        u[n+1] = 1./(m + b*dt/2)*(2*m*u[n] + \
                 (b*dt/2 - m)*u[n-1] + dt**2*(F[n] - s(u[n])))
    return u

Dissection

The array t holds all the time points where we want a solution. The total number of intervals, $N$, is then computed as

N = t.size - 1   # or len(t) - 1

Creating an array of length $N+1$ where we can store the solution is done by

u = zeros(N+1)

For loops over array indices are coded as

for n in range(1, N):
            ...

which generates a sequence of integers from 1 up to N, but not including N.

Usage

The function solve_linear_damping resides in a file solver.py.



In [3]:

    
# The solver module is in src/solver.py; tell Python about that
import sys
sys.path.insert(0, 'src')



In [4]:

    
%matplotlib inline

from solver import solver_linear_damping
from numpy import linspace, zeros, pi

def s(u):
    return 2*u

T = 10*pi      # simulate for t in [0,T]
dt = 0.2
N = int(round(T/dt))
t = linspace(0, T, N+1)
F = zeros(t.size)
I = 1; V = 0
m = 2; b = 0.2
u = solver_linear_damping(I, V, m, b, s, F, t)

from matplotlib.pyplot import *
plot(t, u)
show()