Geographical Research Methods 3: Numerical Modelling

2021-2022

Q1: simulation of a water wave. The baselevel is 5 meters, the equation of the wave was given (sin-wave) and it was located on the righthand side of the domain, the domain: 0-20 km, the wave is located between 19-20 km. The movement of the wave needed to be simulated, a negative advection coefficient is given (wave moved towards the left). Simulate the advection of the wave with the leapfrog method and plot the wave at 0,6,12,18,24 hours, without (periodic) boundary conditions ! so y(1) and y(end) = 5.

Q2: simulate the propegation of the wave now using both advection and diffusion, again plot the wave at 0,6,12,18,24 hours.

Q3: Plot the height of the water level at the 10 km mark during the whole 24 hours.

extra questions if you had the time:

Q4: what was the maximal height the wave reached at the 10 km point during the whole 24 hours.

Q5: ?

Questions on the assignment:

• encapsulations (private, protected, public), what do they mean and what do they do ?

• handle and mixin.heterogeneous, what do they do?

• what are the advantages of OOP and for which other geogoraphic topics could they be used for ?

2020-2021

Part I:

1) Temperature over a linear distance is represented by half of a sine wave with on y-axis the temperature and on x-axis the distance (total distance of 10 km). The temperature is constant and equal to 20°C until at 5 km it starts increasing; at 5.5 km the temperature reaches a maximum of 25°C; and at 6 km it is 20°C again (see fig). The city of Leuven lies between km 5 and km 6 of this linear distance and this wave represents the urban heat island effect. Model advection of temperature (using the leapfrog scheme) over 24 hours, with a wind speed from left to right of 0.3 m/s. You also have to implement (periodic) boundary conditions and the Robert-Asselin filter. (During the exam, you have access to the sript initial.m from the exercises on Toledo which you can use as a starting template).

2) Use the same script but also add a diffusion term with D = 2.5 m²/s. Determine the moment in time when the amplitude of the wave has decreased to 2.5°C (half of the initial amplitude).

3) Use the previous script but adapt it to calculate the temperature evolution for a range in D values. D variates between 1.5 and 10 in steps of 0.5. Calculate the amplitude for each D and then plotting diffusivity on the x-axis and amplitude on the y-axis.

4) This was an extra question if you had time left.

Part II: questions about the OOP assignment (answering in words)

• Motivate why you used encapsulations.

• Explain why it was necessary to use 'handle' in your class script

• Give example of application when making use of OOP is more helpful/efficient.

9/06/2020

Intro: Deepwater Horizon oil spill: environmental disaster, Gulf of Mexico [2010]

“The Deepwater Horizon rig, owned and operated by offshore-oil-drilling company Transocean and leased by oil company BP, was situated in the Macondo oil prospect in the Mississippi Canyon, a valley in the continental shelf. The oil well over which it was positioned was located on the seabed 4,993 feet below the surface and extended approximately 18,000 feet into the rock. On the night of April 20 a surge of natural gas blasted through a concrete core recently installed by contractor Halliburton in order to seal the well for later use. It later emerged through documents released by Wikileaks that a similar incident had occurred on a BP-owned rig in the Caspian Sea in September 2008. Both cores were likely too weak to withstand the pressure because they were composed of a concrete mixture that used nitrogen gas to accelerate curing.”

Goal

Your goal for today will be to simulate the amount of oil pollution (ppt) in the Gulf of Mexico following the explosion at the Deepwater Horizon oil rig. Oil will be spread over the Gulf due to ocean currents and dilution in the sea water. You have to simulate both processes in each of the following exercises.

1) 1D simulation

In a first exercise, you will simulate oil pollution over a transect between New Orleans and Havana (Fig. 2). Consider a one-dimensional case with a domain length of 1000 km. The Deepwater Horizon rig is situated at 100 km from New Orleans. In the next couple of steps, you will calculate the time required to reach a specific oil concentration in the center of the Gulf Of Mexico (1 × 10^5 ppt), 500 km from New Orleans, and 400 km from the rig. You can assume a constant oil concentration of 1 × 10^9 ppt at the Deepwater Horizon rig.

- Write the differential equation you must solve for this problem in full (one equation) in the box below:

Write a script which calculates the time required to obtain a critical oil concentration of 1 × 10^5 ppt at at the center of the transect. Write the script such that the direction of the ocean current can be adjusted from left to right and vice versa, only by adjusting one variable.

