Oplossingen · hoofdstuk 11Trillingen en resonantie

import math
k, m, zeta = 9839.0, 233.5, 0.05            # mast met een extra demper (model)
w0, c, F0 = math.sqrt(k / m), 2 * zeta * math.sqrt(k * m), 100.0

def evenwichtsamplitude(f, duur=60.0, h=0.004):
    """Heun op m x'' + c x' + k x = F0 cos(w t); de grootste uitwijking in de laatste 10 s."""
    w, x, v, t, top = 2 * math.pi * f, 0.0, 0.0, 0.0, 0.0
    a = lambda t, x, v: (F0 * math.cos(w * t) - c * v - k * x) / m
    while t < duur:
        a1 = a(t, x, v)
        xp, vp = x + h * v, v + h * a1
        a2 = a(t + h, xp, vp)
        x, v, t = x + h / 2 * (v + vp), v + h / 2 * (a1 + a2), t + h
        if t > duur - 10:
            top = max(top, abs(x))
    return top

for f in (0.5, 0.8, 0.95, 1.0, 1.03, 1.06, 1.1, 1.3, 2.0):
    r = 2 * math.pi * f / w0
    formule = F0 / k / math.sqrt((1 - r * r)**2 + (2 * zeta * r)**2)
    print(f"f = {f:4.2f} Hz: simulatie {1000 * evenwichtsamplitude(f):5.1f} mm, formule {1000 * formule:5.1f} mm")

f = 0.50 Hz: simulatie  13.2 mm, formule  13.2 mm
f = 0.80 Hz: simulatie  24.9 mm, formule  24.9 mm
f = 0.95 Hz: simulatie  56.5 mm, formule  56.5 mm
f = 1.00 Hz: simulatie  87.9 mm, formule  88.0 mm
f = 1.03 Hz: simulatie 101.7 mm, formule 101.8 mm
f = 1.06 Hz: simulatie  88.2 mm, formule  88.1 mm
f = 1.10 Hz: simulatie  59.5 mm, formule  59.5 mm
f = 1.30 Hz: simulatie  17.0 mm, formule  17.0 mm
f = 2.00 Hz: simulatie   3.7 mm, formule   3.7 mm

import math
g, L = 9.81, 1.00

def periode(th0_graden, h=1e-4):
    """Heun op th'' = -(g/L) sin th vanuit rust; de tijd tot th voor het eerst door nul gaat is een kwart periode."""
    th, w, t = math.radians(th0_graden), 0.0, 0.0
    while th > 0:
        a1 = -g / L * math.sin(th)
        thp, wp = th + h * w, w + h * a1
        a2 = -g / L * math.sin(thp)
        th_nieuw, w = th + h / 2 * (w + wp), w + h / 2 * (a1 + a2)
        if th_nieuw <= 0:                                      # lineair interpoleren tot het nulpunt
            t += h * th / (th - th_nieuw)
            break
        th, t = th_nieuw, t + h
    return 4 * t

T0 = 2 * math.pi * math.sqrt(L / g)
for a in (5, 10, 20, 30, 45, 60, 90):
    print(f"{a:2d} graden: T/T0 = {periode(a) / T0:.4f}   benadering {1 + math.radians(a)**2 / 16:.4f}")

 5 graden: T/T0 = 1.0005   benadering 1.0005
10 graden: T/T0 = 1.0019   benadering 1.0019
20 graden: T/T0 = 1.0077   benadering 1.0076
30 graden: T/T0 = 1.0174   benadering 1.0171
45 graden: T/T0 = 1.0400   benadering 1.0386
60 graden: T/T0 = 1.0732   benadering 1.0685
90 graden: T/T0 = 1.1803   benadering 1.1542