PC deel examen

Bijlagen

📦 Download alle bijlages (ZIP)

# -*- coding: utf-8 -*-
"""
EXAMEN DIFFERENTIAALVERGELIJKINGEN - 1STE ZITTIJD - DEEL 2 - VRAAG 8

Naam student: ...............................

>>> PLAATS UIT COMMENTAAR WAAR NODIG EN VUL VERDER AAN <<<

"""

#%%############################################################################
# Imports
###############################################################################

import numpy as np
import sympy as sp
import matplotlib.pyplot as plt
from matplotlib import cm

import matplotlib
matplotlib.use("qtagg")
plt.ion()

#%%############################################################################
# (a) Discretiseer de randvoorwaarde (geen code nodig)
###############################################################################

#%%############################################################################
# (b) Duid aan in Figuur 1 (geen code nodig)
###############################################################################

#%%############################################################################
# (c) Benader de oplossing
###############################################################################

### Parameterwaarden
D = 0.9
v = 1.5
k = 0.006
dx = 0.5
dt = 0.2

### x- en t-intervallen
x_span = [0, 2000]
t_span = [0, 1440]

#### Arrays
# x-array
nx = int(...)     # totaal aantal knopen voor x
x = np.linspace(..., ..., ...)

# t-array
nt = int(...)     # aantal tijdstippen
t = np.linspace(..., ..., ...)

# C(x,t)-array
C = np.zeros((nx, nt))

### Dirichlet rand- en/of beginvoorwaarde (indien van toepassing)
C[..., ...] = ...
C[..., ...] = ...

### Eindige differentie
# Loop over de benodige knopen en vind de benadering (dit kan een minuut duren)
for j in range(..., ...):
    for i in range(..., ...):
        C[..., ...] = ...
    # Gediscretiseerde Neumann-randvoorwaarde (indien van toepassing)
    C[..., ...] = ...

# Indien er imaginaire knopen zijn, verwijder die dan
C = C[..., ...]
x = x[...]

#%%############################################################################
# (d) Maak het surface plot en bewaar de figuur
###############################################################################

### Creëer de juiste vorm en types van de data om te plotten
C_T = np.transpose(C)
x_mg, t_mg = np.meshgrid(..., ...)

### Plot de figuur
fig5, ax5 = plt.subplots(num=5, subplot_kw={"projection": "3d"})
surf = ax5.plot_surface(..., ..., ...,
                       cmap=cm.coolwarm,
                       linewidth=0, antialiased=False)

### Benoem de assen en figuur
ax5.set_xlabel("$x$ (meter")
ax5.set_ylabel("$t$ (minuten)")
ax5.set_zlabel("$C$ (mg/L)")
ax5.set_title("Concentratie hypochloriet in tijd en ruimte")

### Sla de figuur op
fig5.savefig("figuur_concentratie_hypochloriet.png")

#%%############################################################################
# (e) Plot de evolutie van de concentratie bij het controlepunt
###############################################################################

### Creëer de array met de evolutie van de concentratie bij het controlepunt
x_controle = 1000
index_controle = ...
C_controle = C[..., ...]

### Plot de figuur
fig6, ax6 = plt.subplots(num=6)
ax6.plot(t, C_controle)
ax6.axhline(0.1, color='k', ls='--', lw=1)

### benoem de assen en figuur
ax6.set_xlabel("$t$ (minuten)")
ax6.set_ylabel("$C$ (mg/L)")
ax6.set_title("Concentratie hypochloriet op controlepunt (1000 meter)")

### Sla de figuur op
fig6.savefig("figuur_controlepunt.png")

### Vind de tijd t' waarvoor de concentratie voor het eerst
### groter wordt dan de drempelwaarde 0.1
C_drempel = 0.1
drempel_index = ...
t_drempel = t[drempel_index]

### Print deze waarde
print(f"Het gezochte tijdstip t' is {np.round(t_drempel,1)} minuten.")

#%%############################################################################
# (f) Vind de concentratie voor x=2000, t=1440
###############################################################################

### vind en print de waarde
C_einde = C[..., ...]
print(f"De gezochte concentratie is {np.round(C_einde,3)} mg/L.")

import numpy as np
import matplotlib.pyplot as plt


# Richtingsveld (voor één 1ste orde DV)
def richtingsveld(f, tintval, yintval, *args, n=25, ax=None, color='k'):
    # Deze functie plot het richtingsveld van y' = f(t,y)
    # Syntax: richtingsveld(f, tintval, yintval, n)
    # Input: f: fuction handle naar het rechterlid van de DV
    #        grenzen van de grafiek:
    #        tintval = [tmin, tmax]
    #        yintval = [ymin, ymax]
    #        *args= extra argumenten voor de functie f
    #        n = precisie van het richtingsveld, default n = 25
    #        ax = axis-object indien verschillende plots in een
    #             figuur, default ax=None
    #        color = kleurcode pijltjes, default color='k' (zwart)
    # Output: richtingsveld van y' = f(t,y)
    tmin = tintval[0]
    tmax = tintval[1]
    ymin = yintval[0]
    ymax = yintval[1]
    T = np.linspace(tmin, tmax, n)
    y = np.linspace(ymin, ymax, n)
    # Vector field
    [X, Y] = np.meshgrid(T, y)
    U = 1
    V = f(T, Y, *args)
    if V.ndim == 1:
        m = V.shape[0]
        V = np.reshape(V, (1, m))
        V = np.repeat(V, repeats=m, axis=0)
    # Normalize arrows
    N = np.sqrt(U**2 + V**2)  
    U = U/N
    V = V/N
    if ax:
        ax.quiver(X, Y, U, V, angles='xy', color=color)    
        ax.set_xlim(tintval)
        ax.set_ylim(yintval)
        ax.set_xlabel('$t$')
        ax.set_ylabel('$y$')
    else:
        plt.quiver(X, Y, U, V, angles='xy')    
        plt.xlim(tintval)
        plt.ylim(yintval)
        plt.xlabel('$t$')
        plt.ylabel('$y$')
        plt.show()


# Methode van Euler voor één 1ste orde DV
def mijnEuler(f, interval, y0, n, *args):
    # Benaderen van y' = f(t, y) over tintval = [t0, tn]
    # met beginwaarde y(t0) = y0 en stapgrootte dt = (tn - t0)/n
    # Syntax: [t, y] = mijnEuler(f, tintval, y0, n)
    # Input:  f: fuction handle naar het rechterlid van de DV
    #         tintval = [tmin, tmax]: grenzen van de grafiek
    #         y0 = beginwaarde
    #         n = aantal stappen
    #         *args= extra argumenten voor de functie f
    # Output: t = tijdsvector
    #         y = benadering van y' = f(t,y)
    y = np.zeros(n+1)      # pre-allocatie van y als array
    y[0] = y0                 # y0 is het eerste element van
    dt = (interval[1] - interval[0])/n
    t = np.linspace(interval[0], interval[1], n+1, endpoint = True)
    for k in range(0, n):     # k = 0, 1, 2, ..., n-1
        y[k+1] = y[k] + dt*f(t[k], y[k], *args)
    return [t, y]