- Set the initial conditions: at Deepwater Horizon the oil concentration equals 1 × 10^9 ppt, over the rest of the domain the concentration is 0 ppm.

- The ocean current equals 0.25 km h-1 (left/right)

- The diffusion constant equals 15 km2 h-1

- Use a spatial resolution of 2 km

- Calculate the maximum timestep guaranteeing a stable solution

- For the boundary conditions you can set the outer part of the domain equal to the value of the inwards neighboring cell:

OIL(1) = OIL(2);

OIL(end) = OIL(end-1);

- Throughout the simulation, the concentration at Deepwater Horizon remains 1 × 10^9 ppt, set the concentration at Deepwater Horizon to this value during each iteration (before calculating transport and dispersion of the particles).

- Plot the intermediate result every n_out iterations, with n_out having a value of your choice

- Plot the concentration (y-axes) at the center of the Gulf (500 km) over time (x-axes)

- Plot the final concentration (y-axes) over the distance of the transect (x-axes)

- Make a single ‘.m’ file with the entire solution, called ‘Name_DH_1D.m’. You do not have to save the figures separately.

HINT:

- Update the solution both after calculating the diffusion and after calculating advection to guarantee a stable solution at CFL 0.9

2016-2017

Reeks 1

2015-2016

Exam matlab 19 January 2015 en 12 January

2014-2015

Exam matlab 16 January 2015

Hurricane Katrina was the deadliest and most destructive Atlantic tropical cyclone of the 2005 Atlantic hurricane season. Around 1,833 people died in the hurricane and subsequent floods; total property damage was estimated at 108 billion USD. In this exercise we are going to get insight in how this event affected the local economy of Louisiana and in general how economies evolve. For simplicity and illustrative purposes we will assume that there is no interaction with the rest of the world and Lousiana act as a local economy.

The income (Y ) can be written as a function of capital (K), labour (L) and technology (A) as follows:

Y = Kα · (A · L) 1−α

The constant α is called the elasticity of the model and can have a value in between 0 and 1.

The initial values (for K, A, L, and α in 1950) are:

K=500

A=100

L=300

α = 0.3

Population growth affects a similar growth in labour L grows with a rate of 1.2% following:

dL dt = r · L

where r is the growth rate per year (r = 0.012).

exercise 1: Plot the income as a function of time for the period 1950 until 2049. Note that you can refer to 1950 as ’year=1’ and 2049 as ’year=100’, but make sure the axis of your plot is labelled correctly from 1950.

Now the disaster kicks in. Due to people deceasing and moving away from Louisiana, population was reduced by 7% in 2005. In addition, half of the capital was destroyed or incapacitated by Katrina.

exercise 2: Modify your program from the previous exercise so that it takes into account this reduction in labour and capital in the year 2005.

Now we are going to take into account the change in capital due to savings (fraction of output that is devoted to investment; first term in equation below) and the depreciation (the decrease in value of assets; second term in the equation below).

dK dt = s · Y − δ · K

where s is the fraction of output that is devoted to investment (s = 0.1) and δ the depreciation rate (δ = 0.05).

exercise 3: Modify your program from the previous exercise so that it takes into account this change in capital. Instead of income (Y), now plot the income per capita, which is given by:

y = Y/L

Also display the time at which income per capita has reached exactly the average value between 1950 and 2049. In addition, display the time at which the increase in income per capita is largest.

exercise 4: Now we are going to study the effect of technology. To take this into account, we first have to calculate the growth rate of the income per capita over the period 1950-2049. You can do this by using linear regression. Secondly you have to write a function that calculates the growth rate of y as a function of A. In a last step you will vary A between 0 and 200 with steps of 10. Plot now the growth rate of y as a function of technology A.

exercise 5: Now we are simultaneous going to take into account the effect of the magnitude of the disaster. For this, you will vary the reduction of capital by the disaster between 0% and 99%. Modify your program from the previous exercise so that it shows the growth rate of y as a function of both technology A and percentage of capital destroyed by the disaster (3D plot).

You can save the programmes in a folder Katrina_givenname_familyname.m

Exam matlab 20 January 2015

In this exercise we are going to use a very simple model of the earth (land only) on which grass grows. The net rate of change in the fraction of the area of the earth covered by grass (A) is given by:

dA/dt = A ((1 − A) · Gf − D)

