Model Validation
Read this material on mathematical modeling. Try to complete the exercise at the end of the chapters. Don't worry if you have not yet learned the math necessary to complete the exercises. Rather, notice how many physical systems can be represented by mathematical models. We can then implement these models using programming languages like Python and R. We will discuss those later. Think about the types of systems you interact with, personally or in business situations. How might you use a mathematical model to help you understand this system?
2. Bounded growth (Phase line, Stable equilibrium)
2.3 Reaction time (Euler's Method)
Problem
Back to the rainbowfish for another modelling cycle. Now we want to know how long it takes to breed them. Do you have to invest years to grow a good-sized population, or can it be done in a few weeks?
In section 1 the number of rainbowfish in the aquarium was modelled with a limited growth rate. At last, the solution is plausible, also after an infinitely long time.
In this section we focus on answering the following question:
How long does it take for the rainbowfish population to get into its equilibrium state?
Do you have to invest years to grow a good-sized population, or is the business up and running in a few weeks?
Mathematical model
To recapitulate, the new problem is:
How long does it take for the rainbowfish population to get into its equilibrium state?
Earlier modelling cycles have resulted in the following information:
- P(t) is the number of rainbowfish in the aquarium, with t in days.
- We start with 30 fish. The unbounded birth rate of the fish is b=0.7 per day.
- The death rate is ignored: d≈0.
- The tank can sustain a population of 750 rainbowfish.
- Each day 20 fish are sold.
The resulting initial value problem is:
\(\dfrac{dP}{dt} = 0.7 P (1 - \dfrac{P}{750} - 20\), \(P(0) = 30\).
The differential equation has two equilibrium solutions, namely an unstable one of approximately 29.75 and a stable one of approximately 720.25. The second one is relevant for the new problem.
How long does it take to reach that equilibrium? Mathematically, the answer is: "infinitely long" ;-). The nearer P gets to the equilibrium, the lower the value of the derivative (but it stays positive), so the growth will keep slowing down, never reaching zero growth, and P will never reach its equilibrium value.
Luckily, living fish come in units, so we rephrase the question:
How long does it take for the rainbowfish population to reach the last integer before the mathematical equilibrium?
This means that we want to find time teq such that P(teq)=720.
The direction field of P and some possible solutions are shown below.
As you can see, many solutions are possible, each of which reaches 720 at a different time. In an exercise of § 2.1 you have estimated the value of teq by imagining one of the solution curves. As you can see, it is rather easy to guess how quickly P rises in the middle. It is much more difficult to guess how long the beginning and the end phase take.
On the following pages, you will learn a method to approximate the solutions in a computer simulation, and then you can estimate much more precisely how long everything will take.
Calculation (Euler's method)
The differential equation for the rainbowfish that we have now, could still be solved by hand. It would give you the analytical solution, which is exact. However, we are not doing that in this course, as we want to focus on the modelling, and not on analytical solution methods. We also do not want to limit our models to equations that can be solved analytically.
So, we will approximate the solutions of the differential equation with a numerical method, and we choose Euler's method. In the next video, this method is explained.
Euler's method
Calculation (with Python)
Python program for Euler's Method
In the previous section Euler's Method was introduced, and now you can approximate the solution of a differential equation. Of course calculations by hand are inefficient. The method is very suitable to implement in a computer program: much faster and probably more accurate. So in this section, you will use Python to simulate the same solutions.
The basic code for Euler's method is:
# Program : Euler's method # Author : MOOC team Mathematical Modelling Basics # Created : April, 2017 import numpy as np import matplotlib.pyplot as plt print("Solution for dP/dt = 0.7*P") # in Python 2.7: use no brackets # Initializations Dt = 0.1 # timestep Delta t P_init = 10 # initial population t_init = 0 # initial time t_end = 5 # stopping time n_steps = int(round((t_end-t_init)/Dt)) # total number of timesteps t_arr = np.zeros(n_steps + 1) # create an array of zeros for t P_arr = np.zeros(n_steps + 1) # create an array of zeros for P t_arr[0] = t_init # add the initial P to the array P_arr[0] = P_init # add the initial t to the array # Euler's method for i in range (1, n_steps + 1): P = P_arr[i-1] t = t_arr[i-1] dPdt = 0.7*P # calculate the derivative P_arr[i] = P + Dt*dPdt # calculate P on the next time step t_arr[i] = t + Dt # adding the new t-value to the list # Plot the results fig = plt.figure() # create figure plt.plot(t_arr, P_arr, linewidth = 4) # plot population vs. time plt.title('dP/dt = 0.7P, P(0)=10', fontsize = 25) plt.xlabel('t (in days)', fontsize = 20) plt.ylabel('P(t)', fontsize = 20) plt.xticks(fontsize = 15) plt.yticks(fontsize = 15) plt.grid(True) # show grid plt.axis([0, 5, 0, 200]) # define the axes plt.show() # show the plot # save the figure as .jpg fig.savefig('Rainbowfish.jpg', dpi=fig.dpi, bbox_inches = "tight")
Run it, and see what it does!
Basic_Euler.py: header and initialisations
# Program : Euler's Method
# Author : MOOC team Mathematical Modelling Basics
# Created : April, 2017
import numpy as np
import matplotlib.pyplot as plt
print("Solution for dP/dt = 0.7*P")
Dt = 0.1 # timestep Delta t
P_init = 10 # initial population
t_init = 0 # initial time
t_end = 5 # stopping time
n_steps = int(round((t_end-t_init)/Dt)) # total number of timesteps
P_arr = np.zeros(n_steps + 1) # create an array of zeros for P
t_arr = np.zeros(n_steps + 1) # create an array of zeros for t
t_arr[0] = t_init # add the initial t to the array
P_arr[0] = P_init # add the initial P to the array
Basic_Euler.py: Iterations Euler's Method
for i in range (1, n_steps + 1):
P = P_arr[i-1]
t = t_arr[i-1]
dPdt = 0.7*P # calculate the derivative
P_arr[i] = P + Dt*dPdt # calculate P on the next time step
t_arr[i] = t + Dt # adding the new t-value to the list
Basic_Euler.py: Plotting the results
fig = plt.figure() # create figure
plt.plot(t_arr, P_arr, linewidth = 4) # plot population vs. time
plt.title('dP/dt = 0.7P, P(0)=10', fontsize = 25)
plt.xlabel('t (in days)', fontsize = 20)
plt.ylabel('P(t)', fontsize = 20)
plt.xticks(fontsize = 15)
plt.yticks(fontsize = 15)
plt.grid(True) # show grid
plt.axis([0, 5, 0, 200]) # define the axes
plt.show() # show the plot
# save the figure as .jpg
fig.savefig('Rainbowfish.jpg', dpi=fig.dpi, bbox_inches = "tight")
As you have seen in the previous section, a large time step of 7 days leads to results that are not valid for the rainbowfish population. In the next video you will see what happens when the stepsize is reduced.