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

Innehåll

Bakgrund

Fall utan luftmotstånd

Fall med luftmotstånd, utan begynnelsehastighet

Fall med luftmotstånd, med begynnelsehastighet

Kastbanor med luftmotstånd

Harmonisk svängning i fjäder

Matematisk pendel

Referenser

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")

3.3.2 Integrerande kod

# Funktionen som definierar själva differentialekvationen
def trajectory(Y, t):
    # De fyra startvillkoren ger fyra Y-värden. Dock är det inte frågan om en fjärdederivator;
    # det är frågan om två andraderivatorer och två förstaderivatorer.
    sx, sy, vx, vy = Y
    # Alternativ tilldelning, eftersom sx och sy inte används:
    # vx,vy = Y[2:]
    # men det blir inte lika tydligt!

    # Differentialekvationerna definieras enligt rutan ovan. Observera att sx och sy
    # som sagt inte ingår.
    dvxdt = ax - k/m*vx*np.abs(vx)
    dvydt = ay - k/m*vy*np.abs(vy)

    # Slutligen integreras allt 
    return [vx, vy, dvxdt, dvydt]

# Anropet till odeint() sker på samma sätt som tidigare med alla startvillkor y0
y = odeint(trajectory, y0, t)
# där y[0] kommer att innehålla en lista med x-koordinater
# och y[1] en lista med korresponderande y-koordinater.
# y[3] och y[4] kommer att innehålla korresponderande hastigheter i x- respektive y-led,
# men dessa är vi intresserade av att hantera här.

3.3.3 Lösning

In [10]:
def trajectory(Y, t):
    sx, sy, vx, vy = Y
    dvxdt = ax - k/m*vx*np.abs(vx)
    dvydt = ay - k/m*vy*np.abs(vy)
    return [vx, vy, dvxdt, dvydt]

# Nedanstående funktion omvandlar polära koordinater till rektangulära.
# Behövs för att enkelt kunna ange fart och riktning, medan den funktion som odeint()
# anropar är designad för rektangulära koordinater.
def pol2rect(r, theta, deg = True):
    if(deg):
        theta = theta * 2 * math.pi / 360
    x = r * math.cos(theta)
    y = r * math.sin(theta)
    return[x, y]

t = np.linspace(0, 3, 45)
sx0 = 0
sy0 = 0
ax = 0
ay = -9.82
m = 1

speed = 10
elevation = 60
vx0, vy0 = pol2rect(r = speed, theta = elevation)
y0 = [sx0, sy0, vx0, vy0]

for k in [0., 0.05, 0.1, 0.2, 0.4, 0.8]:
    y = odeint(trajectory, y0, t)
    graphs = plt.plot(y[:,0], y[:,1],'o-', ms = 13, mec = 'w', mew = 3, label = "k = {0:.2f}".format(k));
plt.title("Kastbanor, luftmotståndskoefficienten varierar", y=1.0)
plt.xlabel("Läge i x-led [m]")
plt.ylabel("Läge i y-led [m]")
plt.legend();
plt.xlim(left = 0, right = 10)
plt.ylim(bottom = 0, top = 4)
plt.grid();
plt.text(8.05, 2.2, "Utkastfart: {0:.0f} m/s\nUtkastvinkel: {1:.0f}°".format(speed, elevation),
         fontsize=14, horizontalalignment='left');
In [11]:
speed = 10
elevation = 30
vx0, vy0 = pol2rect(r = speed, theta = elevation)
y0 = [sx0, sy0, vx0, vy0]

for k in [0., 0.05, 0.1, 0.2, 0.4, 0.8]:
    y = odeint(trajectory, y0, t)
    graphs = plt.plot(y[:,0], y[:,1],'o-', ms = 13, mec = 'w', mew = 3, label = "k = {0:.2f}".format(k));
plt.title("Kastbanor, luftmotståndskoefficienten varierar", y=1.0)
plt.xlabel("Läge i x-led [m]")
plt.ylabel("Läge i y-led [m]")
plt.legend();
plt.xlim(left = 0, right = 10)
plt.ylim(bottom = 0, top = 4)
plt.grid();
plt.text(8.05, 2.2, "Utkastfart: {0:.0f} m/s\nUtkastvinkel: {1:.0f}°".format(speed, elevation),
         fontsize=14, horizontalalignment='left');
plt.text(3, 3.05, "Diagrammet har avsiktligt samma skala som det\nsom visar banor med utkastvinkeln 60°.\nI övrigt är det identiska data, vilket gör dem jämförbara.");
In [12]:
y0 = [sx0, sy0, vx0, vy0]
k=0.4
elevation = 30