where D is the death rate (a constant of 0.1 per 10 millon years).

The growth rate of the grass is 0.4 per 10 millon years (Gf ).

exercise 1: Discretise this equation (A(i + 1) = ...). Use the discretised equation to calculate A as a function of time. Every time step in the program corresponds to a period of 10 millon years. Run the model for 200 time units (is 2 billion years). Use a starting value for A of 0.001. Write to the screen the time at which the growth stabilizes (here defined as that the change over one timestep becomes smaller than 1% of the difference between A at the time under consideration compared to the initial value of A.).

exercise 2: Now we are going to modify the growth and death rate. Vary both Gf and D between 0 and 1 and make a 2.5D plot of the value of A after 2 billion years as a function of both Gf and D.

Now we are going to extend the model to include the temperature on earth and how growing conditions are affected by this temperature. First we are going to have a linear increase in temperature from a value of 270 K at the start of the model run to a value of 290 K at the end of the model run. The growth rate of the grass per 10 millon years (Gf ) is a function of the temperature: it only grow in-between a temperature range of 273.15 and 283.15 and at a temperature of 278.15 K, the growth rate is maximal:

Gf = 1 − 0.04 · (278.15 − Te) ² for 273.15 ≤ Te ≤ 283.15

Gf = 0 for Te < 273.15 and Te > 283.15

exercise 3: Make a script that displays the fraction of the earth covered by grass (A), taking into account the temperature increase and varying growth rate. Again run the model for 200 time units (is 2 billion years) and use a starting value for A of 0.001.

In a last step, we are also taking into account the change in the reflectivity of the earth due to the growth of grass and how this affects temperature. The equilibrium temperature (Te) of the earth is given by:

Te = (S0 · (1 − αplanet))/ (4 · σ))^ 0.25

where S0 is the solar constant (W m−2 ) and σ the StefanBoltzmann constant.

The planetary albedo of the earth (αplanet) is a function of the fraction of the area of the earth (A) covered by grass, the albedo of the grass (αgrass) and the albedo of the bare surface (αbare):

αplanet = (1 − A) · αbare + A · αgrass

exercise 4: The solar constant S0 increases linearly with time: At the beginning of, the sun provides an S0 of 2600 W m−2 and at the end of (after 2 billion uears) S0 is 3500 W m−2 . Mare a script that produces one graph in which the evolution of the fraction of the area covered by grass is given as a function of time and one graph in which the evolution of the equilibrium temperature of the earth is given as a function of time for a period of 2 billion years (corresponding with 200 time steps). Write the time when the area covered with grass is at its maximal value to the screen. Do the same for the time when the area covered with daisies exceeds 30% of the total area of the planet.

exercise 5: Now let the albedo of grass increase from 0.25 to 0.6 with steps of 0.01. Plot the time when the area of daisies reaches its maximum value as a function of the albedo of the daisies.

You can save the program under grass_givenname_familyname.m

σ = 5.67 · 10^−8 W m^−2 K^−4

αbare = 0.5

αgrass = 0.25

D = 0.1 per 10 millon years

2013 - 2014

22 januari 2014:

Modelling a crowd in a corridor

Modelling a crowd in a corridor has been used to study optimal escape route in case of emergency (e.g. fires). In this exercise we are going to model te behavior of one and more people in a corridor with a length of 300 m. For simplicity reasons, we will model this model as a one dimensional line.

1. As humans tend to walk at about 1.4 m/s, we use this speed (v0) to calculate the position of one person in a corridor. The man starts at position 0 and we use the following equation to calculate the position (x) as a function of time: dx/dt = v0. We do all model runs using a timestep dt of 1 s a,d we model in total the position of the man for 5 minutes. When te man leaves the corridor, the position of the man will obtain a NaN value, which means Not-a-Number (x=NaN). Make a script that calculates the position of the man as a function of time. The script should generate a plot that shows the position of the man as a function of time.

2. Now we are going to model the behavior of more people in the corridor. The exercise is the same as the previous one, but now at every timestep one person enters the corridat and starts at position 0. Make a script that calculates the position of each person as a function of the time. The script should generate a 2D or 2.5D plot that shows the position as a function of time (x-axis) and the number of the person (y-axis).

3. The crowd gets larger and this exercise is the same as the previous exercise except that now every timestep 10 people enter the corridor. The speed of the people varies from person to person with a minimum of 0 and a maximum of 2.5 m/s. You can get this into your model by using the function 'rand', which generates random values between 0 and 1. You will have to find a way to use this function for generating values from the uniform distribution on the interval [0-2.5]. Note that for one person the speed does not vary with time. Make a script that calculates and generates a 2D or 2.5D plot that shows the position as a function of time (x-axis) and the number of the person (y-axis).

4. Make a script that calculates the position of the slowest person and write this to the screen. Also write to the screen the number of people that have left the corridor and the moment at which the fastest person has left the corridor. In addition, calculate the average speed of all people that are in the corridor at the end of the simulation.

Challenge:

5. Now the speed is going to vary as a function of the person, but also as a function of the number of people nearby. I.e. persons will walk slower if there are more people nearby. To take this behavior into account, we need the people density for each pedestrian, which we define here as the total number of people from one meter behind to one meter ahead of the pedestrian. If a person is closer than one meter from the entrance or exit of the corridor, you take the first or the last two meters of the corridor to calculate this density. Each pedestrian has a maximum speed (vmax), which varies from person to person with a minimum of 0 and a maximum of 2.5 m/s. Again, you can get this into your model by generating values from the uniform distribution on the interval [0-2.5] (see part 3). If the density is smaller or equal to 10 people than the speed of the pedestrian is vmax. If the density is large or esual to 100 people, than the speed drops to zero. In the density range from 10 to 100 people. The speed (v) is given bij v = - ((vmax*density)/90) - ((vmax*100)/90). Make a script that calculates this and generates a 2D or 2.5D plot that shows the position as a function of time (x-axis) and the number of the pedestrian (y-axis). In addition a plot should be generated that shows the density of people a a function of time (x-axis) and the number of pedestrian (y-axis).

6. Now the exit gets blocked and only 8 people can leave the corridor in one timestep. The others stay at their position. Make the same plot as in exercise 5 to study the behavior of the crowd as a function of time.

13 januari 2014:

Verkeer op een ringweg

Het modelleren van verkeersstromen is een methode die gebruikt wordt om te helpen bij het ontwerpen van transportinfrastructuur. In deze oefeningen gaan we verkeersstromen modelleren. We beginnen met een één-dimensionaal model van de auto's die bewegen langs een gesloten ringweg. U kunt dit modelleren door de weg onder te verdelen in discrete posities en nummers aan deze posities toe te wijzen van 1 tot N, waarbij N=1000. De auto gaat één positie vooruit per tijdsstap. We kunnen een ringweg modelleren met behulp van 'periodieke randvoorwaarden'. Dit betekent dat de autoop positie N devolgende tijdsstap zal verplaatsen naar positie 1.

Laten we beginnen met één auto op de ringweg!

1. Elke tijdsstap gaat de auto een positie verder op deze ringweg. Maak een script dat de positie van de auto berekent als functie van de tijd, met de periodieke randvoorwaarden. Het script moet een plot produceren die de positie van de auto geeft als functie van de tijd.

2. Nu zijn er verkeerslichten op positie 10, 70, 250 en 500 en kan de auto alleen vooruitgaan op tijdstippen die u kunt delen door 10 (dus 10, 20, 30, 40, ...). Maak een script dat een plot genereert, die de positie geeft van de auto als functie van de tijd.

In een tweede reeks oefeningen zijn er meer auto's op de weg. Er zijn in totaal M auto's. De dichtheid p van de auto's wordt gedefinieerd als M/N.

3. We beginnen met een dichtheid van 0,5, dat wil zeggen 500 auto's op de weg. Al deze auto's moeten een initiele positie krijgen. Maak een script dat elke auto een willekeurige beginpositie geeft. Gebruik de functie randi( imax, [M,N]); deze geeft een m bij n matrix met willekeurige gehele getallen getrokken uit de discrete verdeling met interval [1, imax]. Zorg er wel voor dat er zich maximaal maar 1 auto op elke positie bevindt. Maak een script dat een scatterplot van de uitgangspositie als functie van het autonummer creëert.

4. De auto zal alleen verplaatsen wanneer de weg niet geblokkeerd is, dus als er geen auto aanwezig is op de voorliggende positie. Maak een script dat de positie van auto nummer 1 en nummer 10 als functie van de tijd berekend. Het script moet een plot genereren dat de positie van deze twee auto's toont als functie van de tijd voor p= 0,5

5. Bereken en plot de gemiddelde snelheid van de auto's als functie van de tijd voor p = 0,5. De snelheid ligt tussen de waarden 0 en 1 posities per tijdstap.

6. Nu gaan we het gedrag van de auto's als functie van de dichtheid bestuderen. Bereken en plot (2D of 2.5D plot) de gemiddelde snelheid van de auto als functie van de dichtheid (op de x-as) en de tijd (op de y-as).

+ 2 extra oefeningen waarbij je moet simuleren voor een 2-baans ringweg.

8 januari 2014

Het klimaatsysteem van onze aarde is comples, en daarom is het moeilijk de evolutie van het klimaatysteem te doorgronden.Door gebruik te maken van vereenvoudigde modellen kunnen we het begrip van het klimaatsysteem verbeteren. Hier werken we met een vereenvoudigde versie van het zgn. "daisyworld" model bedacht door James Lovelock: we beschouwen een planeet zonder atmosfeer waar alleen licht en donkergrijze madeliefjes (daisies) kunnen groeien. Deze bloemen groeien slechts binnen een goed gedefinieerd temperatuursinterval. Om de evolutie van de temperatuur in daisyworld te kunnen simuleren gebruiken we de volgende formules:

De evenwichtstemperatuur (Te) van de planeet wordt gegeven door:

Te = ((So*(1-alfaplanet))/(4*sigma))^0.25

waar So de zonneconstante (W m^-2) en sigma de StefanBoltzmann constante is. Het planetair albedo van daisyworld (alfaplanet) is een functie van de fractie van het oppervlak van de planeet (A) bedekt met madeliefjes, het albedo van de madeliefjes (alfadaisy) en het albedo van het kale oppervlak (alfabare):

alfaplanet = (1-A)*alfabare+A*alfadaisy

De nettoverandering van A per tijdseenheid wordt gegeven door:

dA/dt = A((1-A)*Gt-D)

waar D het sterftecijfer is (een constante met een waarde van 0.1 per 10 miljoen jaar). Het groeicijfer van de madeliefjes per 10 miljoen jaar (Gf) is een functie van de temperatuur: zij kunnen enkel groeien in een temperatuursinterval tussen de 273.15 en 283.15 K en bij een temperatuur van 278.15 K is het groeicijfer maximaal:

Gf = 1-0.04*(278.15-Te)^2 voor 273.15<=Te<=283.15

Gf = 0 voor Te<273.15 en Te>283.15

Om de differentiaal vergelijking te kunnen oplossen, kan een eindige verschillen methode (Euler) worden gebruikt:

A(t+1) = A(t) + delta t * A(t) * ((1 - A(t)) * Gf(t) - D)

Het idee is nu om een programma te schrijven waarin de temperatuur (Te) van daisyworld en de fractie van het oppervlak bedekt met bloemen worden berekend als een functie van de tijd. Elke tijdstap in het model correspondeert met een tijdsperiode van 10 miljoen jaar. Het model moet gerund worden voor 200 tijdseenheden (2miljard jaar). De zonneconstante So neemt lineair toe met de tijd. In het begin van daisyworld levert de zon een So van 2600 W m^-2 en op het einde (na 2 miljard jaar) is So 3500 W m^-2.

Je programma zal dus twee grafieken moeten produceren: eentje waarin de evolutie van de fractie van het oppervlak bedekt met bloemen (A) wordt weergegeven als functie van de tijd en eentje waar de evolutie van de temperatuur (Te) wordt weergegeven als functie van de tijd voor een periode van 2 miljard jaar (corresponderende met 200 tijdsstappen). De waarde voor het albedo van de madeliefjes (oftewel 0.4 oftewel 0.6) moet kunnen worden ingegeven via het scherm. Noteer dat een minimale waarde voor A nodig is om de madeliefjes te kunnen laten groeien. Immers uit vgl. 3 blijkt dat de bloemen enkel kunnen groeien als A niet gelijk is aan 0. Voeg daarom een regel toe in je programma die ervoor zorgt dat A nooit kleiner wordt dan 0.001.

Een tweede opdracht is om een optie aan je programma toe te voegen die ervoor zorgt dat het tijdstip wanneer het oppervlak bedekt met bloemen maximaal is wordt weggeschreven naar het scherm.

Je kan je programma bewaren onder daisyworld1_voornaam_achternaam.m

sigma=5.67*10^-8 W m^-2 K^-4

alfabare=0.5

alfadaisy= 0.4 (voor donkergrijze madeliefjes) of 0.6 (voor lichgrijze madeliefjes)

D=0.1 per 10 miljoen jaar

Als je nog van een uitdaging houdt kan je ook het onderstaande proberen op te lossen: laat nu het albedo van de madeliefjes varieren tussen 0.4 en 0.6 met stappen van 0.01. Maak een grafiek waarin het tijdstip wanneer A maximaal is wordt weergegeven als een functie van het albedo van de bloemen. Je kan je programma bewaren onder daisyworld2_voornaam_achternaam.m

2007 – 2008

Matlab:

We gaan het hebben over hoe een bevolking in een bepaald gebied met haar hulpbronnen omspringt. We modelleren dat met differentiaalvergelijkingen.

Voor de bevolking (N) geldt dat ze exponentieel groeit, dus:

dN/dt=rh*N(1-S/m)

oftwel:

dN/dt= (bh-dh)*N(1-S/m)

waarbij b=het basisgeboortecijfer en d= het basissterftecijfer bij afwezigheid van de hulpbron: dat zijn cijfers per hoofd van de bevolking, S=subsistance rate, zijnde de hoeveelheid voedsel nodig om te overlegen, m is de actuele beschikbare voedselhoeveelheid, m wordt berekend per huishouden als:

m=e*theta*q*x

We nemen aan dat de bevolking een dierenpopulatie bejaagt, de jachtbuit wordt bepaald door de volgende vgl:

mtot=e*theta*q*x*(N/p)

Waarbij e=de hunting effort per persoon, theta=een technologiecoefficient, q=de catchability en x= aantal dieren, uitgedrukt in Animal Units(AU), p=aantal personen per huishouden

De dierenpopulatie groeit volgens een zogenaamde logistische groeifunctie:

dx/dt=(-dx+bx)*x*(1-(x/K)-mtot

Waarbij K de maximale is, x is dus de actuele grootte.

Aan jullie de volgende vraag: maak een programma dat de evolutie van bevolking N en de dierenpopulatie x kan voorspellen. Om dat te doen gebruik je voor het oplossen van de differentiaalvgl een eindige-verschillenmethode (die van Euler). Heel concreet betekent dat voor de bevolkingsgroei (D=delta):

DN/Dt=rNt

of:

N(t+1)-N(t)/t+1-t=rN(t)

of:

N(t+1)=N(t)*(1+r)

Je kan zelf de analoge vgl schrijven voor de hulpbron.

De bedoeling is dat je programma de interactie tussen menselijke populatie en haar hulpbronnen gaat beschrijven ifv de tijd door telkens de waarden voor N en x te berkenen voor tijdstip t+1 op basis van de waarden van tijdstip t. Ieder tijdstip t in het programma dient overeen te komen met 1 jaar. Hieronder vind je een tabel met de nodige beginwaarden/parameterwaarden.

N0=400

theta=0.15

rh=bh-dh=0.0443

K=6E7 AU

S=2.9 AU

bx=0.30

dx=0.05

q=1E-10

e=6.43E3

p=4

Je programma dient een grafiek te produceren die de evolutie van N en x stimuleert over een periode van 1500 jaar, ervan uitgaand dat K=x op tijdstip 0. Waarden voor de groeisnelheid van de populatie, rh dienen via het scherm ingegeven te kunnen worden zodat verschillende resultaten snel kunnen vergeleken worden. Verdfer dient je programma de maximale bevolkingswaarde die bereikt wordt en het tijdstip waarop de bevolking 97% van deze waarde bereikt te rapporteren op het scherm. Sla dit programma op en bewaar het onder animal1_jevoornaam_jeachternaam.

In een tweede stap wijzig je je programma zodanig dat je het automatisch kan laten lopen voor verschillende waarden van rx=bx-dx, variërend van 0.05 tot 0.5 in stappen van 0.05. Het programma rapporteert dan voor elke simulatie de maximale bevolking en het tijdstip vwaarop de bevolking 97% bedraagt van die maximale waarde en stockeert deze informatie in specifieke variabelen. Op het einde worden 2 plots gemaakt die de maximale bevolking en het tijdstip waarop 97% van deze maximale bevolking bereikt wordt weergeeft ifv rx. Sla dit programma op onder animal2_jevoornaam_je achternaam