# # Methode van Euler voor één 1ste orde DV
# def mijnEuler(f, interval, y0, n, *args):
#     # Benaderen van y' = f(t, y) over tintval = [t0, tn]
#     # met beginwaarde y(t0) = y0 en stapgrootte dt = (tn - t0)/n
#     # Syntax: [t, y] = mijnEuler(f, tintval, y0, n)
#     # Input:  f: fuction handle naar het rechterlid van de DV
#     #         tintval = [tmin, tmax]: grenzen van de grafiek
#     #         y0 = beginwaarde
#     #         n = aantal stappen
#     #         *args= extra argumenten voor de functie f
#     # Output: t = tijdsvector
#     #         y = benadering van y' = f(t,y)
#     y = np.zeros(n+1)     # pre-allocatie van y als array
#     y[0] = y0                # y0 is het eerste element van
#                              # de oplossingsarray y
#     dt = (interval[1] - interval[0])/n
#     t = ...                  # maak hier de array t aan
#     for k in ...:            # k = 0, 1, 2, ..., n-1
#         # schrijf hier een for-lus die de waarden y[k+1] berekent
#         # op basis van t[k] en y[k]
#         ...
#     return [t, y]


# Midpoint methode voor één 1ste orde DV
def mijnMidpoint(f, interval, y0, n, *args):
    # Benaderen van y' = f(t, y) over tintval = [t0, tn]
    # met beginwaarde y(t0) = y0 en stapgrootte dt = (tn - t0)/n
    # Syntax: [t, y] = mijnMidpoint(f, tintval, y0, n)
    # Input:  f: fuction handle naar het rechterlid van de DV
    #         tintval = [tmin, tmax]: grenzen van de grafiek
    #         y0 = beginwaarde
    #         n = aantal stappen
    #         *args= extra argumenten voor de functie f
    # Output: t = tijdsvector
    #         y = benadering van y' = f(t,y)
    y = np.zeros(n+1)
    y[0] = y0
    dt = (interval[1] - interval[0])/n
    t = np.linspace(interval[0], interval[1], n+1, endpoint = True)
    for k in range(0, n): # k = 0, 1, 2, ..., n-1
        # >>> alpha = 1 <<<
        # m1 = f(t[k], y[k])
        # m2 = f(t[k] + dt, y[k] + dt*m1)
        # y[k+1] = y[k] + dt/2*(m1 + m2)
        # >>> alpha = 1/2 <<<
        m1 = f(t[k], y[k], *args)
        m2 = f(t[k] + 1/2*dt, y[k] + 1/2*dt*m1, *args)
        y[k+1] = y[k] + dt*m2
    return [t, y]


# Methode van Runge-Kutta voor één 1ste orde DV
def mijnRK(f, interval, y0, n, *args):
    # Benaderen van y' = f(t, y) over tintval = [t0, tn]
    # met beginwaarde y(t0) = y0 en stapgrootte dt = (tn - t0)/n
    # Syntax: [t, y] = mijnRK(f, tintval, y0, n)
    # Input:  f: fuction handle naar het rechterlid van de DV
    #         tintval = [tmin, tmax]: grenzen van de grafiek
    #         y0 = beginwaarde
    #         n = aantal stappen
    #         *args= extra argumenten voor de functie f
    # Output: t = tijdsvector
    #         y = benadering van y' = f(t,y)
    y = np.zeros(n+1)
    y[0] = y0
    dt = (interval[1] - interval[0])/n
    t = np.linspace(interval[0], interval[1], n+1, endpoint = True)
    for k in range(n):
        m1 = f(t[k], y[k], *args)
        m2 = f(t[k] + dt/2, y[k] + dt/2*m1, *args)
        m3 = f(t[k] + dt/2, y[k] + dt/2*m2, *args)
        m4 = f(t[k] + dt,   y[k] + dt*m3, *args)
        y[k+1] = y[k] + dt*(m1 + 2*m2 + 2*m3 + m4)/6
    return [t, y]


# Methode van Runge-Kutta voor één 1ste orde DV
def mijnRK4(f, interval, y0, n, *args):
    # Benaderen van y' = f(t, y) over tintval = [t0, tn]
    # met beginwaarde y(t0) = y0 en stapgrootte dt = (tn - t0)/n
    # Syntax: [t, y] = mijnRK4(f, tintval, y0, n)
    # Input:  f: fuction handle naar het rechterlid van de DV
    #         tintval = [tmin, tmax]: grenzen van de grafiek
    #         y0 = beginwaarde
    #         n = aantal stappen
    #         *args= extra argumenten voor de functie f
    # Output: t = tijdsvector
    #         y = benadering van y' = f(t,y)
    y = np.zeros(n+1)
    y[0] = y0
    dt = (interval[1] - interval[0])/n
    t = np.linspace(interval[0], interval[1], n+1, endpoint = True)
    for k in range(n):
        m1 = f(t[k], y[k], *args)
        m2 = f(t[k] + dt/2, y[k] + dt/2*m1, *args)
        m3 = f(t[k] + dt/2, y[k] + dt/2*m2, *args)
        m4 = f(t[k] + dt,   y[k] + dt*m3, *args)
        y[k+1] = y[k] + dt*(m1 + 2*m2 + 2*m3 + m4)/6
    return [t, y]


# Methode van Runge-Kutta voor één 1ste orde DV
def mijnAB(f, interval, y0, n, *args):
    # Benaderen van y' = f(t, y) over tintval = [t0, tn]
    # met beginwaarde y(t0) = y0 en stapgrootte dt = (tn - t0)/n
    # Syntax: [t, y] = mijnAB(f, tintval, y0, n)
    # Input:  f: fuction handle naar het rechterlid van de DV
    #         tintval = [tmin, tmax]: grenzen van de grafiek
    #         y0 = beginwaarde
    #         n = aantal stappen
    #         *args= extra argumenten voor de functie f
    # Output: t = tijdsvector
    #         y = benadering van y' = f(t,y)
    y = np.zeros(n+1)
    y[0] = y0
    dt = (interval[1] - interval[0])/n
    t = np.linspace(interval[0], interval[1], n+1, endpoint = True)
    m1 = f(t[0], y[0], *args)
    m2 = f(t[0] + 1/2*dt, y[0] + 1/2*dt*m1, *args)
    y[1] = y[0] + dt*m2
    for k in range(1, n):
        y[k+1] = y[k] + (dt/2)*( 3*f(t[k], y[k], *args) \
                                - f(t[k-1], y[k-1], *args) )
    return [t, y]


# Functieverschil (voor één 1ste orde DV)
def functieverschil(t, y, exact, *args):
    # Deze functie vergelijkt benadering y en de exacte oplossing
    # y(t) in de punten t(k), en geeft de absolute waarde van de
    # verschillen in deze punten terug in een kolomvector.
    # Syntax: v = functieverschil(t, y, exact)
    # Input:  t = tijdsvector
    #         y = benadering van y' = f(t,y)
    #         exact: function handle naar de exacte oplossing
    #                van de DV
    #         *args= extra argumenten voor de functie exact
    # Output: v = |benadering - exact|
    n = len(t)
    fout = np.zeros(n)
    for k in range(n):
        fout[k] = abs(y[k] - exact(t[k], *args))
    return fout