for speed in [10, 15, 20, 30, 40]:
    vx0, vy0 = pol2rect(r = speed, theta = elevation)
    y0 = [sx0, sy0, vx0, vy0]
    y = odeint(trajectory, y0, t)
    plt.plot(y[:,0], y[:,1],'o-', ms = 13, mec = 'w', mew = 3, label = "$v_0$ = {0:.0f} m/s".format(speed));
plt.title("Kastbanor, utkastfarten varierar", y=1.0)
plt.xlabel("Läge i x-led [m]")
plt.ylabel("Läge i y-led [m]")
plt.legend();
plt.xlim(left = 0, right = 10)
plt.ylim(bottom = 0, top = 4)
plt.grid();
plt.text(7.4, 2.1, "Luftmotståndskoefficient: {0:.1f}\nUtkastvinkel: {1:.0f}°".format(k, elevation),
         fontsize=14, horizontalalignment='left');

3.3.4 Kommentar

Tidpunkten för respektive läge kan inte direkt avläsas i graferna, men avståndet mellan två "prickmarkeringar" är ett och samma tidsintervall (1/15 sekund). Det syns t ex att prickarna ligger tätare i vändpunkten just därför att föremålet inte hinner lika långt under tidsintervallet när det färdas långsammare. Luftmotståndet påverkar såklart mot aktuell rörelseriktning, medan tyngdaccelerationen enbart påverkar nedåt.

I simuleringarna syns det att både luftmotståndskoefficienten $k$ och hastigheten påverkar banans utseende. I samtiga fall används modellen där luftmotståndet är proportionellt mot hastigheten i kvadrat. Notera att modellen utan luftmotstånd ($k$ = 0) ger samma kastvidd vid vinklar vars summa är 90° och samma utkastfart; så är inte fallet med luftmotstånd (här bör läsaren ställa sig frågan varför föremålet kommer längre då det kastas ut med vinkeln 30° jämfört med 60° för en och samma luftmotståndskoefficient).

4 Harmoniska (och nästan harmoniska) svängningar

4.1 Harmonisk svängning i fjäder

Att ett föremål beskriver en harmonisk svängning innebär att det finns en återställande kraft som är proportionell mot avståndet till jämviktsläget. I gymnasiefysiken introduceras det ofta med hjälp av en vikt som hänger i en fjäder: enligt Hookes lag kommer kraften från fjädern att verka proportionellt mot förlängningen.

4.1.1 Differentialekvationer

Differentialekvationen som beskriver läget $y$ vid olika tidpunkter $t$ för ett föremål i harmonisk svängning (i en dimension) blir

$$\frac{d^2y}{dt^2}=-\frac{k}{m}y$$

med föremålets massa $m$, läget i förhållande till jämviktsläget $y$ och fjäderkonstanten $k$.

