1.5 Numerisk modellering av bevegelse

Hvorfor trenger vi numeriske metoder?

I de tidligere kapitlene har vi løst bevegelseslikninger analytisk – det vil si at vi har funnet eksakte formler for posisjon og hastighet som funksjon av tid. For eksempel lærte vi at for konstant akselerasjon gjelder:

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

Men i virkeligheten er akselerasjonen sjelden konstant. Tenk på:
- Luftmotstand som øker med hastigheten
- Fjærkrefter som avhenger av posisjon
- Gravitasjon som varierer med avstand
- Magnetiske krefter som avhenger av både posisjon og hastighet

For slike problemer finnes det ofte ingen analytisk løsning. Da bruker vi i stedet numeriske metoder – algoritmer som beregner tilnærmede verdier steg for steg.

Numeriske metoder:
- Deler opp tiden i små intervaller $\Delta t$
- Beregner posisjon og hastighet trinnvis
- Gir tilnærmede, men ofte svært nøyaktige, løsninger
- Kan implementeres på datamaskiner

Anvendelser:
- Simulering av romfart og satellittbaner
- Værvarsling og klimamodeller
- Dataspill og animasjoner
- Partikkelfysikk og molekylsimuleringer
- Strukturanalyse i ingeniørfag

Eulers metode – matematisk bakgrunn

Eulers metode (oppkalt etter matematikeren Leonhard Euler, 1707-1783) er den enkleste numeriske metoden for å løse differensiallikninger.

Den deriverte som stigningstall

Husk at den deriverte av en funksjon $f(t)$ er definert som:

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

Dette er stigningstallet til tangenten ved punktet $t$.

Eulers approksimasjon

Hvis vi ikke tar grensen, men bruker en liten, men endelig $\Delta t$, får vi en tilnærming:

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

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

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

Geometrisk tolkning

Euler-metoden følger tangentlinjen i hvert punkt. I stedet for å følge den buede kurven eksakt, tar vi et lite skritt langs tangenten. Jo mindre $\Delta t$, jo bedre følger vi den virkelige kurven.

Viktig innsikt:
Euler-metoden erstatter den buede løsningskurven med en rekke rette linjesegmenter (tangentene). Dette gir en tilnærmet løsning som blir bedre jo mindre tidssteg vi bruker.

Implementasjon av Eulers metode

Algoritmen steg for steg

1. Sett startverdier: Velg $x_0$, $v_0$, og $t_0 = 0$
2. Velg tidsstegg: Bestem $\Delta t$ (mindre = mer nøyaktig, men flere beregninger)
3. Gjenta for hvert tidssteg:
- Beregn akselerasjonen $a_n = F_n / m$ (fra kreftene)
- Oppdater hastigheten: $v_{n+1} = v_n + a_n \cdot \Delta t$
- Oppdater posisjonen: $x_{n+1} = x_n + v_n \cdot \Delta t$
- Oppdater tiden: $t_{n+1} = t_n + \Delta t$
4. Fortsett til ønsket 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

Valg av tidsstegg

Tommelregel: $\Delta t$ bør være mye mindre enn den karakteristiske tidsskalaen for systemet.

- Fritt fall (noen sekunder): $\Delta t \approx 0.01 \text{ s}$ - $0.1 \text{ s}$
- Pendel (perioder på 1-2 s): $\Delta t \approx 0.001 \text{ s}$ - $0.01 \text{ s}$
- Raske svingninger: Bruk enda mindre $\Delta t$

Feilkilder og nøyaktighet

Diskretiseringsfeil (avkuttingsfeil)

Den største feilkilden i Euler-metoden er diskretiseringsfeil. Denne oppstår fordi vi erstatter den glatte løsningskurven med rette linjesegmenter.

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

Dette betyr at hvis vi halverer tidsstegget, halverer vi den totale feilen. Euler-metoden er derfor en førsteordens metode.

Akkumulering av feil

For hvert steg legges det til en liten feil. Over mange steg kan disse feilene akkumuleres og vokse:
- Feil i posisjon påvirker beregning av krefter
- Feil i hastighet påvirker beregning av posisjon
- Feilene kan forsterke hverandre

Numerisk stabilitet

For noen problemer kan Euler-metoden bli ustabil – feilene vokser eksponentielt. Dette skjer typisk når:
- Tidsstegget er for stort
- Systemet har raske svingninger
- Kreftene er sterkt avhengige av posisjon

Stabilitetsregel: $\Delta t$ må være mye mindre enn systemets karakteristiske tidsperiode.

Forbedrede metoder

Det finnes mange bedre metoder enn enkel Euler:
- Euler-Cromer (symplektisk Euler): Bruker $v_{n+1}$ i stedet for $v_n$ i posisjonsoppdateringen
- Midtpunktmetoden: Bruker hastigheten midt i intervallet
- Runge-Kutta: Mer nøyaktige, men mer komplekse metoder

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

Denne enkle endringen gir mye bedre energibevaring for periodiske systemer.

Python-eksempel: Fritt fall

Her er et 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")

Forventet resultat:
- Ballen treffer bakken etter ca. 4.52 sekunder
- Sluthastighet: 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-eksempel: 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 observasjoner:
1. Hastigheten øker raskt i starten
2. Økningen avtar når luftmotstanden blir betydelig
3. Hastigheten nærmer seg terminalhastigheten asymptotisk
4. Den numeriske løsningen matcher den analytiske svært godt

Python-eksempel: Matematisk pendel

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

For små vinkler ($\theta < 15°$) kan vi bruke $\sin\theta \approx \theta$ og få en analytisk løsning. Men for store vinkler må vi løse 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 observasjoner:
1. Euler-Cromer bevarer energien godt (viktig for oscillerende systemer!)
2. For store vinkler er perioden litt lengre enn liten-vinkel-formelen gir
3. Vanlig Euler-metode ville gitt økende energi over tid (ustabilt)

Anvendelser av numerisk modellering

1. Skrått kast med luftmotstand

Uten luftmotstand har vi den kjente parabelbanen. Med luftmotstand:
- Banen blir asymmetrisk (brattere nedstigning)
- Rekkevidden reduseres betydelig
- Optimal kastevinkel er under 45°

2. Fallskjermhopping

Simulering med varierende luftmotstandskoeffisient:
- Fritt fall: lav $b$, hastighet ca. 200 km/t
- Fallskjerm utløst: høy $b$, hastighet ca. 20 km/t
- Brå endring i $b$ gir rask oppbremsing

3. Planetbaner

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

Må løses numerisk fordi akselerasjonen avhenger av posisjon.

4. Bilkjøring og bremsing

- Bremskraft avhengig av hastighet (ABS-bremser)
- Friksjon som varierer med underlaget
- Luftmotstand ved høye hastigheter

5. Sportsanalyse

- Golfball med spin (Magnus-effekt)
- Fotball som svinger i lufta
- Skihopperes luftmotstand og oppdrift

6. Værmodellering

- Atmosfæriske bevegelser
- Havsirkulasjon
- Klimaendringer over tid

Oppsummering

Viktige konsepter:

1. Numeriske metoder brukes når analytiske løsninger ikke eksisterer eller er for komplekse

2. Eulers metode er den enkleste 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 en forbedring som bruker ny hastighet i posisjonsoppdateringen:
- $v_{n+1} = v_n + a_n \cdot \Delta t$
- $x_{n+1} = x_n + v_{n+1} \cdot \Delta t$

4. Tidsstegg $\Delta t$: Mindre steg gir bedre nøyaktighet, men krever flere beregninger

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

6. Anvendelser:
- Bevegelse med luftmotstand
- Pendler og svingninger
- Planetbaner
- Komplekse systemer med variable krefter

Programmeringstips:
- Bruk Euler-Cromer for oscillerende systemer (bevarer energi bedre)
- Velg $\Delta t$ mye mindre enn systemets karakteristiske tidsskala
- Kontroller resultater mot analytiske løsninger når mulig
- Sjekk energibevaring for å verifisere stabilitet

Neste kapittel:
Vi skal se på gravitasjon og planetbevegelse, der numeriske metoder er essensielle for å beregne satellittbaner og interplanetære reiser.