3.4. Oefenopgaven

3.4.1. Differentiëren

Differentiëren met lijsten: Probeer het nu zelf!

In de sectie Differentiëren vind je een stukje code waarin de numerieke afgeleide \(\frac{\partial{y}}{\partial{x}}\) gegeven twee sets waarden y en x.

Voer de volgende opdrachten uit

1. Ga na of je begrijpt wat er in de functie afgeleide(x, y) gebeurt. Zie je hier de elementen van de wiskundige afgeleide terug?

2. De lijsten die teruggegeven worden door afgeleide(x, y) zijn één element korter dan de lijsten die we erin stoppen. Zie je waarom dit zo is?

3. Bereken zelf (met de hand) de afgeleide van de functie die voor y_waarden gebruikt wordt en plot deze functie in de grafiek erbij. Komt jouw antwoord overeen met het numerieke antwoord?

4. Verander het aantal elementen in de lijst x_waarden naar een (veel) kleiner getal. Wat gebeurt er met de afwijking tussen de numerieke afgeleide en de afgeleide die je zelf berekend hebt?

Gelukt?

Als het goed is heb je gezien dat de afwijking tussen de numerieke afgeleide en de afgeleide die je zelf analytisch gedaan hebt groter wordt wanneer je minder elementen in je lijst hebt. Dit is een direct effect van de discrete getallenlijn waarop we met de computer rekenen. Hoe groter het aantal elementen in de lijst x_waarden, hoe kleiner de stapgrootte \(h\) wordt en hoe dichter we dus bij de daadwerkelijke limiet kunnen komen.

Differentiëren met NumPy: Controleer het zelf

We zeggen in de sectie Differentiëren met NumPy wel dat afgeleiden berekenen met numpy nauwkeuriger is, maar is dat wel echt zo?

1. Plot in deze code ook de analytisch bepaalde afgeleide, als het goed is komt deze overeen met het resultaat van np.gradient

2. Verander het aantal stappen van x_waarden. Welk verschil valt nu op?

3. Bedenk zelf een functie om de afgeleide van te nemen en kijk of het daarvoor ook werkt. Kies hiervoor iets anders dan een polynoom, probeer bijvoorbeeld eens een exponentiële of goniometrische functie.

Gelukt?

Bij deze functie zien we dat de afgeleide een stuk betrouwbaarder is, ook bij kleine stapgroottes. Door de methode waarop numpy.gradient de afgeleide berekent zullen de afwijkingen niet over het hele bereik voorkomen, zoals eerder, maar zullen deze juist bij het begin en eind van het bereik zich laten zien.

Daarnaast zal een numerieke afgeleide altijd problemen hebben wanneer een functie snel verandert. Probeer bijvoorbeeld maar eens een afgeleide van \(\tan(x)\) te nemen.

3.4.2. Integreren

Riemann Som: Probeer het zelf!

Gebruik de code uit de uitleg om te onderzoeken hoe je het resultaat van de Riemann Som dichter bij de daadwerkelijke integraal kan krijgen.

1. Pas het aantal rechthoeken aan (dit kan door het aantal waarden in x_sample te veranderen). Wat voor effect heeft dit op de Riemann Som?

2. In het voorbeeld wordt de zogenaamde ‘linker Riemann Som’ gebruikt. Dit betekent dat de linker bovenhoek van de rechthoek op de grafiek van \(y(x)\) gelegd wordt. Herschrijf de code zodat je de ‘rechter Rieman Som’ gebruikt. Bij deze som ligt de rechterbovenhoek op de grafiek van \(y(x)\).

   Heeft dit effect op het resultaat?

3. Een andere, en vaak accuratere, versie van de Riemann Som legt het midden van de rechthoek op de grafiek van \(f(x)\). Dat betekent dus dat de hoogte van de rechthoek gelijk is aan \(f(x_i + (x_{i + 1} - x_i) / 2)\). Pas dit aan in de code en vergelijk het resultaat met de eerdere twee varianten.

Integreren met numpy: De Trapeziumregel

Laten we eens onderzoeken hoe numpy.trapz() de integraal bepaald. Dit gebeurt volgens de trapeziumregel. Pas de code van het eerste voorbeeld aan zodat je in een figuur kan zien hoe deze methode werkt.

• Gebruik hiervoor in plaats van Rectangles de klasse Polygon uit matplotlib.patches.

• Als hoekpunten van de Polygon kies je \((x_i, 0), (x_{i + 1}, 0), (x_{i + 1}, y(x_{i + 1}))\) en \((x_i, y(x_i))\).

Kan je aan de hand van de resulterende figuur uitleggen waarom deze methode in de meeste gevallen een beter resultaat oplevert dan de Riemann Som die we eerder gebruikten?

3.4.3. Differentiaalvergelijkingen

Differentiaalvergelijkingen: Zwaartekracht of Gravitatie?

Dicht bij de aarde nemen we aan dat de zwaartekracht constant is met een constante \(g = 9.81\frac{\text{m}}{\text{s}^2}\). Voor een raket zal de afstand tot de aarde echter al snel zo groot worden dat deze aanname niet meer geldt. Daarom moeten we een andere formulering van de zwaartekracht gebruiken:

\[ F_z = \frac{GMm}{r^2}\]