Denna differentialekvation är förhållandevis enkel att lösa (såväl analytiskt och med Pythons:s odeint(). Jag tänkte ta den ytterligare ett steg längre: genom att derivera ytterligare en gång kan vi i lösningen även få fram accelerationen, utöver hastigheten och läget. Dessa tre storheter kan plottas i ett och samma diagram. Den aktuella differentialekvationen vi ska låta odeint() lösa blir då

$$\frac{d^3y}{dt^3}=-\frac{k}{m}\frac{dy}{dt}$$

I just detta fall tänker jag också att det intressanta är att se hur dessa tre storheter kommer att förhålla sig till varandra rent fasmässigt, varför jag ansätter $\frac{k}{m}=1$. Det ger oss alltså att lösa ekvationssystemet

$$\frac{d^3y}{dt^3}=\frac{d^2v}{dt^2}=\frac{da}{dt}=-\frac{dy}{dt}$$

Tre integreringar ger i tur och ordning acceleration, hastighet och läge. Startvillkoren sätts till $y(0)=-1$, $v(0)=0$ och $a(0)=1$.

4.1.2 Integrerande kod

# Funktionen som definierar själva differentialekvationen
def spring(Y, t):
    y, dydt, d2ydt2 = Y # Variabler för läge, hastighet och acceleration ansätts...
    return [dydt, d2ydt2, -dydt] # ...och integreras

# Funktionsanropet sker på vanligt sätt med odeint()
y = odeint(spring, y0, t)
# där y[0] blir en lista på lägen, y[1] på hastigheter och y[2] på accelerationer över tiden.

4.1.3 Lösning

In [13]:
def spring(Y, t):
    y, dydt, d2ydt2 = Y
    return [dydt, d2ydt2, -dydt]

t = np.linspace(0, 10, 100)

y_init = -1
v_init = 0
a_init = 1

y0 = [y_init, v_init, a_init]
y = odeint(spring, y0, t)
plt.grid()
# Y = np.c_[y, -y[:,0]]
graphs = plt.plot(t, y);
plt.legend(graphs, ('Läge', 'Hastighet', 'Acceleration'),loc='upper center',
           bbox_to_anchor=(1.01, 1.05),
           fancybox=True, shadow=True, ncol=1);
plt.ylabel("Relativa värden för respektive storhet");
plt.xlabel("Tid [s]")
plt.title("Relation mellan läge, hastighet och acceleration i harmonisk svängning", y = 1.05);

4.1.4 Kommentar

Snyggt! Här syns t ex att där farten är maximal är accelerationen noll och att det sker i jämviktsläget!

4.2 Matematisk pendel

En matematisk pendel är en tyngd som hänger i ett snöre utan massa. Den oscillerar i ett plan i gravitationsfältet. Den beskriver en approximativt harmonisk svängningsrörelse; den återställande kraften är proportionell mot sinusvärdet för det aktuella utslaget. För små vinklar gäller att $\sin(\theta)\approx \theta$. Men ju större vinkeln är och ju längre tiden går, desto större blir felet med den harmoniska modellen. Den harmoniska modellen är trevlig just därför att den är enkel. Modellen som tar hänsyn till vinkeln är betydligt svårare att lösa analytiskt. Om man dessutom lägger på ett luftmotstånd så blir det riktigt jobbigt. Men som du säkert redan gissat; vi låter här odeint() göra jobbet!

4.2.1 Differentialekvationer

Den förenklade, linjära, modellen tecknas enligt

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

medan den ickelinjära modell som tar hänsyn till vinkel och luftmotstånd tecknas enligt

$$\frac{d^2\theta}{dt^2}=-k\omega^2-\frac{g}{r}\sin\theta$$

där $\theta$ är det aktuella vinkelutslaget, $k$ är luftmotståndskoefficienten, $\omega$ är den aktuella vinkelhastigheten (tidsderivatan av $\theta$), $g$ är tyngdacceleratonen och slutligen $r$ som är pendelns konstanta avstånd till sin upphängningspunkt. Det vi vill se är hur vinkeln $\theta$ förändras med tiden och jämföra den förenklade modellen med den ickelinjära modellen för olika startvärden på vinkelutslaget. Låt oss nu titta på det centrala i koden.

4.2.2 Integrerande kod

# Förenklad modell som utgår från att pendeln beskriver
# en harmonisk svängning.
def simpleModel(Y, t):
    theta, omega = Y
    # Jag väljer r till 0.2487 m, ty det ger periodtiden 1 s i den förenklade modellen.
    retVal = [omega, - g/r*theta]
    return(retVal)

# Någorlunda korrekt modell med hänsyn taget
# till luftmotstånd och vinkel som gör att
# svängningen inte blir harmonisk.
def properModel(Y, t):
    theta, omega = Y
    retVal = [omega, -k*omega*np.abs(omega) - g/r*np.sin(theta)]
    return(retVal)

# Och slutligen anropen till respektive funktion, där y0 är startvillkoret,
# som innehåller startvärden för θ (varierar) och ω (sätts till noll i samtliga fall).
y1 = odeint(simpleModel, y0, t) # Förenklad modell
y2 = odeint(properModel, y0, t) # Bättre modell
# Vi kommer att vara intresserade av y1[0] respektive y2[0] som
# innehåller värden för θ i respektive modell.

4.2.3 Lösning

In [14]:
# Förenklad modell som utgår från att pendeln beskriver
# en harmonisk svängning.
def simpleModel(Y, t):
    θ, ω = Y
    retVal = [ω, - g/r*θ]
    return(retVal)

# Någorlunda korrekt modell med hänsyn taget
# till luftmotstånd och vinkel som gör att
# svängningen inte blir harmonisk.
def properModel(Y, t):
    θ, ω = Y
    retVal = [ω, -k*ω*np.abs(ω) - g/r*np.sin(θ)]
    return(retVal)

def deg2rad(angle):
    return 2* np.pi * angle/360

def rad2deg(angle):
    return angle/(2*np.pi)*360

k = 0.02
r = 0.2487435
g = 9.82
ω0 = 0
t = np.linspace(0, 40, 1000)

θ0 = (deg2rad(5))
y0 = [θ0, ω0]
y1 = odeint(simpleModel, y0, t) # Korrekt modell
y2 = odeint(properModel, y0, t) # Förenklad modell

θ0 = (deg2rad(25))
y0 = [θ0, ω0]
y3 = odeint(simpleModel, y0, t) # Korrekt modell
y4 = odeint(properModel, y0, t) # Förenklad modell

θ0 = (deg2rad(45))
y0 = [θ0, ω0]
y5 = odeint(simpleModel, y0, t) # Korrekt modell
y6 = odeint(properModel, y0, t) # Förenklad modell
In [15]:
fig, ax = plt.subplots(3, 2)
plt.suptitle("Matematisk pendel, längd: 24,87 cm. Förenklad vs ickelinjär modell.", y = 0.98)

# Kolumn 1
ax1 = plt.subplot(321)
plt.title("↓  $0 \leq t \leq 4$ sekunder ↓", y = 1.04)
plt.plot(t[0:100], rad2deg(y1[:,0][0:100]), label = "Förenklad pendel-\nmodell");
plt.plot(t[0:100], rad2deg(y2[:,0][0:100]), label = "Modell som tar hänsyn\ntill vinkel och luftmotstånd");
plt.text(2, -4.5, "$\Theta_0=5^{\circ}$", fontsize=20, horizontalalignment='center');

plt.legend(loc='upper center', bbox_to_anchor=(0.9, 1.05),
          fancybox=True, shadow=False, ncol=1, framealpha=0.5);

plt.setp(ax1.get_xticklabels(), visible=False)
plt.setp(ax1.get_yticklabels(), visible=True)
plt.grid()
plt.ylabel("Vinkelutslag [°]")

ax2 = plt.subplot(323)

plt.plot(t[0:100], rad2deg(y3[:,0][0:100]));
plt.plot(t[0:100], rad2deg(y4[:,0][0:100]));
plt.text(2, -20, "$\Theta_0=25^{\circ}$", fontsize=20, horizontalalignment='center');

plt.ylabel("Vinkelutslag [°]")
plt.setp(ax2.get_xticklabels(), visible=False)
plt.setp(ax2.get_yticklabels(), visible=True)
plt.grid()

ax3 = plt.subplot(325)
plt.xlabel("Tid [s]")
plt.plot(t[0:100], rad2deg(y5[:,0][0:100]));
plt.plot(t[0:100], rad2deg(y6[:,0][0:100]));
plt.text(2, -35, "$\Theta_0=45^{\circ}$", fontsize=20, horizontalalignment='center');

plt.ylabel("Vinkelutslag [°]")
plt.setp(ax3.get_xticklabels(), visible=True)
plt.setp(ax3.get_yticklabels(), visible=True)
plt.grid()

# Kolumn 2
ax4 = plt.subplot(322, sharey=ax1)
plt.title("↓ $28 \leq t \leq 32$ sekunder ↓", y = 1.04)

plt.plot(t[700:800], rad2deg(y1[:,0][700:800]), label = "Förenklad pendel-\nmodell");
plt.plot(t[700:800], rad2deg(y2[:,0][700:800]), label = "Modell som tar hänsyn\ntill vinkel och luftmotstånd");
plt.text(30, -4.5, "$\Theta_0=5^{\circ}$", fontsize=20, horizontalalignment='center');

plt.legend(loc='upper center', bbox_to_anchor=(1.01, 1.05),
          fancybox=True, shadow=False, ncol=1, framealpha=0.5);
plt.setp(ax4.get_xticklabels(), visible=False)
plt.setp(ax4.get_yticklabels(), visible=False)
plt.grid()

ax5 = plt.subplot(324, sharey=ax2)
plt.plot(t[700:800], rad2deg(y3[:,0][700:800]));
plt.plot(t[700:800], rad2deg(y4[:,0][700:800]));
plt.text(30, -20, "$\Theta_0=25^{\circ}$", fontsize=20, horizontalalignment='center');

plt.setp(ax5.get_xticklabels(), visible=False)
plt.setp(ax5.get_yticklabels(), visible=False)
plt.grid()

ax6 = plt.subplot(326, sharey=ax3)
plt.xlabel("Tid [s]")
plt.plot(t[700:800], rad2deg(y5[:,0][700:800]));
plt.plot(t[700:800], rad2deg(y6[:,0][700:800]));
plt.text(30, -35, "$\Theta_0=45^{\circ}$", fontsize=20, horizontalalignment='center');

plt.setp(ax6.get_xticklabels(), visible=True)
plt.setp(ax6.get_yticklabels(), visible=False)
plt.grid()

4.2.4 Kommentar

Vi ser hur svängningen utvecklas, observera tidsaxeln i respektive kolumn! Små vinklar ger ett förhållandevis litet fel även när tiden löper, medan såväl periodtidsfelet som amplitudfelet ökar med ökad startvinkel och ökad tid.

5 Referenser

Detta webbdokument är gjort med Jupyter; ett verktyg som integrerar kod, diagram, beräkningar, text och bild i ett och samma dokument. Många idéer och mycket information om hur man går till väga med numerisk lösning av differentialekvationer och ritandet av grafer har jag fått från

In [ ]: