In [1]:
from IPython.display import HTML
HTML('''<script>
code_show=true; 
function code_toggle() {
 if (code_show){
 $('div.input').hide();
 } else {
 $('div.input').show();
 }
 code_show = !code_show
} 
$( document ).ready(code_toggle);
</script>
<form action="javascript:code_toggle()"><input type="submit" value="Visa / dölj källkoden."></form>''')
Out[1]:
In [3]:
# Först importeras och definieras några saker som används i hela dokumentet
import matplotlib.pyplot as plt
import numpy as np
import math as math
from scipy.integrate import odeint

plt.rcParams['figure.figsize'] = (16., 8.0)
plt.rcParams.update({'font.size': 14})
plt.rcParams['lines.linewidth'] = 5
plt.rcParams['font.family'] = 'serif'

g = -9.82 # Tyngdaccelerationen, sätts till negativt värde då läget ökar med höjd över marken.
Nikodemus Karlsson, gymnasielärare i matematik och fysik
2019-04-22

Numeriska lösningar i Python till några ordinära differentialekvationer

1 Bakgrund

Det förra Jupyter-dokument jag skrev utgick från riktiga mätdata från min telefons accelerometer och trycksensor. Utifrån dessa kunde jag processa beräkningar för att bl a rita upp grafer för läge, hastighet och acceleration samtidigt som lufttrycket registrerades under en hissfärd uppåt och nedåt. Det var kul, syftet var att få en inblick i hur man kan hantera och visualisera numeriska data med hjälp av Python.

Nu fortsätter jag med ett annat tema: differentialekvationer. Jag ville ta reda på hur man numeriskt kan lösa några av de differentialekvationer som elever kommer i kontakt med i gymnasiets fysikkurser med hjälp av Python. Jag kommer inte att skapa en egen stegmetod för att differentialekvationerna numeriskt. Istället är utgångspunkten en funktion som finns tillgänglig i biblioteket scipy: nämligen odeint(). Här kommer jag att utgå från befintliga modeller och inte från faktiska mätdata. Det gör lösningarna, i absoluta tal, till differentialekvationerna inte nödvändigtvis behöver relatera till några verkliga förhållanden. Däremot så framgår likheter och skillnader mellan de använda modellerna samtidigt som jag påvisar och kommenterar den Python-kod som löser differentialekvationerna.


2 Fall utan luftmotstånd

2.1 Fritt fall

Den kanske första differentialekvationen gymnasieelever kommer i kontakt med presenteras i det sammanhanget inte som en differentialekvaton. Det har att göra med acceleration, och en tidig tillämpning är ett föremåls fria fall (då ett föremål påverkas av enbart en kraft; dess tyngd). Under den första fysikkursen tas iallafall upp att föremål som faller med tyngdaccelerationen ökar sin fart med 9.82 m/s varje sekund.

2.1.1 Differentialekvation

Accelerationen definieras som hastighetsförändring per tidsenhet. Givet att det är enkom föremålets tyngd som ger upphov till accelerationen kan vi teckna

$$\frac{dv}{dt}=\frac{F}{m}=\frac{mg}{m}=g$$

Även om som sagt kanske det första eleverna kommer i kontakt med inte är denna differentialekvation, så ingår absolut dess lösning, $v=v_0+gt$, redan i första fysikkursen. Själva differentialekvationen säger också något viktigt: den definierar (tyngd)accelerationen som förändring av hastighet över tid och det appliceras på Newtons andra lag.

I detta dokument kommer jag att lösa de aktuella differentialekvationerna numeriskt med Pythons funktion odeint(). Denna genererar en lista med $y$-värden (eller $v$-värden i detta fall) som sedan kan ritas upp som en graf.

2.1.1 Integrerande kod

# Funktionen som definierar själva differentialekvationen
def freeFall(Y, t):
    # Variabeln Y som kommer med odeint():s anrop, används inte i detta exempel, dock senare.
    dvdt = g # Själva differentialekvationen
    # konstanten g måste naturligvis vara definierad sedan tidigare (g = -9.82)

    # Här sker det magiska: dvdt integreras och värdena återbördas till
    # som en lista [v0, v1, v2,...,v100]
    return(dvdt)

# Anropet till freeFall() ovan sköts av odeint(). odeint() måste dock ha 
# ett startvärde och en lista med värden mellan vilken förändringen sker.
# Den aktuella t-listan i koden nedan innehåller 100 jämnt fördelade värden
# mellan 0 och 4, det ger att dt ≈ 0.04. Ju tätare värden (dt närmare noll),
# desto bättre lösning. I just detta fall behövs egentligen inte så täta värden
# eftersom lösningen blir en rät linje.
v = odeint(freeFall, v0, t) # v innehåller nu värdena [v0, v1, v2,...,v99]

2.1.2 Lösning

För denna vår första differentialekvation blir några lösningar (beroende på begynnelsevärde)

In [4]:
def freeFall(Y, t):
    dvdt = g
    return(dvdt)

t = np.linspace(0, 4, 100)
for v0 in [-10, 0, 10]:
    v = odeint(freeFall, v0, t)
    plt.plot(t, v, label = "v0 = {0:.1f} m/s".format(v0)) # En graf plottas för varje v0

plt.title("Hastighet hos fritt fallande föremål, olika begynnelsehastigheter", y=1.0)
plt.xlabel("Tid [sekunder]")
plt.ylabel("Hastighet [m/s]")
plt.grid()
plt.legend();

2.1.3 Kommentar

I grafen ovan, liksom i fortsättningen, har jag satt positiv riktning uppåt. Det innebär jag t ex i fallet med begynnelsehastigheten 10 m/s tänker mig att ett föremål kastas uppåt med farten 10 m/s för att sedan falla ned igen (i detta fall även till en lägre nivå än ursprungsnivån).


3 Fall med luftmotstånd

3.1 Ingen begynnelsehastighet

Det föremål som faller i lösningen ovan faller fritt, det påverkas det enbart av gravitationskraften. I våra vanliga miljöer är det i praktiken så att även luftmotståndet bidrar till den resulterande kraften (och därmed accelerationen) på det fallande föremålet. Luftmotståndet är berorende av hastigheten, ju högre hastighet desto större luftmotstånd. Gymnasiefysiken lär oss från början att massan på föremålet inte inverkar på accelerationen. Det är sant enbart då gravitationen bidrar till accelerationen. Dock blir det betydligt svårare att räkna på hur hastigheten påverkas då en kraft som i sin tur beror på hastigheten påverkar just hastigheten. En sådan modell ingår inte i gymnasiefysiken att hantera algebraiskt, även om själva differentialekvationen inte blir speciellt mycket svårare att teckna (men väl svårare att lösa analytiskt).

3.1.1 Differentialekvation

Vår differentialekvation som ska beskriva hur den resulterande kraften på ett föremål förändras när när det påverkas av såväl gravitationskraft och luftmotstånd, där luftmotståndet i sin tur beror på hastigheten, blir

$$m\frac{dv}{dt}=mg-kv^2, \text{där } {v(0)=0}$$

där $k$ är en konstant som säger något om luftmotståndet (den beror på föremålets form och atmosfäriska förhållanden).

Det ger accelerationen

$$\frac{dv}{dt}=g-\frac{k}{m}v^2$$

Här ser ett tränat öga att accelerationen minskar med ökad fart, vi ser även att massan fatiskt ingår som en parameter för accelerationen. I det här fallet(!) så används att luftmotståndet är proportionellt mot farten i kvadrat, det är en modell som stämmer tillräckligt bra överens med verkliga förhållanden (även om man kan gå längre och studera hur farten påverkar luftströmmar runt om föremålet som i sin tur skapar andra krafter...) för att skapa oss en bild av hur farten och luftmotstånd hänger ihop.

3.1.2 Integrerande kod

# Funktionen som definierar själva differentialekvationen
def fallWithDrag(Y, t):
    # Y är den variabel som integrationen skett på, dvs hastigheten i det här fallet.
    # Derivatan är ju dvdt i detta fall, efter integration blir det Y...
    v = Y
    # ...som ansätts till v för tydligare hantering.

    # Y (eller v) används (om aktuellt) som funktionsvariablen i differentialekvationen.
    # I vårt fall så får vi
    dvdt = g - k/m*v*np.abs(v)
    return(dvdt) # Let the magic happen!

# Slutligen så anropas fallWithDrag() från odeint() enligt föregående exempel
v = odeint(fallWithDrag, v0, t)

3.1.3 Lösning

In [5]:
def fallWithDrag(Y, t):
    v = Y 
    dvdt = g - k/m*v*np.abs(v) # Derivatan av v (eller Y)
    return(dvdt)

k = 0.05  # Luftmotståndskoefficient
m = 1 # Massa
v0 = 0 # Startvillkor / begynnelsehastighet
t = np.linspace(0, 4, 100) # 100 jämnt fördelade värden mellan 0 och 4.
v = odeint(fallWithDrag, v0, t) # Skapar en lista med ett v-värde för respektive t-värde

for m in [0.2, 0.4, 0.6, 0.8, 1.0]:
    v = odeint(fallWithDrag, v0, t)
    plt.plot(t, v, label = "m = {0:.1f} kg".format(m))

plt.title("Hastighet hos fallande föremål med luftmotstånd,\n ingen begynnelsehastighet", y=1.03)
plt.xlabel("Tid [sekunder]")
plt.ylabel("Hastighet [m/s]")
plt.grid()
plt.legend();

3.1.4 Kommentar

Se där! Vi ser att föremål med större massa får en högre fart än de med mindre. Det stämmer bra överens med det vi är vana vid att se i naturen, t ex att kottar från träd kan komma ned riktigt fort medan löven singlar ned långsamt.

Nu går vi över och studerar hur ett föremåls läge, hastighet och acceleration kommer att utvecklas över tiden det får en begynnelsehastighet, t ex genom ett kast uppåt.

3.2 Med begynnelsehastighet

Jag tänker här åskådliggöra såväl höjd, hastighet som acceleration då ett föremål kastas uppåt varpå enbart tyngdkraft och luftmotståndet verkar på föremålet. För att kunna få fram alla dessa tre rörelsegrafer utgår jag från derivatan av accelerationen och låter odeint() integrera tre gånger enligt nedan.

3.2.1 Differentialekvationer

In [6]:
# Koden nedan fogar enbart in en bild. Finns enklare sätt, men det här sättet gav en skarp bild
# som dessutom hänger med i HTML-export.
# Denna lösning på att visa en bild med matplotlib fann jag på https://stackoverflow.com/a/42314798
def display_image_in_actual_size(im_path):

    dpi = 80
    im_data = plt.imread(im_path)
    height, width, depth = im_data.shape

    # What size does the figure need to be in inches to fit the image?
    figsize = width / float(dpi), height / float(dpi)

    # Create a figure of the right size with one axes that takes up the full figure
    fig = plt.figure(figsize=figsize)
    ax = fig.add_axes([0, 0, 1, 1])

    # Hide spines, ticks, etc.
    ax.axis('off')

    # Display the image.
    ax.imshow(im_data, cmap='gray')

    plt.show()

display_image_in_actual_size("./ode_avy_1d.png")

# Annan variant, ger inte lika skarp bild
# import matplotlib.image as mpimg
# plt.rcParams['figure.figsize'] = (16., 15.0)
# image = mpimg.imread("Diffar.png")
# plt.axis("off")
# plt.imshow(image)
# plt.show()

3.2.2 Integrerande kod

# Funktionen som definierar själva differentialekvationen
def fallWithSquaredDrag(Y, t, k):
    # Y kommer här att innehålla tre komponenter pga de tre startvillkoren...
    y, v, a = Y
    # ...som utgör de sökta storheterna och tilldelas y, v och a.

    # v-faktorn nedan ska alltid vara positiv i modellen nedan, varpå hela faktorn blir negativ,
    # eftersom luftmotståndet verkar *mot* rörelseriktningen.
    dadt = [v, a, -2*k/m*(np.abs(v)*a)]

    # odeint() tar hand om derivatorerna i samband med nedanstående return()...
    return(dadt)
    # ...varpå dessa integreras och återbördas (bakom kulisserna) till det som har anropat odeint()
    # Skillnaden mot tidigare exempel är att det är hela *systemet* av differentialekvationer
    # som integreras!

# Funktionen ovan anropas av odeint() (nedan) med tillhörande parametrar
# init_squared är en vektor som innehåller de tre begynnelsevillkoren [y0, v0, a0]
# t är listan med tidpunkter, ungefär [0, 0.01, 0.02,...,10]
# medan k_squared är det specifika k-värdet.
# y_squared kommer efter lyckat anrop att vara en lista med listor:
# y_squared[0] utgörs av lägeskoordinater vid respektive t-värde
# y_squared[1] utgörs av hastighetskoordinater vid respektive t-värde
# y_squared[2] utgörs av accelerationskoordinater vid respektive t-värde
y_squared = odeint(fallWithSquaredDrag, init_squared, t, args=(k_squared,))

3.2.3 Lösning

In [7]:
# Denna kod skapar funktioner för acceleration, hastighet och läge som kan visas i grafer
def fallWithLinearDrag(Y, t, k):
    y, v, a = Y
    dadt = [v, a, -k/m*a]
    return(dadt)

def fallWithSquaredDrag(Y, t, k):
    y, v, a = Y
    dadt = [v, a, -2*k/m*(np.abs(v)*a)]
    return(dadt)

# Gemensamt för modellerna
y0 = 0
v0 = 20
m = 2

# Parametrar som förändras med modellerna
k_noDrag = 0
k_linear = 1
k_squared = 0.05
a0_noDrag = g
a0_linear = g - k_linear/m*v0
a0_squared = g - k_squared/m*v0*np.abs(v0)
init_noDrag = [y0, v0, a0_noDrag]
init_linear = [y0, v0, a0_linear]
init_squared = [y0, v0, a0_squared]

# Skapar listan med tider...
t = np.linspace(0, 10, 100)

# ...och värden för fallet utan luftmotstånd...
y_noDrag = odeint(fallWithLinearDrag, init_noDrag, t, args=(k_noDrag,))

# ...värden för fallet med linjärt beroende luftmotstånd...
y_linear = odeint(fallWithLinearDrag, init_linear, t, args=(k_linear,))

# ...och värden för fallet med kvadratiskt luftmotstånd.
y_squared = odeint(fallWithSquaredDrag, init_squared, t, args=(k_squared,))
In [8]:
# Nedanstående kod ritar enbart upp graferna efter de listor som tagits fram ovan.
# Det görs med matplotlib:s subplot(), som kan lägga upp flera diagram i ett rutnät. Tekniskt sett
# är alltså nedanstående nio diagram endast ett diagram med respektive graf utlagd!

fig, ax = plt.subplots(3, 3)
plt.suptitle("Föremåls rörelse i g-fält, med och utan luftmotstånd\nGraferna visar en och samma situation utifrån olika modeller", y = 1.0)

# Kolumn 1
ax1 = plt.subplot(331)
plt.title("↓ Utan luftmotstånd ↓", y = 1.04)
plt.plot(t, y_noDrag[:,0])
ax1.set_ylim(-120, 25)
plt.setp(ax1.get_xticklabels(), visible=False)
plt.setp(ax1.get_yticklabels(), visible=True)
plt.grid()
plt.ylabel("Höjd [m]")

ax2 = plt.subplot(334)
ax2.set_ylim(-20,22)
plt.plot(t, y_noDrag[:,1])
plt.ylabel("Hastighet [m/s]")
plt.setp(ax2.get_xticklabels(), visible=False)
plt.setp(ax2.get_yticklabels(), visible=True)
plt.grid()

ax3 = plt.subplot(337)
plt.xlabel("Tid [s]")
plt.plot(t, y_noDrag[:,2])
plt.ylabel("Acceleration [m/s$^2$]")
plt.setp(ax3.get_xticklabels(), visible=True)
plt.setp(ax3.get_yticklabels(), visible=True)
plt.grid()

# Kolumn 2
ax4 = plt.subplot(332, sharex=ax1, sharey=ax1)
plt.title("↓ Luftmotstånd, linjär modell ↓", y = 1.04)
plt.ylim(bottom = -120, top = 25)
plt.plot(t, y_linear[:,0])
plt.setp(ax4.get_xticklabels(), visible=False)
plt.setp(ax4.get_yticklabels(), visible=False)
plt.grid()


ax5 = plt.subplot(335, sharex=ax2, sharey=ax2)
plt.plot(t, y_linear[:,1])
plt.setp(ax5.get_xticklabels(), visible=False)
plt.setp(ax5.get_yticklabels(), visible=False)
plt.grid()

ax6 = plt.subplot(338, sharex=ax3, sharey=ax3)
plt.xlabel("Tid [s]")
plt.plot(t, y_linear[:,2]) 
plt.setp(ax6.get_xticklabels(), visible=True)
plt.setp(ax6.get_yticklabels(), visible=False)
plt.grid()

# Kolumn 3
ax7 = plt.subplot(333, sharex=ax1, sharey=ax1)
plt.ylim(bottom = -120, top = 25)
plt.title("↓ Luftmotstånd, kvadratisk modell ↓", y = 1.04)
plt.plot(t, y_squared[:,0])
plt.setp(ax7.get_xticklabels(), visible=False)
plt.setp(ax7.get_yticklabels(), visible=False)
plt.grid()

ax8 = plt.subplot(336, sharex=ax2, sharey=ax2)
plt.grid()
plt.setp(ax8.get_xticklabels(), visible=False)
plt.setp(ax8.get_yticklabels(), visible=False)
plt.plot(t, y_squared[:,1])

ax9 = plt.subplot(339, sharex=ax3, sharey=ax3)
plt.xlabel("Tid [s]")
plt.plot(t, y_squared[:,2])
plt.setp(ax9.get_xticklabels(), visible=True)
plt.setp(ax9.get_yticklabels(), visible=False)
plt.grid()

#plt.show()
plt.subplots_adjust(top=0.85) # Hyfsa till läget så att rubriken syns vid export
plt.savefig("air_drag.png", transparent=True)

3.2.4 Kommentar

I första kolumnen av graferna nedan är föremålet i fritt fall (enbart påverkat av tyngdkraften) medan de båda andra kolumnerna ger en jämförelse mellan luftmotstånd som påverkar föremålet direkt proportionellt mot farten (det gör att man kommer betydligt närmare verkliga förhållanden och beräkningen bör kunna göras inom ramen för gymnasiets matematik) och kvadratiskt proportionellt (avsevärt svårare att beräkna en lösning analytiskt, men gör också att lösningen kommer ännu närmare verkliga förhållanden).

Utifrån graferna kan läsaren ställa sig åtminstone två frågor: Vad finns det för skillnaden mellan den linjära modellen och den kvadratiska och hur kan det komma sig att accelerationen startar på ett värde som är lägre än $g$ i modellerna med luftmotstånd?

3.3 Kastbanor i luftmotstånd

Skillnaden mellan de tidigare graferna och de kommande kastbanorna är att de tidigare visade ett läge som funktion av tiden. Här ska vi istället ta fram en funktion som visar läget i $x$-led som funktion av läget i $y$-led. Vi kommer alltså att övergå till ett snett kast snarare än kastet rakt uppåt i avsnitt 3.2.

3.3.1 Differentialekvationer

In [9]:
display_image_in_actual_size("./ode_vy_2d.png")