# DIT IS DE NIEUWE fasevlak FUNCTIE !!!        
# Fasevlak voor een stelsel van twee 1ste orde DVn
def fasevlak(f, tintval, xintval, yintval, *args, n=25, ax=None, color='b'):
    # Deze functie plot het richtingsveld van 
    # een stelsel van twee eerste-orde DVn in het fasevlak
    # Syntax: fasevlak(f, xinterval, yinterval, tinterval, *args, 
    #                                           n=25, ax=ax, color=color):
    # Input:  f: function handle naar het rl. van het stelsel DVn
    #         grenzen van de grafiek:
    #         tintval = [tmin, tmax]
    #         xintval = [xmin, xmax]
    #         yintval = [ymin, ymax]
    #         *args = eventueel bijkomende argumenten voor de functie f
    #         n = precisie van het fasevlak
    #         ax = handle naar een axis-object
    #         color = kleur van de pijltjes
    # Output: fasevlak van het stelsel eerste-orde DVn
    tmin = tintval[0]
    tmax = tintval[1]
    xmin = xintval[0]
    xmax = xintval[1]
    ymin = yintval[0]
    ymax = yintval[1]
    T = np.linspace(tmin, tmax, n)
    x = np.linspace(xmin, xmax, n)
    y = np.linspace(ymin, ymax, n)
    [X, Y] = np.meshgrid(x, y)
    U = np.zeros(X.shape)
    V = np.zeros(Y.shape)
    NI, NJ = Y.shape
    if ax:
        for t in T:
            ax.cla()
            for i in range(NI):
                for j in range(NJ):
                    x = X[i, j]
                    y = Y[i, j]
                    yprime = f(t, [x, y], *args)
                    U[i,j] = yprime[0]
                    V[i,j] = yprime[1]
            N = np.sqrt(U**2 + V**2)
            N[N == 0] = 1.0
            U=U/N
            V=V/N
            ax.quiver(X, Y, U, V, angles='xy', color=color, headwidth=5)
            tijd = np.fix(t*100)/100
            ax.set_title('Fasevlak bij t = ' + str(tijd))
            ax.set_xlabel('$y_1$')
            ax.set_ylabel('$y_2$')
            ax.set_xlim(xintval)
            ax.set_ylim(yintval)
            plt.pause(0.1)
            plt.show()
    else:
        for t in T:
            plt.clf()
            for i in range(NI):
                for j in range(NJ):
                    x = X[i, j]
                    y = Y[i, j]
                    yprime = f(t, [x, y], *args)
                    U[i,j] = yprime[0]
                    V[i,j] = yprime[1]
            N = np.sqrt(U**2 + V**2)
            N[N == 0] = 1.0
            plt.quiver(X, Y, U/N, V/N, angles='xy', color=color, headwidth=5)
            tijd = np.fix(t*100)/100
            plt.title('Fasevlak bij t = ' + str(tijd))
            plt.xlabel('$y_1$')
            plt.ylabel('$y_2$')
            plt.xlim(xintval)
            plt.ylim(yintval)
            plt.pause(0.1)
            plt.show()


def fasebaan(f, tintval, xintval, yintval, Y0, *args,\
             n=25, ax=None, color='b'):
    # Deze functie benadert en plot de fasebaan doorheen het
    # punt Y0 = [y10, y20] in het fasevlak.
    # Hiervoor wordt gebruik gemaakt van de RK4-benadering in n
    # stappen over het tijdsinterval tintval.
    # Syntax: fasevlak(f, tintval, xintval, yintval, *args, 
    #                  n=25, ax=ax, color=color):
    # Input:  f: function handle naar het rl. van het stelsel DVn
    #         grenzen van de grafiek:
    #         tintval = [tmin, tmax]
    #         xintval = [xmin, xmax]
    #         yintval = [ymin, ymax]
    #         Y0 = beginvoorwaarden [y10, y20]
    #         *args = eventueel bijkomende argumenten voor f
    #         n = precisie van het fasevlak
    #         ax = handle naar een axis-object
    #         color = kleur van de pijltjes
    # Output: fasebaan doorheen Y0 = [y10, y20] in het fasevlak
    tmin = tintval[0]
    tmax = tintval[1]
    xmin = xintval[0]
    xmax = xintval[1]
    ymin = yintval[0]
    ymax = yintval[1]
    T = np.linspace(tmin, tmax, n)
    x = np.linspace(xmin, xmax, n)
    y = np.linspace(ymin, ymax, n)
    [X, Y] = np.meshgrid(x, y)
    U = np.zeros(X.shape)
    V = np.zeros(Y.shape)
    NI, NJ = Y.shape
    for t in T:
        ax.cla()
        for i in range(NI):
            for j in range(NJ):
                x = X[i, j]
                y = Y[i, j]
                yprime = f(t, [x, y], *args)
                U[i,j] = yprime[0]
                V[i,j] = yprime[1]
        N = np.sqrt(U**2 + V**2)
        N[N == 0] = 1.0
        U=U/N
        V=V/N
        ax.quiver(X, Y, U, V, angles='xy', color=color, headwidth=5)
        ax.plot(Y0[0], Y0[1], 'ko', mfc="None")
        [t_rk, Y_rk] = mijnRK4Stelsel(f, [tmin, t], Y0, n)
        ax.plot(Y_rk[:,0], Y_rk[:,1], 'k-')
        ax.plot(Y_rk[-1,0], Y_rk[-1,1], 'ko')
        plt.show()
        tijd = np.fix(t*100)/100
        ax.set_title("Fasevlak bij t = {:0.2f}".format(tijd))
        ax.set_xlabel('$y_1$')
        ax.set_ylabel('$y_2$')
        ax.set_xlim(xintval)
        ax.set_ylim(yintval)
        plt.pause(0.1)
        plt.show()


# Methode van Euler voor een stelsel 1ste orde DVn
def mijnEulerStelsel(f, tintval, Y0, n, *args):
    # Benaderen van het stelsel y' = f(t, y) over
    # intval = [t0, tn] met beginwaarden y(t0) = Y0
    # stapgrootte dt = (tn - t0)/n
    # Syntax: [t, Y] = mijnEulerStelsel(f, tintval, Y0, n)
    # Input:  f: function handle naar het rl. van het stelsel DVn
    #         tintval = [tmin, tmax]: grenzen van de grafiek
    #         Y0 = beginwaarden
    #         n = aantal stappen
    #         *args= extra argumenten voor de functie f
    # Output: t = tijdsvector
    #         Y = benadering (elke kolom van Y is een benadering
    #             van een DV van het stelsel DVn)
    m = len(Y0)
    Y = np.zeros((n+1, m))
    Y[0,:] = Y0
    dt = (tintval[1] - tintval[0])/n
    t = np.linspace(tintval[0], tintval[1], n+1, endpoint=True)
    for k in range(n):
        m1 = f(t[k], Y[k,:], *args)
        Y[k+1,:] = Y[k,:] + dt*m1
    return [t, Y]


# Methode van Runge-Kutta voor een stelsel 1ste orde DVn
def mijnRKStelsel(f, tintval, Y0, n, *args):
    # Benaderen van het stelsel y' = f(t, y) over
    # intval = [t0, tn] met beginwaarden y(t0) = Y0
    # stapgrootte dt = (tn - t0)/n
    # Syntax: [t, Y] = mijnRKStelsel(f, tintval, Y0, n)
    # Input:  f: function handle naar het rl. van het stelsel DVn
    #         tintval = [tmin, tmax]: grenzen van de grafiek
    #         Y0 = beginwaarden
    #         n = aantal stappen
    #         *args= extra argumenten voor de functie f
    # Output: t = tijdsvector
    #         Y = benadering (elke kolom van Y is een benadering
    #             van een DV van het stelsel DVn)
    m = len(Y0)
    Y = np.zeros((n+1, m))
    Y[0,:] = Y0
    dt = (tintval[1] - tintval[0])/n
    t = np.linspace(tintval[0], tintval[1], n+1, endpoint=True)
    for k in range(n):
        m1 = f(t[k], Y[k,:], *args) 
        m2 = f(t[k]+dt/2, Y[k,:]+dt*m1/2, *args)
        m3 = f(t[k]+dt/2, Y[k,:]+dt*m2/2, *args)
        m4 = f(t[k]+dt, Y[k,:]+dt*m3, *args)
        Y[k+1,:] = Y[k,:] + dt*(m1/6 + m2/3 + m3/3 + m4/6)
    t = np.transpose(t)
    return [t, Y]


# Methode van Runge-Kutta voor een stelsel 1ste orde DVn
def mijnRK4Stelsel(f, tintval, Y0, n, *args):
    # Benaderen van het stelsel y' = f(t, y) over
    # intval = [t0, tn] met beginwaarden y(t0) = Y0
    # stapgrootte dt = (tn - t0)/n
    # Syntax: [t, Y] = mijnRK4Stelsel(f, tintval, Y0, n)
    # Input:  f: function handle naar het rl. van het stelsel DVn
    #         tintval = [tmin, tmax]: grenzen van de grafiek
    #         Y0 = beginwaarden
    #         n = aantal stappen
    #         *args= extra argumenten voor de functie f
    # Output: t = tijdsvector
    #         Y = benadering (elke kolom van Y is een benadering
    #             van een DV van het stelsel DVn)
    m = len(Y0)
    Y = np.zeros((n+1, m))
    Y[0,:] = Y0
    dt = (tintval[1] - tintval[0])/n
    t = np.linspace(tintval[0], tintval[1], n+1, endpoint=True)
    for k in range(n):
        m1 = f(t[k], Y[k,:], *args) 
        m2 = f(t[k]+dt/2, Y[k,:]+dt*m1/2, *args)
        m3 = f(t[k]+dt/2, Y[k,:]+dt*m2/2, *args)
        m4 = f(t[k]+dt, Y[k,:]+dt*m3, *args)
        Y[k+1,:] = Y[k,:] + dt*(m1/6 + m2/3 + m3/3 + m4/6)
    t = np.transpose(t)
    return [t, Y]

# -*- coding: utf-8 -*-
"""
EXAMEN DIFFERENTIAALVERGELIJKINGEN - 1STE ZITTIJD - DEEL 2 - VRAAG 7

Naam student: ...............................

>>> PLAATS UIT COMMENTAAR WAAR NODIG EN VUL VERDER AAN <<<

"""

#%%############################################################################
# Imports
###############################################################################

import numpy as np
import sympy as sp
import matplotlib.pyplot as plt
from functies import fasevlak
from scipy.integrate import solve_ivp

import matplotlib
matplotlib.use("qtagg")
plt.ion()

#%%############################################################################
# Functiedefinities
###############################################################################

# functie voor het eerste deel
def rabbitsfoxes(t, Y):
    R = ...
    F = ...
    Rdot = ...
    Fdot = ...
    return np.array([..., ...])

# functies voor het tweede deel
def step(t, t1, v0, v1):
    return v0 + np.heaviside(t-t1, 0)*(v1-v0)

def rabbitsfoxesadapt(t, Y):
    R = ...
    F = ...
    A = ...
    s = step(t, t1, s0, s0*(1+p))
    Rdot = ...
    Fdot = ...
    Adot = ...
    return np.array([..., ..., ...])

#%%############################################################################
# Instructies
###############################################################################

# Parameterwaarden konijnen
r = 0.54
Kr = 400
d = 0.04
a = 0.011
# Parameterwaarden vossen
f = 0.1
Kf = 80
m = 0.31
y = 0.254
# Parameterwaarden adaptatie
s0 = 0.0008
l = 0.14
p = 0.25

# Beginvoorwaarden
R0 = 370           # Aantal konijnen bij t=0
F0 = 6             # Aantal vossen bij t=0
A0 = 0.0           # Adaptatieniveau bij t=0

# Tijdstip waarop de konijnen de drugs vinden
t1 = 40

# %% (a) Implementatie RL van stelsel DVn

# Zie hierboven

# %% (b) Evenwichtspunten - alle oplossingen

Re, Fe = sp.symbols("R_e, F_e")
eq1 = sp.Eq(..., ...)
eq2 = sp.Eq(..., ...)
evpn = sp.solve([..., ...], [..., ...])
evpn

# %% (b) Evenwichtspunten - strikt positieve oplossingen
Re_val = float(evpn[...][...])
Fe_val = float(evpn[...][...])

print(f"Evenwichtswaarde Re = {Re_val:.0f} konijnen")
print(f"Evenwichtswaarde Fe = {Fe_val:.0f} vossen")

# %% (c) Fasevlak

fig1, ax1 = plt.subplots(num=1)
fasevlak(..., [96, 96], ..., ..., n=28, ax=ax1)
ax1.set_xlabel("R")
ax1.set_ylabel("F")
fig1.savefig('figuur_fasevlak_R_F.png')

# Stabiliteit: ...

# %% (d) Oplossingsbenadering

t_span = ...
n = ...
atol = 1e-6; rtol = 1e-9
t_eval = ...

y0 = [..., ...]
sol = solve_ivp(..., ..., ..., t_eval=..., atol=..., rtol=...)
t_arr = ...      # tijd array
R_arr = ...      # benadering voor R
F_arr = ...      # benadering voor F

fig2, ax2 = plt.subplots(num=2, figsize=(6, 6))
ax2.plot(..., ..., 'k-', label=r"$R$")
ax2.plot(..., ..., 'r-', label=r"$F$")
ax2.set_ylim([0, 400])
ax2.set_xlabel("Tijd (maanden)")
ax2.set_ylabel("Populatiegrootte")
ax2.legend()
ax2.grid(True)
fig2.savefig('figuur_oplossing_R_F.png')

# %% (e) Implementatie RL van stelsel DVn

# Zie hierboven

# %% (f) Oplossingsbenadering

y0 = [..., ..., ...]
sol = solve_ivp(..., ..., ..., t_eval=..., atol=atol, rtol=rtol)
t_arr = ...   # tijd array
R_arr = ...   # benadering voor R
F_arr = ...   # benadering voor F
A_arr = ...   # benadering voor A

fig3, (ax31, ax32) = plt.subplots(num=3, ncols=1, nrows=2, figsize=(6, 6))
ax31.plot(..., ..., 'k-', label=r"$R$")
ax31.plot(..., ..., 'r-', label=r"$F$")
ax31.set_xlabel("Tijd (maanden)")
ax31.set_ylabel("Populatiegrootte")
ax31.legend()
ax31.grid(True)
ax32.plot(..., ..., 'b-', label=r"$A$")
ax32.set_xlabel("Tijd (maanden)")
ax32.set_ylabel("Aanpasbaarheid")
ax32.legend()
ax32.grid(True)
fig3.tight_layout()
fig3.savefig('figuur_oplossing_R_F_A_1.png')

# %% (g) Zoek iteratief naar de minimale p-waarde opdat het finaal
#        aantal konijnen (op t = 96) minstens een waarde van 160 heeft

dp = 0.001    # stapgrootte

sol = solve_ivp(..., ..., ..., t_eval=..., rtol=rtol, atol=atol)
R_arr = ...   # benadering voor R
while ...:
    p = ...
    sol = solve_ivp(..., ..., ..., t_eval=..., rtol=rtol, atol=atol)
    R_arr = ...   # benadering voor R

F_arr = ...   # benadering voor F
A_arr = ...   # benadering voor A

p_min = p     # gevonden minimum waarde van p
p = 0.25      # reset p naar oorspronkelijke waarde 

print(f"p = {p_min:.3f}")     # print de gevonden minimum waarde van p

fig4, (ax41, ax42) = plt.subplots(num=4, ncols=1, nrows=2, figsize=(6, 6))
ax41.plot(..., ..., 'k-', label=r"$R$")
ax41.plot(..., ..., 'r-', label=r"$F$")
ax41.set_xlabel("Tijd (maanden)")
ax41.set_ylabel("Populatiegrootte")
ax41.legend()
ax41.grid(True)
ax42.plot(..., ..., 'b-', label=r"$A$")
ax42.set_xlabel("Tijd (maanden)")
ax42.set_ylabel("Aanpasbaarheid")
ax42.legend()
ax42.grid(True)
fig4.tight_layout()
fig4.savefig('figuur_oplossing_R_F_A_2.png')

PC deel examen

functies.py is bijlage bestand om de andere correct te kunnen uitvoeren. Hierin moet geen code worden geschreven

Bijlagen

📦 Download originele bestanden (ZIP)

import numpy as np
import matplotlib.pyplot as plt


# Richtingsveld (voor één 1ste orde DV)
def richtingsveld(f, tintval, yintval, *args, n=25, ax=None, color='k'):
    # Deze functie plot het richtingsveld van y' = f(t,y)
    # Syntax: richtingsveld(f, tintval, yintval, n)
    # Input: f: fuction handle naar het rechterlid van de DV
    #        grenzen van de grafiek:
    #        tintval = [tmin, tmax]
    #        yintval = [ymin, ymax]
    #        *args= extra argumenten voor de functie f
    #        n = precisie van het richtingsveld, default n = 25
    #        ax = axis-object indien verschillende plots in een
    #             figuur, default ax=None
    #        color = kleurcode pijltjes, default color='k' (zwart)
    # Output: richtingsveld van y' = f(t,y)
    tmin = tintval[0]
    tmax = tintval[1]
    ymin = yintval[0]
    ymax = yintval[1]
    T = np.linspace(tmin, tmax, n)
    y = np.linspace(ymin, ymax, n)
    # Vector field
    [X, Y] = np.meshgrid(T, y)
    U = 1
    V = f(T, Y, *args)
    if V.ndim == 1:
        m = V.shape[0]
        V = np.reshape(V, (1, m))
        V = np.repeat(V, repeats=m, axis=0)
    # Normalize arrows
    N = np.sqrt(U**2 + V**2)  
    U = U/N
    V = V/N
    if ax:
        ax.quiver(X, Y, U, V, angles='xy', color=color)    
        ax.set_xlim(tintval)
        ax.set_ylim(yintval)
        ax.set_xlabel('$t$')
        ax.set_ylabel('$y$')
    else:
        plt.quiver(X, Y, U, V, angles='xy')    
        plt.xlim(tintval)
        plt.ylim(yintval)
        plt.xlabel('$t$')
        plt.ylabel('$y$')
        plt.show()


# Methode van Euler voor één 1ste orde DV
def mijnEuler(f, interval, y0, n, *args):
    # Benaderen van y' = f(t, y) over tintval = [t0, tn]
    # met beginwaarde y(t0) = y0 en stapgrootte dt = (tn - t0)/n
    # Syntax: [t, y] = mijnEuler(f, tintval, y0, n)
    # Input:  f: fuction handle naar het rechterlid van de DV
    #         tintval = [tmin, tmax]: grenzen van de grafiek
    #         y0 = beginwaarde
    #         n = aantal stappen
    #         *args= extra argumenten voor de functie f
    # Output: t = tijdsvector
    #         y = benadering van y' = f(t,y)
    y = np.zeros(n+1)      # pre-allocatie van y als array
    y[0] = y0                 # y0 is het eerste element van
    dt = (interval[1] - interval[0])/n
    t = np.linspace(interval[0], interval[1], n+1, endpoint = True)
    for k in range(0, n):     # k = 0, 1, 2, ..., n-1
        y[k+1] = y[k] + dt*f(t[k], y[k], *args)
    return [t, y]


# Methode van Euler voor één 1ste orde DV
# def mijnEuler(f, interval, y0, n, *args):
#     # Benaderen van y' = f(t, y) over tintval = [t0, tn]
#     # met beginwaarde y(t0) = y0 en stapgrootte dt = (tn - t0)/n
#     # Syntax: [t, y] = mijnEuler(f, tintval, y0, n)
#     # Input:  f: fuction handle naar het rechterlid van de DV
#     #         tintval = [tmin, tmax]: grenzen van de grafiek
#     #         y0 = beginwaarde
#     #         n = aantal stappen
#     #         *args= extra argumenten voor de functie f
#     # Output: t = tijdsvector
#     #         y = benadering van y' = f(t,y)
#     y = np.zeros(n+1)     # pre-allocatie van y als array
#     y[0] = y0                # y0 is het eerste element van
#                              # de oplossingsarray y
#     dt = (interval[1] - interval[0])/n
#     t = ...                  # maak hier de array t aan
#     for k in ...:            # k = 0, 1, 2, ..., n-1
#         # schrijf hier een for-lus die de waarden y[k+1] berekent
#         # op basis van t[k] en y[k]
#         ...
#     return [t, y]


# Midpoint methode voor één 1ste orde DV
def mijnMidpoint(f, interval, y0, n, *args):
    # Benaderen van y' = f(t, y) over tintval = [t0, tn]
    # met beginwaarde y(t0) = y0 en stapgrootte dt = (tn - t0)/n
    # Syntax: [t, y] = mijnMidpoint(f, tintval, y0, n)
    # Input:  f: fuction handle naar het rechterlid van de DV
    #         tintval = [tmin, tmax]: grenzen van de grafiek
    #         y0 = beginwaarde
    #         n = aantal stappen
    #         *args= extra argumenten voor de functie f
    # Output: t = tijdsvector
    #         y = benadering van y' = f(t,y)
    y = np.zeros(n+1)
    y[0] = y0
    dt = (interval[1] - interval[0])/n
    t = np.linspace(interval[0], interval[1], n+1, endpoint = True)
    for k in range(0, n): # k = 0, 1, 2, ..., n-1
        # >>> alpha = 1 <<<
        # m1 = f(t[k], y[k])
        # m2 = f(t[k] + dt, y[k] + dt*m1)
        # y[k+1] = y[k] + dt/2*(m1 + m2)
        # >>> alpha = 1/2 <<<
        m1 = f(t[k], y[k], *args)
        m2 = f(t[k] + 1/2*dt, y[k] + 1/2*dt*m1, *args)
        y[k+1] = y[k] + dt*m2
    return [t, y]


# Methode van Runge-Kutta voor één 1ste orde DV
def mijnRK(f, interval, y0, n, *args):
    # Benaderen van y' = f(t, y) over tintval = [t0, tn]
    # met beginwaarde y(t0) = y0 en stapgrootte dt = (tn - t0)/n
    # Syntax: [t, y] = mijnRK(f, tintval, y0, n)
    # Input:  f: fuction handle naar het rechterlid van de DV
    #         tintval = [tmin, tmax]: grenzen van de grafiek
    #         y0 = beginwaarde
    #         n = aantal stappen
    #         *args= extra argumenten voor de functie f
    # Output: t = tijdsvector
    #         y = benadering van y' = f(t,y)
    y = np.zeros(n+1)
    y[0] = y0
    dt = (interval[1] - interval[0])/n
    t = np.linspace(interval[0], interval[1], n+1, endpoint = True)
    for k in range(n):
        m1 = f(t[k], y[k], *args)
        m2 = f(t[k] + dt/2, y[k] + dt/2*m1, *args)
        m3 = f(t[k] + dt/2, y[k] + dt/2*m2, *args)
        m4 = f(t[k] + dt,   y[k] + dt*m3, *args)
        y[k+1] = y[k] + dt*(m1 + 2*m2 + 2*m3 + m4)/6
    return [t, y]


# Methode van Runge-Kutta voor één 1ste orde DV
def mijnRK4(f, interval, y0, n, *args):
    # Benaderen van y' = f(t, y) over tintval = [t0, tn]
    # met beginwaarde y(t0) = y0 en stapgrootte dt = (tn - t0)/n
    # Syntax: [t, y] = mijnRK4(f, tintval, y0, n)
    # Input:  f: fuction handle naar het rechterlid van de DV
    #         tintval = [tmin, tmax]: grenzen van de grafiek
    #         y0 = beginwaarde
    #         n = aantal stappen
    #         *args= extra argumenten voor de functie f
    # Output: t = tijdsvector
    #         y = benadering van y' = f(t,y)
    y = np.zeros(n+1)
    y[0] = y0
    dt = (interval[1] - interval[0])/n
    t = np.linspace(interval[0], interval[1], n+1, endpoint = True)
    for k in range(n):
        m1 = f(t[k], y[k], *args)
        m2 = f(t[k] + dt/2, y[k] + dt/2*m1, *args)
        m3 = f(t[k] + dt/2, y[k] + dt/2*m2, *args)
        m4 = f(t[k] + dt,   y[k] + dt*m3, *args)
        y[k+1] = y[k] + dt*(m1 + 2*m2 + 2*m3 + m4)/6
    return [t, y]


# Methode van Runge-Kutta voor één 1ste orde DV
def mijnAB(f, interval, y0, n, *args):
    # Benaderen van y' = f(t, y) over tintval = [t0, tn]
    # met beginwaarde y(t0) = y0 en stapgrootte dt = (tn - t0)/n
    # Syntax: [t, y] = mijnAB(f, tintval, y0, n)
    # Input:  f: fuction handle naar het rechterlid van de DV
    #         tintval = [tmin, tmax]: grenzen van de grafiek
    #         y0 = beginwaarde
    #         n = aantal stappen
    #         *args= extra argumenten voor de functie f
    # Output: t = tijdsvector
    #         y = benadering van y' = f(t,y)
    y = np.zeros(n+1)
    y[0] = y0
    dt = (interval[1] - interval[0])/n
    t = np.linspace(interval[0], interval[1], n+1, endpoint = True)
    m1 = f(t[0], y[0], *args)
    m2 = f(t[0] + 1/2*dt, y[0] + 1/2*dt*m1, *args)
    y[1] = y[0] + dt*m2
    for k in range(1, n):
        y[k+1] = y[k] + (dt/2)*( 3*f(t[k], y[k], *args) - f(t[k-1], y[k-1], *args) )
    return [t, y]


# Functieverschil (voor één 1ste orde DV)
def functieverschil(t, y, exact, *args):
    # Deze functie vergelijkt benadering y en de exacte oplossing
    # y(t) in de punten t(k), en geeft de absolute waarde van de
    # verschillen in deze punten terug in een kolomvector.
    # Syntax: v = functieverschil(t, y, exact)
    # Input:  t = tijdsvector
    #         y = benadering van y' = f(t,y)
    #         exact: function handle naar de exacte oplossing
    #                van de DV
    #         *args= extra argumenten voor de functie exact
    # Output: v = |benadering - exact|
    n = len(t)
    fout = np.zeros(n)
    for k in range(n):
        fout[k] = abs(y[k] - exact(t[k], *args))
    return fout


# DIT IS DE NIEUWE fasevlak FUNCTIE !!!        
# Fasevlak voor een stelsel van twee 1ste orde DVn
def fasevlak(f, tintval, xintval, yintval, *args, n=25, ax=None, color='b'):
    # Deze functie plot het richtingsveld van 
    # een stelsel van twee eerste-orde DVn in het fasevlak
    # Syntax: fasevlak(f, xinterval, yinterval, tinterval, *args, 
    #                                           n=25, ax=ax, color=color):
    # Input:  f: function handle naar het rl. van het stelsel DVn
    #         grenzen van de grafiek:
    #         tintval = [tmin, tmax]
    #         xintval = [xmin, xmax]
    #         yintval = [ymin, ymax]
    #         *args = eventueel bijkomende argumenten voor de functie f
    #         n = precisie van het fasevlak
    #         ax = handle naar een axis-object
    #         color = kleur van de pijltjes
    # Output: fasevlak van het stelsel eerste-orde DVn
    tmin = tintval[0]
    tmax = tintval[1]
    xmin = xintval[0]
    xmax = xintval[1]
    ymin = yintval[0]
    ymax = yintval[1]
    T = np.linspace(tmin, tmax, n)
    x = np.linspace(xmin, xmax, n)
    y = np.linspace(ymin, ymax, n)
    [X, Y] = np.meshgrid(x, y)
    U = np.zeros(X.shape)
    V = np.zeros(Y.shape)
    NI, NJ = Y.shape
    if ax:
        for t in T:
            ax.cla()
            for i in range(NI):
                for j in range(NJ):
                    x = X[i, j]
                    y = Y[i, j]
                    yprime = f(t, [x, y], *args)
                    U[i,j] = yprime[0]
                    V[i,j] = yprime[1]
            N = np.sqrt(U**2 + V**2)
            N[N == 0] = 1.0
            U=U/N
            V=V/N
            ax.quiver(X, Y, U, V, angles='xy', color=color, headwidth=5)
            tijd = np.fix(t*100)/100
            ax.set_title('Fasevlak bij t = ' + str(tijd))
            ax.set_xlabel('$y_1$')
            ax.set_ylabel('$y_2$')
            ax.set_xlim(xintval)
            ax.set_ylim(yintval)
            plt.pause(0.1)
            plt.show()
    else:
        for t in T:
            plt.clf()
            for i in range(NI):
                for j in range(NJ):
                    x = X[i, j]
                    y = Y[i, j]
                    yprime = f(t, [x, y], *args)
                    U[i,j] = yprime[0]
                    V[i,j] = yprime[1]
            N = np.sqrt(U**2 + V**2)
            N[N == 0] = 1.0
            plt.quiver(X, Y, U/N, V/N, angles='xy', color=color, headwidth=5)
            tijd = np.fix(t*100)/100
            plt.title('Fasevlak bij t = ' + str(tijd))
            plt.xlabel('$y_1$')
            plt.ylabel('$y_2$')
            plt.xlim(xintval)
            plt.ylim(yintval)
            plt.pause(0.1)
            plt.show()


def fasebaan(f, tintval, xintval, yintval, Y0, *args,\
             n=25, ax=None, color='b'):
    # Deze functie benadert en plot de fasebaan doorheen het
    # punt Y0 = [y10, y20] in het fasevlak.
    # Hiervoor wordt gebruik gemaakt van de RK4-benadering in n
    # stappen over het tijdsinterval tintval.
    # Syntax: fasevlak(f, tintval, xintval, yintval, *args, 
    #                  n=25, ax=ax, color=color):
    # Input:  f: function handle naar het rl. van het stelsel DVn
    #         grenzen van de grafiek:
    #         tintval = [tmin, tmax]
    #         xintval = [xmin, xmax]
    #         yintval = [ymin, ymax]
    #         Y0 = beginvoorwaarden [y10, y20]
    #         *args = eventueel bijkomende argumenten voor f
    #         n = precisie van het fasevlak
    #         ax = handle naar een axis-object
    #         color = kleur van de pijltjes
    # Output: fasebaan doorheen Y0 = [y10, y20] in het fasevlak
    tmin = tintval[0]
    tmax = tintval[1]
    xmin = xintval[0]
    xmax = xintval[1]
    ymin = yintval[0]
    ymax = yintval[1]
    T = np.linspace(tmin, tmax, n)
    x = np.linspace(xmin, xmax, n)
    y = np.linspace(ymin, ymax, n)
    [X, Y] = np.meshgrid(x, y)
    U = np.zeros(X.shape)
    V = np.zeros(Y.shape)
    NI, NJ = Y.shape
    for t in T:
        ax.cla()
        for i in range(NI):
            for j in range(NJ):
                x = X[i, j]
                y = Y[i, j]
                yprime = f(t, [x, y], *args)
                U[i,j] = yprime[0]
                V[i,j] = yprime[1]
        N = np.sqrt(U**2 + V**2)
        N[N == 0] = 1.0
        U=U/N
        V=V/N
        ax.quiver(X, Y, U, V, angles='xy', color=color, headwidth=5)
        ax.plot(Y0[0], Y0[1], 'ko', mfc="None")
        [t_rk, Y_rk] = mijnRK4Stelsel(f, [tmin, t], Y0, n)
        ax.plot(Y_rk[:,0], Y_rk[:,1], 'k-')
        ax.plot(Y_rk[-1,0], Y_rk[-1,1], 'ko')
        plt.show()
        tijd = np.fix(t*100)/100
        ax.set_title("Fasevlak bij t = {:0.2f}".format(tijd))
        ax.set_xlabel('$y_1$')
        ax.set_ylabel('$y_2$')
        ax.set_xlim(xintval)
        ax.set_ylim(yintval)
        plt.pause(0.1)
        plt.show()


# Methode van Euler voor een stelsel 1ste orde DVn
def mijnEulerStelsel(f, tintval, Y0, n, *args):
    # Benaderen van het stelsel y' = f(t, y) over
    # intval = [t0, tn] met beginwaarden y(t0) = Y0
    # stapgrootte dt = (tn - t0)/n
    # Syntax: [t, Y] = mijnEulerStelsel(f, tintval, Y0, n)
    # Input:  f: function handle naar het rl. van het stelsel DVn
    #         tintval = [tmin, tmax]: grenzen van de grafiek
    #         Y0 = beginwaarden
    #         n = aantal stappen
    #         *args= extra argumenten voor de functie f
    # Output: t = tijdsvector
    #         Y = benadering (elke kolom van Y is een benadering
    #             van een DV van het stelsel DVn)
    m = len(Y0)
    Y = np.zeros((n+1, m))
    Y[0,:] = Y0
    dt = (tintval[1] - tintval[0])/n
    t = np.linspace(tintval[0], tintval[1], n+1, endpoint=True)
    for k in range(n):
        m1 = f(t[k], Y[k,:], *args)
        Y[k+1,:] = Y[k,:] + dt*m1
    return [t, Y]


# Methode van Runge-Kutta voor een stelsel 1ste orde DVn
def mijnRKStelsel(f, tintval, Y0, n, *args):
    # Benaderen van het stelsel y' = f(t, y) over
    # intval = [t0, tn] met beginwaarden y(t0) = Y0
    # stapgrootte dt = (tn - t0)/n
    # Syntax: [t, Y] = mijnRKStelsel(f, tintval, Y0, n)
    # Input:  f: function handle naar het rl. van het stelsel DVn
    #         tintval = [tmin, tmax]: grenzen van de grafiek
    #         Y0 = beginwaarden
    #         n = aantal stappen
    #         *args= extra argumenten voor de functie f
    # Output: t = tijdsvector
    #         Y = benadering (elke kolom van Y is een benadering
    #             van een DV van het stelsel DVn)
    m = len(Y0)
    Y = np.zeros((n+1, m))
    Y[0,:] = Y0
    dt = (tintval[1] - tintval[0])/n
    t = np.linspace(tintval[0], tintval[1], n+1, endpoint=True)
    for k in range(n):
        m1 = f(t[k], Y[k,:], *args) 
        m2 = f(t[k]+dt/2, Y[k,:]+dt*m1/2, *args)
        m3 = f(t[k]+dt/2, Y[k,:]+dt*m2/2, *args)
        m4 = f(t[k]+dt, Y[k,:]+dt*m3, *args)
        Y[k+1,:] = Y[k,:] + dt*(m1/6 + m2/3 + m3/3 + m4/6)
    t = np.transpose(t)
    return [t, Y]


# Methode van Runge-Kutta voor een stelsel 1ste orde DVn
def mijnRK4Stelsel(f, tintval, Y0, n, *args):
    # Benaderen van het stelsel y' = f(t, y) over
    # intval = [t0, tn] met beginwaarden y(t0) = Y0
    # stapgrootte dt = (tn - t0)/n
    # Syntax: [t, Y] = mijnRK4Stelsel(f, tintval, Y0, n)
    # Input:  f: function handle naar het rl. van het stelsel DVn
    #         tintval = [tmin, tmax]: grenzen van de grafiek
    #         Y0 = beginwaarden
    #         n = aantal stappen
    #         *args= extra argumenten voor de functie f
    # Output: t = tijdsvector
    #         Y = benadering (elke kolom van Y is een benadering
    #             van een DV van het stelsel DVn)
    m = len(Y0)
    Y = np.zeros((n+1, m))
    Y[0,:] = Y0
    dt = (tintval[1] - tintval[0])/n
    t = np.linspace(tintval[0], tintval[1], n+1, endpoint=True)
    for k in range(n):
        m1 = f(t[k], Y[k,:], *args) 
        m2 = f(t[k]+dt/2, Y[k,:]+dt*m1/2, *args)
        m3 = f(t[k]+dt/2, Y[k,:]+dt*m2/2, *args)
        m4 = f(t[k]+dt, Y[k,:]+dt*m3, *args)
        Y[k+1,:] = Y[k,:] + dt*(m1/6 + m2/3 + m3/3 + m4/6)
    t = np.transpose(t)
    return [t, Y]

# -*- coding: utf-8 -*-
"""
EXAMEN DIFFERENTIAALVERGELIJKINGEN - 1STE ZITTIJD - DEEL 2 - VRAAG 7

Naam student: ...............................

>>> PLAATS UIT COMMENTAAR WAAR NODIG EN VUL VERDER AAN <<<

"""

#%%############################################################################
# Imports
###############################################################################

import numpy as np
import sympy as sp
import matplotlib.pyplot as plt
from functies import fasevlak
from scipy.integrate import solve_ivp

import matplotlib
matplotlib.use("qtagg")
plt.ion()

#%%############################################################################
# Functiedefinities
###############################################################################

def step(t, t1, v0, v1):          # DEZE FUNCTIE NIET WIJZIGEN!
    return v0 + np.heaviside(t-t1, 0)*(v1-v0)

def trees_parasites(t, Y):
    r = 0.16; K = 860; d = 2.0e-4
    b = 2.7e-1; c = 1.2e-4

    T = ...
    P = ...
    e = step(t, t1, e0, e1)       # DEZE INSTRUCTIE NIET WIJZIGEN!
    Tdot = ...
    Pdot = ...
    return np.array([..., ...])

#%%############################################################################
# Instructies
###############################################################################

# Parameterwaarden
r = 0.16; K = 860; d = 2.0e-4
b = 2.7e-1; c = 1.2e-4

t1 = 60.0         # tijdstip t=60
e0 = 0.5          # waarde van e voor 0 <= t <= 60
e1 = 1.0          # waarde van e voor 60 < t <= 120

delta_e = 0.01    # stapgrootte bij deelvraag (e)

# Beginvoorwaarden
T0 = 100          # initieel dichtheid bomen
P0 = 5            # initiaal dichtheid parasieten

#%% (a) Implementatie RL van stelsel DVn

# Zie hierboven

#%% (b) Evenwichtspunten wanneer e = 0.5
# Gebruik de variable e0 voor e want e0 is gelijk aan 0.5

t1 = 60.0; e0 = 0.5; e1 = 1.0

Te, Pe = sp.symbols("Te, Pe")
eq1 = sp.Eq(..., ...)
eq2 = sp.Eq(..., ...)
evpn = sp.solve([..., ...], [..., ...])
evpn

# Kies deze waarvoor Te en Pe strikt positief zijn:
Te_val = round(float(evpn[...][...]))
Pe_val = round(float(evpn[...][...]))
print("Evenwichtswaarde voor de bomen:      {:d}".format(Te_val))
print("Evenwichtswaarde voor de parasieten: {:d}".format(Pe_val))

#%% (c) Fasevlak op tijdstip 60 in het (T x P)-gebied [0, 800] x [0, 800]

t1 = 60.0; e0 = 0.5; e1 = 1.0

fig, ax = plt.subplots(num=1)
fasevlak(..., [60, 60], ..., ..., n=25, ax=ax)   # [60, 60] betekent tijdstip t=60
ax.set_xlabel(r"$T$")
ax.set_ylabel(r"$P$")
fig.savefig('figuur_fasevlak_T_P.png')

# Stabiliteit: ...

#%% (d) Oplossingsbenadering

t1 = 60.0; e0 = 0.5; e1 = 1.0

t_span = ...
n = 300
y0 = [..., ...]
atol = 1e-6; rtol = 1e-9
t_eval = ...   # array van n+1 tijdstippen in [0, 120]
sol = solve_ivp(..., ..., ..., t_eval=..., rtol=rtol, atol=atol)
t_arr = ...    # tijd array, zelfde als t_eval
T_arr = ...    # benadering voor T
P_arr = ...    # benadering voor P

print("Dichtheid bomen bij t = 120 jaar:      {:5d}".format(round(T_arr[-1])))
print("Dichtheid parasieten bij t = 120 jaar: {:5d}".format(round(P_arr[-1])))

# Maak een figuur van T en P i.f.v. t
fig = plt.figure(num=2)
plt.plot(..., ..., 'b-', label=r"$T(t)$")
plt.plot(..., ..., 'g-', label=r"$P(t)$")
plt.xlabel(r"$t$ [$jaar$]")
plt.ylabel("Dichtheid")
plt.ylim([0, 800])
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.savefig('figuur_evolutie_T_P_1.png')

#%% (e) Minimale waarde van de coefficient e opdat de dichtheid aan parasieten
#       onder de 100 blijft.
# Gebruik hier de variable e1 want e1 is de waarde voor e voor 60 < t <= 120.

t1 = 60.0; e0 = 0.5; e1 = ...

sol = solve_ivp(..., ..., ..., t_eval=..., rtol=rtol, atol=atol)
T_arr = ...    # benadering voor T
P_arr = ...    # benadering voor P

while ...:
    e1 = ...      # wijzig e1 met een stapgrootte delta_e
    sol = solve_ivp(..., ..., ..., t_eval=..., rtol=rtol, atol=atol)
    T_arr = ...    # benadering voor T
    P_arr = ...    # benadering voor P

print("Coefficient e: {:3.2f}".format(e1))
print("Dichtheid bomen bij t = 120 jaar:      {:5d}".format(round(T_arr[-1])))
print("Dichtheid parasieten bij t = 120 jaar: {:5d}".format(round(P_arr[-1])))

# Maak een nieuwe figuur van T en P i.f.v. t, gebruikmakende van de nieuwe
# waarde voor de variabele e1
sol = solve_ivp(..., ..., ..., t_eval=..., rtol=rtol, atol=atol)
t_arr = ...    # tijd array, zelfde als t_eval
T_arr = ...    # benadering voor T
P_arr = ...    # benadering voor P

fig = plt.figure(num=3)
plt.plot(..., ..., 'b-', label=r"$T(t)$")
plt.plot(..., ..., 'g-', label=r"$P(t)$")
plt.xlabel(r"$t$ [$jaar$]")
plt.ylabel("Dichtheid")
plt.ylim([0, 800])
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.savefig('figuur_evolutie_T_P_2.png')

# -*- coding: utf-8 -*-
"""
EXAMEN DIFFERENTIAALVERGELIJKINGEN - 1STE ZITTIJD - DEEL 2 - VRAAG 8

Naam student: ...............................

>>> PLAATS UIT COMMENTAAR WAAR NODIG EN VUL VERDER AAN <<<

"""

#%%############################################################################
# Imports
###############################################################################

import numpy as np
import matplotlib.pyplot as plt
from functies import functieverschil
from matplotlib import cm

import matplotlib
matplotlib.use("qtagg")
plt.ion()


#%%############################################################################
# Instructies
###############################################################################

# fysische constanten
D = 1
r = 0.1

# intervalgrootte
dx = 1
dt = 0.1

# initiële populatiedichtheid
U_init = 0.1

# fysische grenzen
xmin = 0
xmax = 100
tmin = 0
tmax = 100

# arrays maken voor de ruimtedimensie, tijdsdimensie, en de populatiedichtheid
x = np.arange(..., ..., ...) # +2 knopen voor elke kant
t = np.arange(..., ..., ...)

# aantal knopen
nx = int(...)
nt = int(...)

U = np.zeros((nx, nt))

# Dirichlet rand- en/of beginvoorwaarde(n)
U[..., ...] = ...


# %% Eindige differentiemethode (kern van het script)

# Loop over alle knopen (denk aan eventuele Neumann-voorwaarde(n)!)
for j in range(..., ...):
    for i in range(..., ...):
        U[..., ...] = ...
    # Gediscretiseerde Neumann-randvoorwaarde(n) (indien van toepassing)
    U[..., ...] = ...
    U[..., ...] = ...


# %% Plot and save resultaten

# enkel niet-imaginaire knopen bijhouden
U = U[..., ...]
x = x[...]

# oppervlakteplot
U_T = np.transpose(U)
x_mg, t_mg = np.meshgrid(..., ...)

fig, ax = plt.subplots(num=1, subplot_kw={"projection": "3d"})
surf = ax.plot_surface(..., ..., ..., cmap=cm.coolwarm,
                       linewidth=0, antialiased=False)
ax.set_xlabel("x")
ax.set_ylabel("t")
ax.set_zlabel("U")
plt.savefig("fisher-kpp-beek.png")


# %% Populatieverhouding wanneer t = 100

# bereken maximale waarde wanneer t = 100
U_max = ...

# bereken en print de verhouding
teller = ...
noemer = ...
verhouding = round(..., 5)

print(f"Populatieverhouding: {verhouding}.")