Refining and Optimizing Analytical Models
Site: | Saylor Academy |
Course: | BUS250: Introduction to Business Intelligence and Analytics |
Book: | Refining and Optimizing Analytical Models |
Printed by: | Guest user |
Date: | Tuesday, July 1, 2025, 7:29 PM |
Description
This section gives an example of how the model developed earlier could be extended to account for additional knowledge of the system's structure. Again, if you have not yet been exposed to this particular type of mathematical analysis, simply note that all models can be refined as we gain more knowledge about the system. Think back to the model you identified in the last section. How could you extend this model to make it a better reflection of reality?
Interaction (Euler's method for Systems)
Rainbowfish and gouramis
In the second module, we modelled the population size of rainbowfish in an large aquarium using limited growth and harvesting.
In this section, we consider the following problem:
A client of the fish farm wants an aquarium containing both rainbowfish and gouramis. The client is advised that the particular kind of gouramis he wants, can get aggressive, and that they hunt and eat rainbowfish. The client is undounted and still wants this combination. He has a budget for 20 rainbowfish, 5 gouramis and an aquarium which has a maximum capacity for 100 rainbowfish. Is it possible to have both fish in the same aquarium without one of them going extinct?
On the next page, we will start with the construction of a mathematical model for the rainbowfish and gouramis.
Source: Marleen Keijzer, Dennis den Ouden-vander Horst, and Kees Vuik, https://ocw.tudelft.nl/courses/mathematical-modeling-basics/ This work is licensed under a Creative Commons Attribution-NonCommercial-ShareAlike 4.0 License.
Interaction Mathematical model
Rainbowfish and gouramis
The problem is the followingA client of the fish farm wants an aquarium containing both rainbowfish and gouramis. The client is advised that the particular kind of gouramis he wants, can get aggressive, and that they hunt and eat rainbowfish. The client is undounted and still wants this combination. He has a budget for 20 rainbowfish, 5 gouramis and an aquarium which has a maximum capacity for 100 rainbowfish. Is it possible to have both fish in the same aquarium without one of them going extinct?
In the next video, a model is constructed for these two populations.
Interacting populations
Interaction Calculation
The client wants rainbowfish and gouramis in the same aquarium. That might not be a good idea, because the gouramis eat rainbowfish. Will it lead to a stable mix of the two species, or will one of them die out?
On the previous page you have seen the construction of the initial value problem for the rainbowfish population P(t) and the gourami population G(t):
\( \begin{cases}\frac{d P}{d t}=0.7 P(t)-0.007 P^2(t)-0.04 P(t) G(t), & P(0)=20, \\ \frac{d G}{d t}=-0.25 G(t)+0.008 P(t) G(t), & G(0)=5 .\end{cases} \)
In the next video you will learn how you can use Euler's Method for a system of differential equations.
Euler's Method for Systems
In the video you have learned about Euler's Method for systems of differential equations.
For the general differential equation
\( \dfrac{d \vec X}{dt} = \vec F(t, \vec X) \),
the n-th step of Euler's method is given by
in which Δt is some time step you have to choose.
Of course, doing calculations by hand quickly becomes laborious. So again we turn to a computer program. The following file contains an example code in Python for a system of differential equations:
Basic_Euler_program_systems.py
First some background information on arrays in Python.
Arrays in Python
For systems of differential equations, it is convenient to combine the different dependent variables in a vector. Here we use
\(\vec{X}(t)=\left[\begin{array}{l}
P(t) \\
G(t)
\end{array}\right] .\)
Especially when there are many dependent variables, a vector notation makes the system often clearer and neater, both in analytical calculations and in numerical simulations.
In Python, you could try to implement vectors as lists. However, lists are not appropriate to calculate with. For example, when you try: (3.0,3.0)/2.0, you will get an error message. The numpy package however, does have arrays with which you can implement vectors and matrices: np.array([3.0,3.0])/2 returns array([1.5,1.5]).
Python's default mode to calculate with vectors and matrices is by element, so beware. For example, matrix
\( A=\left[\begin{array}{l}2 & 3 \\4 & 5\end{array}\right] \)
can be implemented as A = np.array([[2,3],[4,5]]). When you calculate A*A, or equivalently A**2, Python returns array([[ 4, 9], [16, 25]]), so all four elements have been squared, but this is not the product of the two matrices. To calculate the product of two matrices, you need a special command: np.dot(A,A) which returns the correct array([[ 16, 21], [28, 37]]).
Calculation (phase plane)
On the previous page, you have learned to approximate solutions of a system of differential equations with Euler's method. It is the same as for single differential equations, only now you work with a vector instead of a scalar. A graph of the solutions of the initial value problem is the following.
Both P and G have been plotted as functions of time t. Another way to visualise these solutions is to plot P and G against each other: to draw the solutions in the phase plane.
The phase plane is in this case the PG-plane. At any time t the solutions P and G constitute a point (P(t),G(t)) in the phase plane. The collection of points (P(t),G(t)) for parameter t in an interval form a curve in the phase plane.
A solution curve drawn in the phase plane is called a trajectory. If you take another initial value for P and/or G, you obtain a different solution, and you can add its trajectory to the phase plane as well. A few trajectories in a phase plane constitute a phase portrait of the differential equations. An example of a phase portrait is the following graph:
Interaction Validation
Validation
The problem was
A client of the fish farm wants an aquarium with 20 rainbowfish, and 5 gouramis. Is it possible to have both fish in the same aquarium without one of them going extinct?
After constructing a mathematical model (for which the numbers are admittedly soso), graphs of the resulting solutions for both populations are
So the number of rainbowfish will grow fast from the initial 20 to about 55 in the first week. Because of the plentyfull rainbowfish, the number of gouramis will triple in a week and a half. The gouramis will eat more and more of the rainbowfish, so the number of rainbowfish will decrease again to about 25 after two weeks. With fewer rainbowfish, the number of gouramis decreases again, giving the rainbowfish population the chance to grow again. And so on, and so on. The increases and decreases of both populations become smaller in time, and after a month and a half, both populations stabilize to approximately 31 rainbowfish and 12 gouramis.
Practice problem (Clock!)
Problem
Building a clock for a laptop
The clock speed of a computer is set in a base clock. The frequency of the base clock is multiplied in a clock multiplier, and the result is the clock speed at which the calculations are done by the CPU of the computer. We want to build a simple base clock that generates a signal of 100MHz = 108 oscillations per second. The clock multiplier that we have multiplies that frequency with a factor 16, so the CPU can run at 1.6×109 hertz = 1.6GHz.
A simple base clock is an LC-oscillator. The basis of an LC-oscillator is an electrical circuit that contains an inductor, a capacitor, a resistor and a battery.
We have got a 5V battery: VB=5V, and an inductor with an induction of 0.004 microhenry: L=0.004μH. The resistance is estimated as 0.1 ohm: R=0.1Ω. We want to know what capacitor should be used.
Our problem thus becomes
Find the capacitance C of the capacitor such that the current in the electrical circuit oscillates with a frequency of 100 MHz.
Electrical engineers probably can answer this question off the top of their heads, and we admire them for that! The rest of us however needs a mathematical model, and to calculate solutions to answer this question. So, the next page is for the mathematical model for this circuit, and the calculations will hopefully be interesting for the electrical engineers as well.
On the next page, you will derive the mathematical model needed to solve the problem.
Mathematical model
Course subject(s) 3. Extending the model
For a base clock for a computer, we want to build an LC-oscillator that generates 100MHz oscillations.
The components that are given are:
- a 5V battery: VB=5V,
- an inductor: L=0.004μH,
- the resistance: R=0.1Ω.
Our problem is
Find the capacitance C of the capacitor ([C]= henry) such that the current in the electrical circuit oscillates with a frequency of 100 MHz.
The following physical laws relate the current I(t) running in the circuit (measured in ampere), the charge on (one side of) the capacitor Q(t) (in coulomb) and the voltage drops over the different elements (in volts).
The current, I(t), is related to the charge on the capacitor, Q(t), according to I(t)=dQ/dt.
The voltage drop over the capacitor, V(t), is related to the charge on the capacitor by: V(t)=Q(t)/C.
The voltage drop over the inductor, VL(t), is related to the change in the current: VL(t)=L(dI/dt).
And an expression for voltage drop over the resistor, VR(t) is given by: VR(t)=RI(t).
Kirchhoff's voltage law states that the voltage drop around a loop in a circuit equals zero: V+VL+VR−VB=0.
There are different options for a set of dependent variables to calculate in. We choose to calculate in current I(t) and the voltage drop over the capacitor V(t).
Calculation
The system of differential equation for the LC-oscillator is:
\(c\dfrac{dV}{dt} = I(t)\), \(V(0) = 0\),
\(L \dfrac{dI}{dt} = V_B - RI (t) - (t) \), \(I(0) =0\).
Constants VB,L and R are given: VB=5V, L=0.004μH, R=0.1Ω, and C has to be found such that the circuit oscillates with at frequency of 100MHz =108Hz.
First, we need to choose a time scale, and then we have to scale the problem and choose an appropriate step size.
Validation
Validation
On the last page, you have calculated for what values this LC-oscillator will oscillate with a frequency of 100MHz:
VB=5V, L=0.004μH, R=0.1Ω, and C=0.633nF.
This choosing of the right capacitance for the condensator is called tuning. In older radios this is done mechanically with a variable condensator, also called a tuning condensator.
Rotary capacitor by Zátonyi Sándor CC BY 3.0
As you can see in the graph, the amplitude of the signal dies out, because however small, there is some energy loss because R≠0 (and in reality also some charge will leak from the capacitor). So this basic LC-oscillator will only work for say 500nanoseconds or so, which is not very long! To maintain the oscillation for longer times, other components (e.g. a transistor or an Opamp) are built into the circuit, to give the current little pushes at the right times, and keep the amplitude constant.
Oscillator circuits are not only used to set the clock speed in computers and electronic clocks and to tune radio and televisions, but also to generate those radio and television waves, and to generate the buzz of a door bell and the electronic beeps in video games.
And of course, modern oscillators are way more complicated than the one shown here. For example, often a quartz christal is cleverly used to keep the oscillations generated by the circuit really constant over longer times. Many of the components used nowadays in circuits are non-linear, and their mathematical models contain non-linear terms and the resulting waves are often far from nice, periodic sines or cosines. Very interesting to model!
Changing the resistance R
Untill now, we have kept the resistance R=0.1Ω. In the next exercises, let's investigate what happens when R is changed
Long-term behaviour (Equilibrium, Phase space)
Rainbowfish and gouramis
In section Interaction (Euler's method for Systems), you encountered the following situation:
A client of the fish farm wants an aquarium containing both
rainbowfish and gouramis. He has a budget for 20 rainbowfish, 5 gouramis
and an aquarium which has a maximum capacity for 100 rainbowfish.
A model was constructed and solutions were calculated with Euler's method. According to the model, both populations start out to grow fast to triple their initial sizes. Then first the number of rainbowfish decreases again, and the gouramis follow 6 days later. Both populations keep oscillating in a damped oscillation, and in a month or so both populations have settled in their end states of approximately 31 and 12 fishes respectively.
In this section we are going to investigate this problem further to understand the model better.
What would happen with different initial numbers of both kinds of fish? There are three equilibrium points. One is the end state. What is going to happen with the populations when we start near one of the other equilibrium points? Can we calculate some results analytically, so we can understand better how the parameters we have chosen for the model influence the results? And maybe then we have more information to chose realistic values for the parameters.
The mathematical model for the populations of rainbowfish and gourami is the same as in Interaction (Euler's method for Systems), so on the next page we skip immediately to the calculation part of the modelling cycle.
Calculation (linearisation)
A mathematical model for interacting populations of rainbowfish (P(t)) and gourami (G(t)) in a 100-rainbowfish aquarium is the following system of differential equations.
\(\begin{cases}\frac{d P}{d t}=0.7 P(t)-0.007 P^2(t)-0.04 P(t) G(t), & P(0)=20, \\ \frac{d G}{d t}=-0.25 G(t)+0.008 P(t) G(t), & G(0)=5 .\end{cases}\)
You have learned to simulate solutions with Euler’s method and to visualise these in the phase plane. Here you see a combination of trajectories, each with a different set of initial values.
Can we calculate some results analytically, so we can understand better how the parameters we have chosen for the model influence the results?
In this section you will learn how you to linearise the system in the neighbourhood of an equilibrium point. With the linearisations, you can predict some of the results analytically from the parameters.
The definition of an equilibrium point is:
An equilibrium point of a system of differential equations
\(\dfrac{d \vec X}{dt} = \vec F(t, \vec X)\)
is a point X→0 where
\(\dfrac{d \vec X_0}{dt} = \vec F(t, \vec X_0) = \vec 0\).
In §§3.1 Calculation (phase plane), you have calculated the equilibrium points for this system: (0,0), (100,0) and (31.25,12.03125). Adding these equilibrium points to the phase portrait gives:
You can see special things happening near the equilibrium points. For example near the point (100,0), you can see lines going first towards and then away from the equilibrium. Could we have predicted this behaviour using just the differential equations?
The answer is "Yes", and in the remainder of this section you will learn how you can do this.
The first step in obtaining a prediction of the behavior near the equilibrium point is the linearisation of the differential equations around the equilibrium point.
The linearisation of a function f(P,G) of two variables P and G is the process of approximating this function f in the neighbourhood of the point (P0,G0) by a function L(P,G) with the formula:
\(L(P, G)=f\left(P_0, G_0\right)+\left.\frac{\partial f}{\partial P}\right|_{\left(P_0, G_0\right)}\left(P-P_0\right)+\left.\frac{\partial f}{\partial G}\right|_{\left(P_0, G_0\right)}\left(G-G_0\right) .\)
is the plane tangent to f in the point (P0,G0,f(P0,G0)). A visualisation of such a tangent plane is:
In our case we have two differential equations, and two functions to linearise. In a general form:
\( \begin{cases}\frac{d P}{d t}=f_1 (P, G), \\ \frac{d G}{d t}= f_2 (P, G).\end{cases} \)
We can linearise the two right-hand-side functions f1 and f2 around the same equilibrium point (P0,G0), which leads to the general system
\(\left\{\begin{array}{l}
\frac{d P}{d t}=\left.\frac{\partial f_1}{\partial P}\right|_{\left(P_0, G_0\right)}\left(P-P_0\right)+\left.\frac{\partial f_1}{\partial G}\right|_{\left(P_0, G_0\right)}\left(G-G_0\right), \\
\frac{d G}{d t}=\left.\frac{\partial f_2}{\partial P}\right|_{\left(P_0, G_0\right)}\left(P-P_0\right)+\left.\frac{\partial f_2}{\partial G}\right|_{\left(P_0, G_0\right)}\left(G-G_0\right) .
\end{array}\right.\)
Note that the terms f1(P0,G0) and f2(P0,G0) have been left out, because both are zero in an equilibrium point by definition.
Calculation (saddle points and nodes)
On the last page, some new notation was introduced:
\( \vec M (t) = \vec X (t) - \vec X_0 =\left[\begin{array}{l}P(t) \\G(t)\end{array}\right] -\left[\begin{array}{l}P_0 \\G_0\end{array}\right] \),
where (P0,G0) is an equilibrium point. You have learned that a system of differential equations that is linearised around an equilibrium point, can be written as:
\( \dfrac{d \vec M}{dt} = J (\vec X_0) \vec M (t) \).
When our rainbowfish/gourami system is linearised around (P0,G0)=(100,0), the result is:
\( \dfrac{d \vec M}{dt} \left[\begin{array}{l}-0.7 & -4) \\0 & 0.55\end{array}\right] \vec M(t) \).
In the next video, Dennis explains what this new differential equation can tell you about the original differential equations.
Equilibrium points in a system
In the video you have seen that (100,0) is a saddle point by looking at the solutions of the linearised differential equation. The behaviour of solutions near a saddle point is explained by the eigenvalues of the Jacobian matrix: one is positive, and one is negative.
Of course, the eigenvalues of a 2×2-matrix can also be both negative or both positive. Then both factors eλ1t and eλ2t will either both decrease in time (when λ1<0 and λ2<0) or both increase in time, (when λ1>0 and λ2>0). The equilibrium points are then called nodes.
An equilibrium point X→0 is called a saddle point if the Jacobian matrix J(X→0) has one negative and one positive eigenvalue. A saddle point is unstable because some of the solutions that start near the equilibrium point (here the origin) leave the neighborhood of the origin. A typical sketch of the solutions near a saddle point in the phase plane is given by
An equilibrium point X→0 is called a stable node if the Jacobian matrix J(X→0) has two negative eigenvalues: all solutions that start near the equilibrium point stay near the equilirium point. An example of a phase portrait is
An equilibrium point X→0 is called an unstable node if the Jacobian matrix J(X→0) has two positive eigenvalues. A typical sketch of the solutions near an unstable node in the phase plane is given by
Calculation (spiral points)
On the previous page you learned about saddle points, stable nodes and unstable nodes. On this page you will encounter even more types of equilibrium points.
As you have seen in the last exercise, the eigenvalues of the Jacobian matrix in the point (31.25,12.03) are complex. Complex numbers are not positive or negative, so this equilibrium point is neither a stable node, nor an unstable node nor a saddle point. So what is it then?
In the phase portrait, you see that the solutions near the equilibrium (31.25,12.03) spiral around it. Can we explain this with the complex eigenvalues you found?
Consider the general form of the solution of \( \vec{M} \) from the video:
\(\vec{M}(t)=c_1 e^{\lambda_1 t} \vec{v}_1+c_2 e^{\lambda_2 t} \vec{v}_2\)
The λ's are complex numbers and the solution contains factors \(e^{λt}\). For exponential functions of complex numbers, we can use Euler's formula
\( e^{i \phi} = \cos (\phi) + i \; \sin (\phi) \)
Take the first eigenvalue λ1, which is complex of the form λ1=a+bi. Then
\(e^{\lambda_1 t}=e^{a t+i b t}=e^{a t}(\cos (b t)+i \sin (b t))\)
(\( \lambda_2 = \overline \lambda_1 = a - bi \) gives a similar expression). The solution \( \vec{M}(t) \) is now an expression with the complex i in it many times. Of course, \(\vec{M}\) denotes a vector with (shifted) population sizes, so it should be real. In this course, we will not go into the details of how you can rewrite the solution for \(\vec{M}\) to a real expression, but just state the result:
When the eigenvalues of the Jacobian matrix are complex numbers: λ=a±bi,
the solutions for \(\vec{M} (t)\) contain factors \(e^{at} \cos (bt)\) and \(e^{at} \sin (bt)\).
Looking at these factors, we see that each component of \(\vec{M} (t)\) has oscillating parts: either \(\cos (bt)\), or \(\sin (bt)\) or both. So all terms oscillate with an angular frequency of b radians per day, where b is the imaginary part of the eigenvalues. The graphs of the solutions versus time oscillate up and down. In the phase portrait the trajectories will circle the origin.
Each component of \(\vec{M} (t)\) also contains a factor \(e^{at}\). If a < 0 that factor decreases in time, so the oscillations will be damped and in phase space the trajectories will spiral in towards the origin. If a > 0, the exponential factor will increase in time, making the oscillations larger and larger, and in the phase plane the trajectories will spiral outwards.
This type of behavior around an equilibrium point makes it a spiral point. When the real part a of the eigenvalue is negative, so a < 0, all solutions that start near the origin will end up in the origin, so that spiral point is a stable equilibrium point. When a > 0, solutions that start near the origin will spiral away from the origin, so then the origin is called an unstable spiral point.
An equilibrium point \(\vec{X}_0\) is called a stable spiral point if the Jacobian matrix \(J (\vec{X}_0)\) has two complex eigenvalues \( \lambda = a \pm bi \) with negative real parts: a < 0. A typical sketch of the solutions near a stable spiral point in the phase plane is given by
An equilibrium point \(\vec{X}_0\) is called an unstable spiral point if the Jacobian matrix \((\vec{X}_0)\) has two complex eigenvalues \( \lambda = a \pm bi \) with positive real parts: a > 0. A phase portrait of an unstable spiral point is given in the following graph.
When a=0, the solution does not spiral inwards nor outwards, so the solution will be periodic, and in phase space the trajectory is a closed curve. Even though that curve is most often an ellips, the equilibrium point is called a circle point.
An equilibrium point \(\vec{X}_0\) is called a circle point if the Jacobian matrix \(J(\vec{X}_0)\) has two complex eigenvalues with zero real parts \( \lambda = a \pm bi \):
Besides the six types of stability for equilibrium points discussed here, other types of stability exist, but these transcend the scope of this course.
Validation
We have worked with the system of two differential equations for the rainbowfish and gouramis.
\( \begin{cases}\frac{d P}{d t}=0.7 P-0.007 P^2- \alpha PG, & \text {where α = 0.04} , \\ \frac{d G}{d t}=-0.25 G+ \beta PG, & \text {where β = 0.008}.\end{cases} \)
You have learned how to approximate solutions numerically with Euler's method. In this section, we have used analytical methods to characterize the solutions beforehand. First the equilibrium points (P0,G0) were calculated:
\((0,0), (100, 0)\) and \( (\dfrac{1}{4 \beta} , \dfrac{0.7}{\alpha} (1 - \dfrac{1}{400 \beta}) ) \).
Then the differential equations were linearised around the three equilibrium points:
\(\begin{gathered}
\text { For } \vec{M}(t)=\left[\begin{array}{c}
P(t)-P_0 \\
G(t)-G_0
\end{array}\right], \\
\frac{d \vec{M}}{d t} \approx\left[\begin{array}{cc}
0.7-0.014 P_0-\alpha G_0 & -\alpha P_0 \\
\beta G_0 & -0.25+\beta P_0
\end{array}\right] \vec{M}(t) .
\end{gathered}\)
By calculating the eigenvalues of the three matrices, the types of the three equilibrium points were determined. In the following exercises, you are going to review whether these analytical calculations were usefull.
Practice problem (Clock!)
Problem
Building a clock for a computer
In the previous section we wanted to build a base clock for a computer, an LC-oscillator that generates 100MHz oscillations. Using Euler’s method you have found an approximate value for the capacitance of the capacitor. And after that, the value of the resistance was changed, yielding other types of solutions.
In this module, you do the same calculations, not numerically, but analytically. This will give you more insight in oscillating solutions.
Mathematical model
In Section 3.2 the mathematical model for the circuit was derived. I(t) is the current (in ampere), V(t) is the voltage drop over the capacitor (in volt). The derivatives and the constants were scaled, so the numbers for the simulations were of a reasonable size. For analytical calculations, this scaling is not necessary, so we will use the original variables.
\(\dfrac{dV}{dt}=\dfrac{I}{C}\), \(V(0) = 0\),
\(\dfrac{dI}{dt} = \dfrac{-V}{L} - \dfrac{RI}{L} + \dfrac{V_B}{L}\), \(I(0) = 0\),
where \(V_B=5V\), \(R=0.1Ω\), \(L=4×10^{−9}H\), and for the capacitance, we start with \(C=0.5×10^{−9}F\).
On the next page you will investigate the equilibrium point and its type.
Calculation
The mathematical model is
\(\dfrac{dV}{dt}=\dfrac{I}{C}\), \(V(0) = 0\),
\(\dfrac{dI}{dt} = \dfrac{-V}{L} - \dfrac{RI}{L} + \dfrac{V_B}{L}\), \(I(0) = 0\),
where \(V_B=5V\), \(R=0.1Ω\), \(L=4×10^{−9}H\), and \(C=0.5×10^{−9}F\).
Validation
On the last page, you have calculated for what values this LC-oscillator will oscillate with a frequency of 100 MHZ:
\(V_B = 5V\), \(L = 4nH\), \( R = 0.1 \Omega \), and \(C = 0.633 nF\)
In the next exercises you change the resistance R and calculate what happens analytically.