Waarbij \(G\) de gravitatieconstante is, \(M\) de massa van de aarde, \(m\) de massa van de raket en \(r\) de afstand vanaf het midden van de aarde tot aan de raket.

• Neem de code uit Differentiaalvergelijkingen over en pas deze aan zodat deze formulering van de zwaartekracht daarin gebruikt wordt. Hiervoor moet je enkele waarden op internet opzoeken en de functie F_z aanpassen zodat deze de hoogte van de raket als input argument krijgt.

Zie je een verschil met de eerdere resultaten?

3.4.4. Afsluitende Opdracht

Je kan je antwoorden op deze vragen controleren in Vocareum. Kijk hiervoor op Brightspace.

Beweging van een karretje

In deze opdracht ga je een simpele simulatie waarmee de snelheid van een karretje berekend wordt schrijven.

Het verloop van de snelheid van het karretje kunnen we modelleren aan de hand van de volgende bekende vergelijking:

\[ \sum F = m\cdot a \]

Op het karretje werken drie krachten:

• een constante aandrijfkracht \(F_a = 1{,}5\text{N}\)

• een kracht als gevolg van de luchtweerstand: \(F_\text{wind} = \pm k_w(v - v_w)^2\)

• een kracht als gevolg van de wrijving: \(F_\text{wrijving} = -k_c\cdot v\). Hierbij is \(k_c\) een wrijvingsconstante

In de tweede formule is \(v\) de snelheid van het karretje en \(v_w\) de windsnelheid in de richting van het karretje. Normaal gesproken zal \(F_w\) een negatieve waarde hebben, maar het kan voorkomen dat het karretje zoveel wind mee heeft dat deze juist een positieve waarde aanneemt.

Om de simulatie op te bouwen moeten een aantal stappen gezet worden:

1. Het importeren van de benodigde modules

2. Het inlezen van de gegevens voor \(v_\text{wind}\)

3. Het implementeren van functies voor de snelheidsafhankelijke krachten

4. Het schrijven van een for-loop om iteratief de oplossing te bepalen

5. Het plotten van het resultaat

Stap 1: Het importeren van de modules

Begin je script door pandas, numpy en matplotlib te importeren

Stap 2: Lees de data voor \(v_\text{wind}\) in

Het bestand met de gegevens over de windsnelheid in de richting van het karretje is hier (rechts-klik, Opslaan als…) te vinden.

Schrijf een stukje code dat deze data inleest en vervolgens de inhoud van de kolom met windsnelheid tegen de tijd plot.

Voorbeeld-resultaat

Stap 3: Functies voor de weerstandskrachten

De weerstandskrachten (als gevolg van de luchtweerstand en de wrijving) kan je het makkelijkst met behulp van functies implementeren in de code. Schrijf daarvoor twee functies:

• Een functie F_weerstand, met als inputs de huidige snelheid van het karretje en de huidige windsnelheid. De functie moet als output de luchtweerstandskracht teruggeven. Let erop dat de kracht positief moet zijn als de windsnelheid groter is dan de huidige snelheid van het karretje, en negatief als het karretje tegenwind heeft (of harder gaat dan de meewind). Gebruik hiervoor de volgende formule:

  \[ F_\text{wind} = \pm k_w (v - v_w)^2,\quad k_w = 0.5 \]

• Een functie F_wrijving, met als input de huidige snelheid van het karretje. Deze functie moet de wrijvingskracht teruggeven.

  Hierbij hoort de volgende formule:

  \[ F_\text{wrijving} = -k_c \cdot v,\quad k_c = 1 \]

Stap 4: Iteratie/Integratie

Om een numerieke benadering voaor de oplossing van de differentiaalvergelijking te bepalen maken we gebruik van de volgende recursieve formule, waarbij \(i\) een index aangeeft:

\[v[i + 1] = v[i] + \frac{1}{m}\left(\sum F[i]\right)\cdot \Delta t \]

Bestudeer het voorbeeld uit Differentiaalvergelijkingen en pas die code aan zodat de beweging voor het karretje gesimuleerd wordt. Gebruik hiervoor de volgende gegevens:

Grootheid	Waarde
Massa \(m\)	\(10\,\text{kg}\)
Tijdstap \(\Delta t\)	\(0{,}1\,\text{s}\)

Simuleer de beweging gedurende de 100 seconden waarvoor we wind-data heben. Sla de positie en snelheid van het karretje op in twee numpy arrays met de namen positie en snelheid.

Stap 5: Plot het resultaat

Plot tenslotte het resultaat in een mooie grafiek. Maak hiervoor een figuur aan met twee assenstelsels naast elkaar. Noem de figuur fig en de assenstelsels gezamenlijk ax.

Plot in de linker figuur de positie van het karretje tegen de tijd en in de rechter figuur de snelheid van het karretje tegen de tijd. Voorzie zoals altijd je assenstelsels van nuttige labels.

Subplots naast elkaar

Met het commando plt.subplots() geven de eerste twee getallen aan hoeveel assenstelsels er in de figuur moeten komen. In de documentatie kan je zien hoe je ervoor zorgt dat er twee assenstelsels naast elkaar komen te staan. Onderaan die pagina staan ook verschillende voorbeelden.