1.5 Numerisk modellering av bevegelse

Kvifor treng vi numeriske metodar?

I dei tidlegare kapitla har vi løyst rørslelikningar analytisk – det vil seie at vi har funne eksakte formlar for posisjon og fart som funksjon av tid. Til dømes lærte vi at for konstant akselerasjon gjeld:

$$v(t) = v_0 + at$$
$$x(t) = x_0 + v_0 t + \frac{1}{2}at^2$$

Men i røynda er akselerasjonen sjeldan konstant. Tenk på:
- Luftmotstand som aukar med farten
- Fjørkrefter som avheng av posisjon
- Gravitasjon som varierer med avstand
- Magnetiske krefter som avheng av både posisjon og fart

For slike problem finst det ofte inga analytisk løysing. Då bruker vi i staden numeriske metodar – algoritmar som reknar ut tilnærma verdiar steg for steg.

Numeriske metodar:
- Deler opp tida i små intervall $\Delta t$
- Reknar ut posisjon og fart trinnvis
- Gjev tilnærma, men ofte svært nøyaktige, løysingar
- Kan implementerast på datamaskiner

Bruksområde:
- Simulering av romfart og satellittbaner
- Vêrvarsling og klimamodellar
- Dataspel og animasjonar
- Partikkelfysikk og molekylsimuleringar
- Strukturanalyse i ingeniørfag

Eulers metode – matematisk bakgrunn

Eulers metode (oppkalla etter matematikaren Leonhard Euler, 1707-1783) er den enklaste numeriske metoden for å løyse differensiallikningar.

Den deriverte som stigingstal

Hugs at den deriverte av ein funksjon $f(t)$ er definert som:

$$f'(t) = \lim_{\Delta t \to 0} \frac{f(t + \Delta t) - f(t)}{\Delta t}$$

Dette er stigingstalet til tangenten ved punktet $t$.

Eulers approksimasjon

Viss vi ikkje tek grensa, men bruker ein liten, men endeleg $\Delta t$, får vi ei tilnærming:

$$f'(t) \approx \frac{f(t + \Delta t) - f(t)}{\Delta t}$$

Ved å løyse for $f(t + \Delta t)$ får vi Euler-formelen:

$$f(t + \Delta t) \approx f(t) + f'(t) \cdot \Delta t$$

Geometrisk tolking

Euler-metoden følgjer tangentlinja i kvart punkt. I staden for å følgje den bøygde kurva eksakt, tek vi eit lite steg langs tangenten. Jo mindre $\Delta t$, jo betre følgjer vi den verkelege kurva.

Viktig innsikt:
Euler-metoden erstattar den bøygde løysingskurva med ei rekkje rette linjesegment (tangentane). Dette gjev ei tilnærma løysing som blir betre jo mindre tidssteg vi bruker.

Implementasjon av Eulers metode

Algoritmen steg for steg

1. Set startverdiar: Vel $x_0$, $v_0$, og $t_0 = 0$
2. Vel tidssteg: Fastsett $\Delta t$ (mindre = meir nøyaktig, men fleire utrekningar)
3. Gjenta for kvart tidssteg:
- Rekn ut akselerasjonen $a_n = F_n / m$ (frå kreftene)
- Oppdater farten: $v_{n+1} = v_n + a_n \cdot \Delta t$
- Oppdater posisjonen: $x_{n+1} = x_n + v_n \cdot \Delta t$
- Oppdater tida: $t_{n+1} = t_n + \Delta t$
4. Hald fram til ønskt sluttid er nådd

Pseudokode

x = x0          # Startposisjon
v = v0          # Starthastighet
t = 0           # Starttid
dt = 0.01       # Tidsstegg

while t < t_slutt:
    a = F(x, v, t) / m    # Beregn akselerasjon fra krefter
    v = v + a * dt        # Oppdater hastighet
    x = x + v * dt        # Oppdater posisjon
    t = t + dt            # Oppdater tid

Val av tidssteg

Tommelregel: $\Delta t$ bør vere mykje mindre enn den karakteristiske tidsskalaen for systemet.

- Fritt fall (nokre sekund): $\Delta t \approx 0.01 \text{ s}$ - $0.1 \text{ s}$
- Pendel (periodar på 1-2 s): $\Delta t \approx 0.001 \text{ s}$ - $0.01 \text{ s}$
- Raske svingingar: Bruk endå mindre $\Delta t$

Feilkjelder og nøyaktigheit

Diskretiseringsfeil (avkuttingsfeil)

Den største feilkjelda i Euler-metoden er diskretiseringsfeil. Denne oppstår fordi vi erstattar den glatte løysingskurva med rette linjesegment.

For Eulers metode gjeld:
- Feilen per steg er proporsjonal med $(\Delta t)^2$
- Total feil etter mange steg er proporsjonal med $\Delta t$

Dette tyder at viss vi halverer tidssteget, halverer vi den totale feilen. Euler-metoden er derfor ein førsteordens metode.

Akkumulering av feil

For kvart steg blir det lagt til ein liten feil. Over mange steg kan desse feila akkumulerast og vekse:
- Feil i posisjon påverkar utrekning av krefter
- Feil i fart påverkar utrekning av posisjon
- Feila kan forsterke kvarandre

Numerisk stabilitet

For nokre problem kan Euler-metoden bli ustabil – feila veks eksponentielt. Dette skjer typisk når:
- Tidssteget er for stort
- Systemet har raske svingingar
- Kreftene er sterkt avhengige av posisjon

Stabilitetsregel: $\Delta t$ må vere mykje mindre enn den karakteristiske tidsperioden til systemet.

Forbetra metodar

Det finst mange betre metodar enn enkel Euler:
- Euler-Cromer (symplektisk Euler): Bruker $v_{n+1}$ i staden for $v_n$ i posisjonsoppdateringa
- Midtpunktmetoden: Bruker farten midt i intervallet
- Runge-Kutta: Meir nøyaktige, men meir komplekse metodar

Euler-Cromer-modifikasjon:
$$v_{n+1} = v_n + a_n \cdot \Delta t$$
$$x_{n+1} = x_n + v_{n+1} \cdot \Delta t \quad \text{(bruker ny fart!)}$$

Denne enkle endringa gjev mykje betre energibevaring for periodiske system.

Python-døme: Fritt fall

Her er eit komplett Python-program som simulerer fritt fall med Euler-Cromer-metoden:

# Fritt fall med Euler-Cromer-metoden
import matplotlib.pyplot as plt

# Fysiske konstanter
g = 9.81  # Tyngdeakselerasjon (m/s^2)

# Startverdier
y = 100   # Starthøyde (m)
v = 0     # Starthastighet (m/s)
t = 0     # Starttid (s)

# Numeriske parametere
dt = 0.01  # Tidsstegg (s)
t_slutt = 5  # Sluttid (s)

# Lister for å lagre resultater
t_liste = [t]
y_liste = [y]
v_liste = [v]

# Hovedløkke - Euler-Cromer
while t < t_slutt and y > 0:  # Stopp når ballen treffer bakken
    # Akselerasjon (konstant for fritt fall)
    a = -g

    # Oppdater hastighet først
    v = v + a * dt

    # Deretter posisjon (med ny hastighet - Euler-Cromer)
    y = y + v * dt

    # Oppdater tid
    t = t + dt

    # Lagre verdier
    t_liste.append(t)
    y_liste.append(y)
    v_liste.append(v)

# Plot resultatene
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 8))

ax1.plot(t_liste, y_liste, 'b-', label='Numerisk')
ax1.set_xlabel('Tid (s)')
ax1.set_ylabel('Høyde (m)')
ax1.set_title('Fritt fall - Posisjon')
ax1.grid(True)
ax1.legend()

ax2.plot(t_liste, v_liste, 'r-', label='Numerisk')
ax2.set_xlabel('Tid (s)')
ax2.set_ylabel('Hastighet (m/s)')
ax2.set_title('Fritt fall - Hastighet')
ax2.grid(True)
ax2.legend()

plt.tight_layout()
plt.show()

# Skriv ut sluttresultat
print(f"Ballen traff bakken etter {t:.2f} sekunder")
print(f"Sluthastighet: {v:.2f} m/s")

Forventa resultat:
- Ballen treff bakken etter ca. 4.52 sekund
- Sluttfart: ca. -44.3 m/s

Analytisk kontroll:
$$t = \sqrt{\frac{2h}{g}} = \sqrt{\frac{2 \cdot 100}{9.81}} = 4.52 \text{ s}$$
$$v = gt = 9.81 \cdot 4.52 = 44.3 \text{ m/s}$$

Python-døme: Fall med luftmotstand

# Fall med luftmotstand - kvadratisk drag
import matplotlib.pyplot as plt
import numpy as np

# Fysiske parametere
g = 9.81      # Tyngdeakselerasjon (m/s^2)
m = 80        # Masse (kg)
b = 0.25      # Luftmotstandskoeffisient (kg/m)

# Terminalhastighet (for sammenligning)
v_terminal = np.sqrt(m * g / b)
print(f"Terminalhastighet: {v_terminal:.1f} m/s")

# Startverdier
y = 4000  # Starthøyde (m)
v = 0     # Starthastighet (m/s)
t = 0     # Starttid (s)

# Numeriske parametere
dt = 0.01     # Tidsstegg (s)
t_slutt = 60  # Sluttid (s)

# Lister for resultater
t_liste = [t]
y_liste = [y]
v_liste = [v]

# Hovedløkke - Euler-Cromer
while t < t_slutt and y > 0:
    # Akselerasjon med luftmotstand
    # F = mg - b*v^2 (positiv v betyr nedover)
    a = g - (b/m) * v * abs(v)

    # Euler-Cromer oppdatering
    v = v + a * dt
    y = y - v * dt  # Minus fordi y er høyde, v er positiv nedover
    t = t + dt

    # Lagre
    t_liste.append(t)
    y_liste.append(y)
    v_liste.append(v)

# Analytisk løsning for hastighet
t_analytisk = np.linspace(0, max(t_liste), 100)
v_analytisk = v_terminal * np.tanh(g * t_analytisk / v_terminal)

# Plot
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 8))

ax1.plot(t_liste, y_liste, 'b-')
ax1.set_xlabel('Tid (s)')
ax1.set_ylabel('Høyde (m)')
ax1.set_title('Fall med luftmotstand - Posisjon')
ax1.grid(True)

ax2.plot(t_liste, v_liste, 'b-', label='Numerisk')
ax2.plot(t_analytisk, v_analytisk, 'r--', label='Analytisk')
ax2.axhline(y=v_terminal, color='g', linestyle=':',
            label=f'Terminal: {v_terminal:.1f} m/s')
ax2.set_xlabel('Tid (s)')
ax2.set_ylabel('Hastighet (m/s)')
ax2.set_title('Fall med luftmotstand - Hastighet')
ax2.legend()
ax2.grid(True)

plt.tight_layout()
plt.show()

Viktige observasjonar:
1. Farten aukar raskt i starten
2. Auken avtek når luftmotstanden blir monaleg
3. Farten nærmar seg terminalfarten asymptotisk
4. Den numeriske løysinga matchar den analytiske svært godt

Python-døme: Matematisk pendel

For ein pendel med lengd $L$ og utslag $\theta$ er akselerasjonen:
$$\alpha = -\frac{g}{L}\sin\theta$$

For små vinklar ($\theta < 15°$) kan vi bruke $\sin\theta \approx \theta$ og få ei analytisk løysing. Men for store vinklar må vi løyse numerisk.

# Matematisk pendel - numerisk løsning
import matplotlib.pyplot as plt
import numpy as np

# Fysiske parametere
g = 9.81   # Tyngdeakselerasjon (m/s^2)
L = 1.0    # Pendellengde (m)

# Startverdier
theta = np.radians(45)  # Startvinkel (45 grader)
omega = 0               # Startvinkelhastighet (rad/s)
t = 0

# Numeriske parametere
dt = 0.001    # Tidsstegg (s) - må være lite for pendel!
t_slutt = 10  # Sluttid (s)

# Lister for resultater
t_liste = [t]
theta_liste = [theta]
omega_liste = [omega]

# Beregn total mekanisk energi (for å sjekke bevaring)
def energi(theta, omega):
    E_kin = 0.5 * L**2 * omega**2  # Kinetisk (per masse)
    E_pot = g * L * (1 - np.cos(theta))  # Potensiell (per masse)
    return E_kin + E_pot

E0 = energi(theta, omega)
E_liste = [E0]

# Hovedløkke - Euler-Cromer
while t < t_slutt:
    # Vinkelakselerasjon
    alpha = -(g/L) * np.sin(theta)

    # Euler-Cromer
    omega = omega + alpha * dt
    theta = theta + omega * dt
    t = t + dt

    # Lagre
    t_liste.append(t)
    theta_liste.append(theta)
    omega_liste.append(omega)
    E_liste.append(energi(theta, omega))

# Plot
fig, (ax1, ax2, ax3) = plt.subplots(3, 1, figsize=(10, 10))

ax1.plot(t_liste, np.degrees(theta_liste))
ax1.set_xlabel('Tid (s)')
ax1.set_ylabel('Vinkel (grader)')
ax1.set_title('Pendel - Vinkelutslag')
ax1.grid(True)

ax2.plot(t_liste, omega_liste)
ax2.set_xlabel('Tid (s)')
ax2.set_ylabel('Vinkelhastighet (rad/s)')
ax2.set_title('Pendel - Vinkelhastighet')
ax2.grid(True)

# Energibevaring
ax3.plot(t_liste, np.array(E_liste)/E0)
ax3.set_xlabel('Tid (s)')
ax3.set_ylabel('E/E0')
ax3.set_title('Energibevaring (bør være 1.0)')
ax3.set_ylim(0.99, 1.01)
ax3.grid(True)

plt.tight_layout()
plt.show()

# Beregn periode
perioder = []
forrige_null = None
for i in range(1, len(theta_liste)):
    if theta_liste[i-1] > 0 and theta_liste[i] <= 0:
        if forrige_null is not None:
            perioder.append(t_liste[i] - forrige_null)
        forrige_null = t_liste[i]

if perioder:
    print(f"Målt periode: {np.mean(perioder):.3f} s")
    print(f"Liten-vinkel-formel: {2*np.pi*np.sqrt(L/g):.3f} s")

Viktige observasjonar:
1. Euler-Cromer bevarer energien godt (viktig for oscillerande system!)
2. For store vinklar er perioden litt lengre enn liten-vinkel-formelen gjev
3. Vanleg Euler-metode ville gitt aukande energi over tid (ustabilt)

Bruksområde for numerisk modellering

1. Skrått kast med luftmotstand

Utan luftmotstand har vi den kjende parabelbanen. Med luftmotstand:
- Banen blir asymmetrisk (brattare nedstiging)
- Rekkjevidda blir monaleg redusert
- Optimal kastevinkel er under 45°

2. Fallskjermhopping

Simulering med varierande luftmotstandskoeffisient:
- Fritt fall: låg $b$, fart ca. 200 km/t
- Fallskjerm utløyst: høg $b$, fart ca. 20 km/t
- Brå endring i $b$ gjev rask oppbremsing

3. Planetbaner

For elliptiske baner der gravitasjonen varierer med avstand:
$$F = -\frac{GMm}{r^2}$$

Må løysast numerisk fordi akselerasjonen avheng av posisjon.

4. Bilkøyring og bremsing

- Bremsekraft avhengig av fart (ABS-bremser)
- Friksjon som varierer med underlaget
- Luftmotstand ved høge fartar

5. Sportsanalyse

- Golfball med spinn (Magnus-effekt)
- Fotball som svingar i lufta
- Luftmotstand og oppdrift til skihopparar

6. Vêrmodellering

- Atmosfæriske rørsler
- Havsirkulasjon
- Klimaendringar over tid

Oppsummering

Viktige konsept:

1. Numeriske metodar blir brukte når analytiske løysingar ikkje finst eller er for komplekse

2. Eulers metode er den enklaste numeriske metoden:
- $v_{n+1} = v_n + a_n \cdot \Delta t$
- $x_{n+1} = x_n + v_n \cdot \Delta t$

3. Euler-Cromer-metoden er ei forbetring som bruker ny fart i posisjonsoppdateringa:
- $v_{n+1} = v_n + a_n \cdot \Delta t$
- $x_{n+1} = x_n + v_{n+1} \cdot \Delta t$

4. Tidssteg $\Delta t$: Mindre steg gjev betre nøyaktigheit, men krev fleire utrekningar

5. Feilkjelder:
- Diskretiseringsfeil (avkuttingsfeil)
- Akkumulering av feil over tid
- Numerisk ustabilitet

6. Bruksområde:
- Rørsle med luftmotstand
- Pendlar og svingingar
- Planetbaner
- Komplekse system med variable krefter

Programmeringstips:
- Bruk Euler-Cromer for oscillerande system (bevarer energi betre)
- Vel $\Delta t$ mykje mindre enn den karakteristiske tidsskalaen til systemet
- Kontroller resultat mot analytiske løysingar når det er mogleg
- Sjekk energibevaring for å verifisere stabilitet

Neste kapittel:
Vi skal sjå på gravitasjon og planetrørsle, der numeriske metodar er essensielle for å rekne ut satellittbaner og interplanetære reiser.