1.3 Numeriske metoder og programmering

Numeriske metoder og programmering

Hvorfor trenger vi numeriske metoder?

Mange fysiske problemer kan ikke løses analytisk (med penn og papir). Eksempler:
- Fritt fall med luftmotstand
- Trekroppsproblemer i astronomi
- Værsimuleringer
- Kvanteproblemer

Løsning: Bruk datamaskinen til å beregne tilnærmede løsninger!

Hva er en numerisk metode?

En numerisk metode er en algoritme som:
- Deler opp et problem i små steg
- Beregner en tilnærmet løsning steg for steg
- Gir resultater med kontrollert presisjon

Fordeler:
- Kan løse problemer uten analytisk løsning
- Rask (med moderne datamaskiner)
- Visuell (kan lage grafer og animasjoner)

Ulemper:
- Tilnærmet (ikke eksakt)
- Krever programmering
- Kan akkumulere feil

Programmering i fysikk

Python er det mest populære språket for vitenskapelig programmering:
- Enkelt å lære
- Kraftige biblioteker (NumPy, Matplotlib)
- Gratis og åpen kildekode
- Brukes i forskning og industri

Hva skal vi lære:
1. Grunnleggende Python-syntaks
2. Eulers metode for differensiallikninger
3. Simulering av bevegelse
4. Visualisering av resultater

Differensiallikninger i fysikk

En differensialligning er en likning som inneholder deriverte.

Newtons 2. lov som differensialligning

Newtons 2. lov:
$$F = ma$$

Akselerasjon er den dobbeltderiverte av posisjon:
$$a = \frac{dv}{dt} = \frac{d^2s}{dt^2}$$

Derfor:
$$m\frac{d^2s}{dt^2} = F$$

Dette er en andreordens differensialligning.

Eksempel: Fritt fall uten luftmotstand

Kraft:
$$F = -mg$$

(Negativ fordi tyngdekraften peker nedover)

Differensialligning:
$$m\frac{d^2s}{dt^2} = -mg$$

$$\frac{d^2s}{dt^2} = -g$$

Analytisk løsning:
$$s(t) = s_0 + v_0 t - \frac{1}{2}gt^2$$

Denne kan vi løse med penn og papir! ✓

Eksempel: Fritt fall MED luftmotstand

Krefter:
- Tyngdekraft: $F_g = -mg$
- Luftmotstand: $\displaystyle F_d = \frac{1}{2}C_d \rho A v^2$ (oppover hvis nedover bevegelse)

Total kraft:
$$F = -mg + \frac{1}{2}C_d \rho A v^2$$

Differensialligning:
$$m\frac{dv}{dt} = -mg + \frac{1}{2}C_d \rho A v^2$$

Problem: Denne har INGEN analytisk løsning! ✗

Løsning: Bruk numeriske metoder! ✓

Omskriving til første orden

For å bruke Eulers metode omskriver vi til førsteordens system:

$$\frac{ds}{dt} = v$$

$$\frac{dv}{dt} = -g + \frac{C_d \rho A v^2}{2m}$$

Dette er to koblede førsteordens differensiallikninger.

Eulers metode

Eulers metode er den enkleste numeriske metoden for å løse differensiallikninger.

Grunnidé

Derivert = helling av tangent:
$$\frac{dy}{dx} = f(x, y)$$

Tilnærming:
Hvis vi kjenner $y$ ved $x$, kan vi estimere $y$ ved $x + \Delta x$:

$$y(x + \Delta x) \approx y(x) + \frac{dy}{dx} \cdot \Delta x$$

$$y(x + \Delta x) \approx y(x) + f(x, y) \cdot \Delta x$$

Eulers algoritme

Gitt:
- Differensialligning: $\displaystyle \frac{dy}{dt} = f(t, y)$
- Initialverdi: $y(t_0) = y_0$
- Tidssteg: $\Delta t$ (ofte skrevet $dt$)

Algoritme:

1. Start med $(t_0, y_0)$
2. Beregn $\displaystyle \frac{dy}{dt} = f(t_0, y_0)$
3. Estimér neste verdi: $y_1 = y_0 + f(t_0, y_0) \cdot \Delta t$
4. Oppdatér tid: $t_1 = t_0 + \Delta t$
5. Gjenta fra steg 2 med $(t_1, y_1)$

Generell formel:
$$y_{n+1} = y_n + f(t_n, y_n) \cdot \Delta t$$
$$t_{n+1} = t_n + \Delta t$$

Visualisering

Tenk på Eulers metode som å:
1. Stå på en kurve
2. Se tangentretningen (derivert)
3. Ta et lite steg i tangentretningen
4. Gjenta

Presisjon:
- Mindre $\Delta t$ → Bedre presisjon (men flere beregninger)
- Større $\Delta t$ → Raskere, men mindre nøyaktig

Eksempel: Eksponentiell vekst

Differensialligning:
$$\frac{dy}{dt} = ky$$

med $y(0) = 1$ og $k = 0.5$.

Analytisk løsning:
$$y(t) = e^{kt}$$

Numerisk løsning (Eulers metode):

La $\Delta t = 0.1$.

Steg 1:
- $t_0 = 0$, $y_0 = 1$
- $\displaystyle \frac{dy}{dt} = 0.5 \times 1 = 0.5$
- $y_1 = 1 + 0.5 \times 0.1 = 1.05$
- $t_1 = 0.1$

Steg 2:
- $t_1 = 0.1$, $y_1 = 1.05$
- $\displaystyle \frac{dy}{dt} = 0.5 \times 1.05 = 0.525$
- $y_2 = 1.05 + 0.525 \times 0.1 = 1.1025$
- $t_2 = 0.2$

Osv.

