1.3 Numeriske metoder og programmering

Når penn og papir gir opp

I forrige kapittel så vi at en tennisball møter så mye luftmotstand at vi ikke kan ignorere den. Men da oppstår et nytt problem: hvis du faktisk tar med luftmotstanden, blir ligningen umulig å løse med penn og papir. Du sitter fast.

Dette er ikke et særtilfelle. Veldig mange spennende fysiske problemer mangler en analytisk løsning -- altså en pen formel du kan skrive ned. Fritt fall med luftmotstand, tre planeter som trekker i hverandre, værsimuleringer, kvanteproblemer: alle er for kompliserte for tradisjonell algebra. Løsningen er å la datamaskinen gjøre grovarbeidet.

En numerisk metode er en algoritme som deler opp et problem i mange små steg og beregner en tilnærmet løsning steg for steg, med kontrollert presisjon. Fordelen er enorm: du kan løse problemer som ikke har noen formel-løsning i det hele tatt, du får svaret raskt, og du kan tegne grafer og animasjoner. Ulempen er at svaret er tilnærmet -- aldri helt eksakt -- og at feil kan hope seg opp underveis.

Verktøyet vi bruker er Python, det desidert mest populære språket i vitenskapelig programmering. Det er lett å lære, gratis, og har kraftige biblioteker som NumPy for tallknusing og Matplotlib for grafer. Med litt Python kan du simulere fysikk som Galileo bare kunne drømme om.

Newtons lov, men som differensialligning

For å forstå numeriske metoder må vi se Newtons andre lov med nye øyne. En differensialligning er ganske enkelt en ligning som inneholder deriverte -- altså endringsrater.

Akselerasjon er jo endringsraten til farten, som igjen er endringsraten til posisjonen: $\displaystyle a = \frac{dv}{dt} = \frac{d^2 s}{dt^2}$. Setter vi dette inn i $F = ma$, får vi $\displaystyle m\frac{d^2 s}{dt^2} = F$. Dette er en andreordens differensialligning, og den er selve hjertet i mekanikken.

Noen ganger er den snill. For fritt fall uten luftmotstand er kraften $F = -mg$, så ligningen blir $\displaystyle \frac{d^2 s}{dt^2} = -g$, med den velkjente analytiske løsningen $\displaystyle s(t) = s_0 + v_0 t - \frac{1}{2}gt^2$. Penn og papir holder.

Men legg til luftmotstand, og snillheten forsvinner. Nå er $\displaystyle m\frac{dv}{dt} = -mg + \frac{1}{2}C_d \rho A v^2$, og denne har ingen analytisk løsning. For å angripe den numerisk skriver vi den om til et system av to koblede førsteordens ligninger: $\displaystyle \frac{ds}{dt} = v$ og $\displaystyle \frac{dv}{dt} = -g + \frac{C_d \rho A v^2}{2m}$. Posisjonen endres med farten, og farten endres med akselerasjonen. Det er denne formen datamaskinen elsker.

Eulers metode -- ett lite steg av gangen

Nå kommer den enkle, geniale ideen som lar datamaskinen løse ligningene: Eulers metode. Forestill deg at du står på en kurve og bare kan se i hvilken retning den peker akkurat der du står -- altså den deriverte, hellingen på tangenten. Da kan du gjette deg framover ved å ta et lite skritt i den retningen, sjekke hellingen på nytt, ta et nytt lite skritt, og slik fortsette.

Matematisk: hvis $\displaystyle \frac{dy}{dt} = f(t, y)$ og du kjenner $y$ akkurat nå, kan du anslå verdien et lite tidssteg $\Delta t$ senere som $y(t + \Delta t) \approx y(t) + f(t, y)\cdot \Delta t$. Skrevet som en oppskrift datamaskinen kan gjenta i det uendelige blir det $y_{n+1} = y_n + f(t_n, y_n)\cdot \Delta t$, og $t_{n+1} = t_n + \Delta t$.

La oss kjenne på det med eksponentiell vekst, $\displaystyle \frac{dy}{dt} = ky$ med $y(0) = 1$ og $k = 0{,}5$. Den analytiske løsningen er $y(t) = e^{kt}$, så vi kan sammenligne. Med $\Delta t = 0{,}1$: i første steg er $\displaystyle \frac{dy}{dt} = 0{,}5 \cdot 1 = 0{,}5$, så $y_1 = 1 + 0{,}5 \cdot 0{,}1 = 1{,}05$. I neste steg er $\displaystyle \frac{dy}{dt} = 0{,}5 \cdot 1{,}05 = 0{,}525$, så $y_2 = 1{,}05 + 0{,}525 \cdot 0{,}1 = 1{,}1025$. Ved $t = 1$ gir Euler 1,6289 mot den eksakte verdien 1,6487 -- et avvik på bare 1,2 %. Ikke verst for så enkel matematikk.

Et fint poeng: når akselerasjonen er konstant, blir Euler faktisk eksakt. En bil som akselererer med $a = 2$ m/s² fra hvile har farten $v = at$ -- en rett linje. Eulers metode følger tangenten, som her er en rett linje, så hvert steg $v_{n+1} = v_n + 2 \cdot \Delta t$ treffer perfekt: 0, 2, 4, 6, 8, 10 m/s etter fem sekunder, helt uten feil. Euler er eksakt for lineære problemer; det er bare når kurven faktisk krummer at feilen sniker seg inn.

Fra idé til Python-kode

Nå skal vi oversette tankegangen til kode. Du trenger overraskende lite Python for å komme i gang. Du lager variabler (m = 0.5, g = 9.81), du regner med dem (Ek = 0.5 * m * v**2, der ** betyr opphøyd i), og du lagrer mange verdier i lister eller NumPy-arrays. Du gjentar handlinger med løkker -- en while-løkke kjører så lenge en betingelse holder, perfekt for «kjør til ballen treffer bakken». Og du pakker gjenbrukbar kode i funksjoner med def.

Selve Eulers metode blir nesten en direkte avskrift av formelen. For fritt fall fra 100 m setter du startverdiene, og inne i løkken regner du akselerasjonen a = -g, oppdaterer farten v = v + a*dt, oppdaterer posisjonen s = s + v*dt, og øker tiden t = t + dt -- helt til s blir null. Underveis lagrer du alle verdiene i lister så du kan tegne dem etterpå med Matplotlib.

Den virkelige kraften viser seg når du legger til luftmotstand. Ta en ball på 0,5 kg med radius 0,05 m sluppet fra 100 m. Tverrsnittet er $A = \pi r^2 \approx 0{,}00785$ m². Akselerasjonen er ikke lenger konstant, men $\displaystyle a = -g + \frac{C_d \rho A v^2}{2m}$ -- den avhenger av farten, og det er nettopp derfor ingen pen formel finnes. Datamaskinen bryr seg ikke; den regner bare steg for steg.

Resultatet er lærerikt. Med luftmotstand bruker ballen rundt 4,87 s på fallet og treffer med cirka 42 m/s, mens uten luftmotstand er det 4,52 s og 44 m/s. Luftmotstanden gjør altså fallet både lengre og litt langsommere. Den teoretiske terminalhastigheten -- der luftmotstand og tyngde balanserer -- er $\displaystyle v_{term} = \sqrt{\frac{2mg}{C_d \rho A}} \approx 71$ m/s, men ballen rekker ikke å komme dit på bare 100 m. På en hastighet-tid-graf ser du dette tydelig: kurven med luftmotstand begynner å flate ut og nærme seg terminalhastigheten, mens den uten luftmotstand fortsetter rett opp som en strek.

Tidssteget -- den lille knappen som avgjør alt

Det finnes én skrue du absolutt må forstå når du jobber numerisk: tidssteget $\Delta t$ (i kode ofte kalt dt). Det avgjør hele kvaliteten på simuleringen.

Velger du et for lite tidssteg, blir resultatet svært nøyaktig -- men beregningen blir treg, fordi du må gjøre enormt mange steg. Velger du et for stort tidssteg, går det fort, men nøyaktigheten ryker, og i verste fall blir løsningen direkte ustabil og eksploderer. Den gode framgangsmåten er å starte lite (for eksempel $dt = 0{,}001$), prøve flere verdier, og sjekke konvergens: når svaret nesten ikke endrer seg om du halverer $dt$, er du på trygg grunn. For et fall fra 100 m som tar rundt 4,5 s bør du bruke minst tusen steg, altså $dt$ omkring 0,0045 s eller mindre.

Du må også være klar over at Eulers metode bare er en tilnærming, og at feil akkumuleres. Det kommer av tre ting: avrundingsfeil fordi datamaskinen har begrenset presisjon, diskretiseringsfeil fordi vi erstatter en jevn, kontinuerlig prosess med hakkete små steg, og rene stabilitetsproblemer ved for store tidssteg.

Hvordan demmer du opp? Bruk mindre tidssteg, valider mot analytiske løsninger der de finnes, og test konvergens. For virkelig krevende oppgaver -- som å simulere planetbaner over millioner av år -- er Euler altfor grov; der ville de små feilene hope seg opp til at planeten «flyr ut av solsystemet». Da trenger man mer avanserte metoder som Runge-Kutta, som du møter senere. Men prinsippet er det samme: ett lite, kontrollert steg av gangen.

Oppsummering

Når penn og papir gir opp, tar datamaskinen over. En numerisk metode løser et problem ved å dele det i mange små steg og regne seg fram til en tilnærmet løsning -- uunnværlig for problemer uten analytisk løsning, som fritt fall med luftmotstand.

Nøkkelen er å se Newtons andre lov som en differensialligning, $\displaystyle m\frac{d^2 s}{dt^2} = F$, og skrive den om til et system av førsteordens ligninger, $\displaystyle \frac{ds}{dt} = v$ og $\displaystyle \frac{dv}{dt} = a$. Eulers metode løser dem ved å ta ett lite steg langs tangenten av gangen: $y_{n+1} = y_n + f(t_n, y_n)\cdot \Delta t$. Den er eksakt for lineære problemer som konstant akselerasjon, og en god tilnærming ellers.

I Python blir dette nesten en direkte avskrift av formelen, med variabler, lister, løkker og funksjoner -- og med NumPy og Matplotlib til tallknusing og grafer. Simuleringen av en fallende ball avslørte både lengre falltid og lavere slutthastighet med luftmotstand, og pekte mot en terminalhastighet på rundt 71 m/s. Til slutt lærte vi den viktigste skruen: tidssteget $\Delta t$. For lite gir treghet, for stort gir unøyaktighet og ustabilitet. Sjekk alltid konvergens, vær bevisst på at feil akkumuleres, og husk -- god fysikksimulering handler om mange små, kontrollerte steg.