Sammenligning:

Python-grunnlag for fysikk

Variabler og beregninger

# Variabler
m = 0.5        # masse i kg
v = 10         # hastighet i m/s
g = 9.81       # tyngdeakselerasjon i m/s²

# Beregninger
Ek = 0.5 * m * v**2    # Kinetisk energi
Ep = m * g * h         # Potensiell energi
p = m * v              # Bevegelsesmengde

# Skriv ut
print("Kinetisk energi:", Ek, "J")
print(f"Bevegelsesmengde: {p} kg·m/s")  # f-string

Lister (arrays)

# Liste med tall
tid = [0, 1, 2, 3, 4, 5]
hastighet = [0, 2, 4, 6, 8, 10]

# Tilgang til elementer
print(tid[0])         # 0 (første element)
print(tid[-1])        # 5 (siste element)
print(hastighet[2])   # 4

# Legge til elementer
tid.append(6)
hastighet.append(12)

# Lengde av liste
n = len(tid)          # 7

Løkker

# For-løkke
for i in range(5):
    print(i)          # Skriver 0, 1, 2, 3, 4

# While-løkke
t = 0
while t < 5:
    print(t)
    t = t + 1         # Eller: t += 1

# Løkke over liste
hastigheter = [0, 2, 4, 6]
for v in hastigheter:
    Ek = 0.5 * m * v**2
    print(f"v = {v} m/s, Ek = {Ek} J")

Betingelser

v = 10

if v > 0:
    print("Beveger seg fremover")
elif v < 0:
    print("Beveger seg bakover")
else:
    print("I ro")

Funksjoner

# Definere funksjon
def kinetisk_energi(m, v):
    """Beregner kinetisk energi."""
    return 0.5 * m * v**2

# Bruke funksjon
Ek = kinetisk_energi(0.5, 10)
print(Ek)  # 25.0

NumPy for numeriske beregninger

import numpy as np

# Array (mer effektiv enn liste)
tid = np.array([0, 1, 2, 3, 4, 5])
hastighet = np.array([0, 2, 4, 6, 8, 10])

# Elementvise operasjoner
Ek = 0.5 * m * hastighet**2  # Beregner for alle verdier

# Nyttige funksjoner
print(np.max(hastighet))      # 10 (største verdi)
print(np.min(hastighet))      # 0 (minste verdi)
print(np.mean(hastighet))     # 5.0 (gjennomsnitt)
print(np.sum(hastighet))      # 30 (sum)

Matplotlib for plotting

import matplotlib.pyplot as plt

# Data
tid = [0, 1, 2, 3, 4, 5]
hastighet = [0, 2, 4, 6, 8, 10]

# Plott
plt.plot(tid, hastighet, 'b-')  # Blå linje
plt.xlabel('Tid (s)')
plt.ylabel('Hastighet (m/s)')
plt.title('Hastighet som funksjon av tid')
plt.grid(True)
plt.show()

Eulers metode i Python

Mal for Eulers metode

import numpy as np
import matplotlib.pyplot as plt

# Parametere
t0 = 0          # Starttid
y0 = 1          # Startverdi
t_slutt = 10    # Sluttid
dt = 0.1        # Tidssteg

# Derivert (differensialligning)
def f(t, y):
    """Returnerer dy/dt."""
    return -0.5 * y  # Eksempel: eksponentiell nedgang

# Initialiser lister
tid = [t0]
y = [y0]

# Eulers metode
t = t0
y_verdi = y0
while t < t_slutt:
    # Beregn derivert
    dydt = f(t, y_verdi)

    # Euler-steg
    y_verdi = y_verdi + dydt * dt
    t = t + dt

    # Lagre
    tid.append(t)
    y.append(y_verdi)

# Plott
plt.plot(tid, y)
plt.xlabel('Tid')
plt.ylabel('y')
plt.title('Løsning av differensialligning')
plt.grid(True)
plt.show()

Eksempel: Fritt fall

import numpy as np
import matplotlib.pyplot as plt

# Parametere
g = 9.81        # m/s²
s0 = 100        # Starthøyde (m)
v0 = 0          # Starthastighet (m/s)
dt = 0.01       # Tidssteg (s)

# Initialiser
tid = [0]
s = [s0]
v = [v0]

# Eulers metode
t = 0
s_verdi = s0
v_verdi = v0

while s_verdi > 0:  # Stopp når objektet treffer bakken
    # Akselerasjon
    a = -g

    # Euler-steg
    v_verdi = v_verdi + a * dt
    s_verdi = s_verdi + v_verdi * dt
    t = t + dt

    # Lagre
    tid.append(t)
    s.append(s_verdi)
    v.append(v_verdi)

# Resultat
print(f"Falltid: {t:.3f} s")
print(f"Slutthastighet: {v_verdi:.3f} m/s")

# Analytisk løsning
t_analytisk = np.sqrt(2 * s0 / g)
v_analytisk = -g * t_analytisk
print(f"Analytisk tid: {t_analytisk:.3f} s")
print(f"Analytisk hastighet: {v_analytisk:.3f} m/s")

# Plott
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))

# Høyde
ax1.plot(tid, s, 'b-', label='Numerisk')
ax1.set_xlabel('Tid (s)')
ax1.set_ylabel('Høyde (m)')
ax1.set_title('Høyde som funksjon av tid')
ax1.grid(True)
ax1.legend()

# Hastighet
ax2.plot(tid, v, 'r-', label='Numerisk')
ax2.set_xlabel('Tid (s)')
ax2.set_ylabel('Hastighet (m/s)')
ax2.set_title('Hastighet som funksjon av tid')
ax2.grid(True)
ax2.legend()

plt.tight_layout()
plt.show()