Figure 17: Illustration of the compaction of geological layers (with different colors) through time.
This section presents many mathematical models that all end up with ODEs of the type \( u'=-au+b \). The applications are taken from biology, finance, and physics, and cover population growth or decay, compound interest and inflation, radioactive decay, cooling of objects, compaction of geological media, pressure variations in the atmosphere, and air resistance on falling or rising bodies.
Real applications of a model \( u'=-au+b \) will often involve a lot of parameters in the expressions for \( a \) and \( b \). It can be quite a challenge to find relevant values of all parameters. In simple problems, however, it turns out that it is not always necessary to estimate all parameters because we can lump them into one or a few dimensionless numbers by using a very attractive technique called scaling. It simply means to stretch the \( u \) and \( t \) axis is the present problem - and suddenly all parameters in the problem are lumped one parameter if \( b\neq 0 \) and no parameter when \( b=0 \)!
Scaling means that we introduce a new function \( \bar u(\bar t) \), with $$ \bar u = \frac{u - u_m}{u_c},\quad \bar t = \frac{t}{t_c},$$ where \( u_m \) is a characteristic value of \( u \), \( u_c \) is a characteristic size of the range of \( u \) values, and \( t_c \) is a characteristic size of the range of \( t_c \) where \( u \) varies significantly. Choosing \( u_m \), \( u_c \), and \( t_c \) is not always easy and often an art in complicated problems. We just state one choice first: $$ u_c = I,\quad u_m = b/a,\quad t_c = 1/a\tp$$ Inserting \( u=u_m + u_c\bar u \) and \( t=t_c\bar t \) in the problem \( u'=-au + b \), assuming \( a \) and \( b \) are constants, results after some algebra in the scaled problem $$ \frac{d\bar u}{d\bar t} = -\bar u,\quad \bar u(0)=1 - \beta,$$ where \( \beta \) is a dimensionless number $$ \beta = \frac{b}{Ia}\tp$$ That is, only the special combination of \( b/(Ia) \) matters, not what the individual values of \( b \), \( a \), and \( I \) are. Moreover, if \( b=0 \), the scaled problem is independent of \( a \) and \( I \)! In practice this means that we can perform one numerical simulation of the scaled problem and recover the solution of any problem for a given \( a \) and \( I \) by stretching the axis in the plot: \( u=I\bar u \) and \( t =\bar t/a \). For \( b\neq 0 \), we simulate the scaled problem for a few \( \beta \) values and recover the physical solution \( u \) by translating and stretching the \( u \) axis and stretching the \( t \) axis.
The scaling breaks down if \( I=0 \). In that case we may choose \( u_m=0 \), \( u_c=b/a \), and \( t_c=1/b \), resulting in a slightly different scaled problem: $$ \frac{d\bar u}{d\bar t} = 1 -\bar u,\quad \bar u(0)=0\tp$$ As with \( b=0 \), the case \( I=0 \) has a scaled problem with no physical parameters!
It is common to drop the bars after scaling and write the scaled problem as \( u'=-u \), \( u(0)=1-\beta \), or \( u'=1-u \), \( u(0)=0 \). Any implementation of the problem \( u'=-au+b \), \( u(0)=I \), can be reused for the scaled problem by setting \( a=1 \), \( b=0 \), and \( I=1-\beta \) in the code, if \( I\neq 0 \), or one sets \( a=1 \), \( b=1 \), and \( I=0 \) when the physical \( I \) is zero. Falling bodies in fluids, as described in the section Vertical motion of a body in a viscous fluid, involves \( u'=-au+b \) with seven physical parameters. All these vanish in the scaled version of the problem if we start the motion from rest!
Let \( N \) be the number of individuals in a population occupying some spatial domain. Despite \( N \) being an integer in this problem, we shall compute with \( N \) as a real number and view \( N(t) \) as a continuous function of time. The basic model assumption is that in a time interval \( \Delta t \) the number of newcomers to the populations (newborns) is proportional to \( N \), with proportionality constant \( \bar b \). The amount of newcomers will increase the population and result in to $$ N(t+\Delta t) = N(t) + \bar bN(t)\tp $$ It is obvious that a long time interval \( \Delta t \) will result in more newcomers and hence a larger \( \bar b \). Therefore, we introduce \( b=\bar b/\Delta t \): the number of newcomers per unit time and per individual. We must then multiply \( b \) by the length of the time interval considered and by the population size to get the total number of new individuals, \( b\Delta t N \).
If the number of removals from the population (deaths) is also proportional to \( N \), with proportionality constant \( d\Delta t \), the population evolves according to $$ N(t+\Delta t) = N(t) + b\Delta t N(t) - d\Delta t N(t)\tp $$ Dividing by \( \Delta t \) and letting \( \Delta t \rightarrow 0 \), we get the ODE $$ \begin{equation} N' = (b-d)N,\quad N(0)=N_0\tp \end{equation} $$ In a population where the death rate (\( d \)) is larger than then newborn rate (\( b \)), \( a>0 \), and the population experiences exponential decay rather than exponential growth.
In some populations there is an immigration of individuals into the spatial domain. With \( I \) individuals coming in per time unit, the equation for the population change becomes $$ N(t+\Delta t) = N(t) + b\Delta t N(t) - d\Delta t N(t) + \Delta t I\tp $$ The corresponding ODE reads $$ \begin{equation} N' = (b-d)N + I,\quad N(0)=N_0 \tp \end{equation} $$
Some simplification arises if we introduce a fractional measure of the population: \( u=N/N_0 \) and set \( r=b-d \). The ODE problem now becomes $$ \begin{equation} u' = ru + f,\quad u(0)=1, \tag{65} \end{equation} $$ where \( f=I/N_0 \) measures the net immigration per time unit as the fraction of the initial population. Very often, \( r \) is approximately constant, but \( f \) is usually a function of time.
The growth rate \( r \) of a population decreases if the environment has limited resources. Suppose the environment can sustain at most \( N_{\max} \) individuals. We may then assume that the growth rate approaches zero as \( N \) approaches \( N_{\max} \), i.e., as \( u \) approaches \( M=N_{\max}/N_0 \). The simplest possible evolution of \( r \) is then a linear function: \( r(t)=r_0(1-u(t)/M) \), where \( r_0 \) is the initial growth rate when the population is small relative to the maximum size and there is enough resources. Using this \( r(t) \) in (65) results in the logistic model for the evolution of a population (assuming for the moment that \( f=0 \)): $$ \begin{equation} u' = r_0(1-u/M)u,\quad u(0)=1 \tp \tag{66} \end{equation} $$ Initially, \( u \) will grow at rate \( r_0 \), but the growth will decay as \( u \) approaches \( M \), and then there is no more change in \( u \), causing \( u\rightarrow M \) as \( t\rightarrow\infty \). Note that the logistic equation \( u'=r_0(1-u/M)u \) is nonlinear because of the quadratic term \( -u^2r_0/M \).
Say the annual interest rate is \( r \) percent and that the bank adds the interest once a year to your investment. If \( u^n \) is the investment in year \( n \), the investment in year \( u^{n+1} \) grows to $$ u^{n+1} = u^n + \frac{r}{100}u^n \tp $$ In reality, the interest rate is added every day. We therefore introduce a parameter \( m \) for the number of periods per year when the interest is added. If \( n \) counts the periods, we have the fundamental model for compound interest: $$ \begin{equation} u^{n+1} = u^n + \frac{r}{100 m}u^n \tp \tag{67} \end{equation} $$ This model is a difference equation, but it can be transformed to a continuous differential equation through a limit process. The first step is to derive a formula for the growth of the investment over a time \( t \). Starting with an investment \( u^0 \), and assuming that \( r \) is constant in time, we get $$ \begin{align*} u^{n+1} &= \left(1 + \frac{r}{100 m}\right)u^{n}\\ &= \left(1 + \frac{r}{100 m}\right)^2u^{n-1}\\ &\ \ \vdots\\ &= \left(1 +\frac{r}{100 m}\right)^{n+1}u^{0} \end{align*} $$ Introducing time \( t \), which here is a real-numbered counter for years, we have that \( n=mt \), so we can write $$ u^{mt} = \left(1 + \frac{r}{100 m}\right)^{mt} u^0\tp $$ The second step is to assume continuous compounding, meaning that the interest is added continuously. This implies \( m\rightarrow\infty \), and in the limit one gets the formula $$ \begin{equation} u(t) = u_0e^{rt/100}, \end{equation} $$ which is nothing but the solution of the ODE problem $$ \begin{equation} u' = \frac{r}{100}u,\quad u(0)=u_0 \tp \tag{68} \end{equation} $$ This is then taken as the ODE model for compound interest if \( r>0 \). However, the reasoning applies equally well to inflation, which is just the case \( r < 0 \). One may also take the \( r \) in (68) as the net growth of an investemt, where \( r \) takes both compound interest and inflation into account. Note that for real applications we must use a time-dependent \( r \) in (68).
Introducing \( a=\frac{r}{100} \), continuous inflation of an initial fortune \( I \) is then a process exhibiting exponential decay according to $$ u' = -au,\quad u(0)=I\tp $$
An atomic nucleus of an unstable atom may lose energy by emitting ionizing particles and thereby be transformed to a nucleus with a different number of protons and neutrons. This process is known as radioactive decay. Actually, the process is stochastic when viewed for a single atom, because it is impossible to predict exactly when a particular atom emits a particle. Nevertheless, with a large number of atoms, \( N \), one may view the process as deterministic and compute the mean behavior of the decay. Below we reason intuitively about an ODE for the mean behavior. Thereafter, we show mathematically that a detailed stochastic model for single atoms leads the same mean behavior.
Suppose at time \( t \), the number of the original atom type is \( N(t) \). A basic model assumption is that the transformation of the atoms of the original type in a small time interval \( \Delta t \) is proportional to \( N \), so that $$ N(t+\Delta t) = N(t) - a\Delta t N(t),$$ where \( a>0 \) is a constant. Introducing \( u=N(t)/N(0) \), dividing by \( \Delta t \) and letting \( \Delta t\rightarrow 0 \) gives the following ODE: $$ \begin{equation} u' = -au,\quad u(0)=1 \tp \end{equation} $$ The parameter \( a \) can for a given nucleus be expressed through the half-life \( t_{1/2} \), which is the time taken for the decay to reduce the initial amount by one half, i.e., \( u(t_{1/2}) = 0.5 \). With \( u(t)=e^{-at} \), we get \( t_{1/2}=a^{-1}\ln 2 \) or \( a=\ln 2/t_{1/2} \).
We have originally \( N_0 \) atoms. Each atom may have decayed or survived at a particular time \( t \). We want to count how many original atoms that are left, i.e., how many atoms that have survived. The survival of a single atom at time \( t \) is a random event. Since there are only two outcomes, survival or decay, we have a Bernoulli trial. Let \( p \) be the probability of survival (implying that the probability of decay is \( 1-p \)). If each atom survives independently of the others, and the probability of survival is the same for every atom, we have \( N_0 \) statistically Bernoulli trials, known as a binomial experiment from probability theory. The probability \( P(N) \) that \( N \) out of the \( N_0 \) atoms have survived at time \( t \) is then given by the famous binomial distribution $$ P(N) = \frac{N_0!}{N! (N_0-N)!}p^N (1-p)^{N_0-N}\tp $$ The mean (or expected) value \( \E{P} \) of \( P(N) \) is known to be \( N_0p \).
It remains to estimate \( p \). Let the interval \( [0,t] \) be divided into \( m \) small subintervals of length \( \Delta t \). We make the assumption that the probability of decay of a single atom in an interval of length \( \Delta t \) is \( \tilde p \), and that this probability is proportional to \( \Delta t \): \( \tilde p = \lambda\Delta t \) (it sounds natural that the probability of decay increases with \( \Delta t \)). The corresponding probability of survival is \( 1-\lambda\Delta t \). Believing that \( \lambda \) is independent of time, we have, for each interval of length \( \Delta t \), a Bernoulli trial: the atom either survives or decays in that interval. Now, \( p \) should be the probability that the atom survives in all the intervals, i.e., that we have \( m \) successful Bernoulli trials in a row and therefore $$ p = (1-\lambda\Delta t)^m\tp$$ The expected number of atoms of the original type at time \( t \) is $$ \begin{equation} \E{P} = N_0p = N_0(1-\lambda\Delta t)^m,\quad m=t/\Delta t\tp \end{equation} $$
To see the relation between the two types of Bernoulli trials and the ODE above, we go to the limit \( \Delta t\rightarrow t \), \( m\rightarrow\infty \). One can show that $$ p = \lim_{m\rightarrow\infty} (1-\lambda\Delta t)^m = \lim_{m\rightarrow\infty} \left(1-\lambda\frac{t}{m}\right)^m = e^{-\lambda t} $$ This is the famous exponential waiting time (or arrival time) distribution for a Poisson process in probability theory (obtained here, as often done, as the limit of a binomial experiment). The probability of decay, \( 1-e^{-\lambda t} \), follows an exponential distribution. The limit means that \( m \) is very large, hence \( \Delta t \) is very small, and \( \tilde p=\lambda\Delta t \) is very small since the intensity of the events, \( \lambda \), is assumed finite. This situation corresponds to a very small probability that an atom will decay in a very short time interval, which is a reasonable model. The same model occurs in lots of different applications, e.g., when waiting for a taxi, or when finding defects along a rope.
With \( p=e^{-\lambda t} \) we get the expected number of original atoms at \( t \) as \( N_0p=N_0e^{-\lambda t} \), which is exactly the solution of the ODE model \( N'=-\lambda N \). This gives also an interpretation of \( a \) via \( \lambda \) or vice versa. Our important finding here is that the ODE model captures the mean behavior of the underlying stochastic model. This is, however, not always the common relation between microscopic stochastic models and macroscopic "averaged" models.
Also of interest is to see that a Forward Euler discretization of \( N'=-\lambda N \), \( N(0)=N_0 \), gives \( N^m = N_0(1-\lambda\Delta t)^m \) at time \( t_m=m\Delta t \), which is exactly the expected value of the stochastic experiment with \( N_0 \) atoms and \( m \) small intervals of length \( \Delta t \), where each atom can decay with probability \( \lambda\Delta t \) in an interval.
A fundamental question is how accurate the ODE model is. The underlying stochastic model fluctuates around its expected value. A measure of the fluctuations is the standard deviation of the binomial experiment with \( N_0 \) atoms, which can be shown to be \( \Std{P}=\sqrt{N_0p(1-p)} \). Compared to the size of the expectation, we get the normalized standard deviation $$ \frac{\sqrt{\Var{P}}}{\E{P}} = N_0^{-1/2}\sqrt{p^{-1}-1} = N_0^{-1/2}\sqrt{(1-e^{-\lambda t})^{-1}-1}\approx (N_0\lambda t)^{-1/2}, $$ showing that the normalized fluctuations are very small if \( N_0 \) is very large, which is usually the case.
When a body at some temperature is placed in a cooling environment, experience shows that the temperature falls rapidly in the beginning, and then the changes in temperature levels off until the body's temperature equals that of the surroundings. Newton carried out some experiments on cooling hot iron and found that the temperature evolved as a "geometric progression at times in arithmetic progression", meaning that the temperature decayed exponentially. Later, this result was formulated as a differential equation: the rate of change of the temperature in a body is proportional to the temperature difference between the body and its surroundings. This statement is known as Newton's law of cooling, which can be mathematically expressed as $$ \begin{equation} {dT\over dt} = -k(T-T_s), \tag{80} \end{equation} $$ where \( T \) is the temperature of the body, \( T_s \) is the temperature of the surroundings, which may be time-dependent, \( t \) is time, and \( k \) is a positive constant. Equation (80) is primarily viewed as an empirical law, valid when heat is efficiently convected away from the surface of the body by a flowing fluid such as air at constant temperature \( T_s \). The heat transfer coefficient \( k \) reflects the transfer of heat from the body to the surroundings and must be determined from physical experiments.
The cooling law (80) needs an initial condition \( T(0)=T_0 \).
Vertical equilibrium of air in the atmosphere is governed by the equation $$ \begin{equation} \frac{dp}{dz} = -\varrho g \tp \tag{70} \end{equation} $$ Here, \( p(z) \) is the air pressure, \( \varrho \) is the density of air, and \( g=9.807\hbox{ m/s}^2 \) is a standard value of the acceleration of gravity. (Equation (70) follows directly from the general Navier-Stokes equations for fluid motion, with the assumption that the air does not move.)
The pressure is related to density and temperature through the ideal gas law $$ \begin{equation} \varrho = \frac{Mp}{R^*T}, \tag{71} \end{equation} $$ where \( M \) is the molar mass of the Earth's air (0.029 kg/mol), \( R^* \) is the universal gas constant (\( 8.314 \) Nm/(mol K)), and \( T \) is the temperature. All variables \( p \), \( \varrho \), and \( T \) vary with the height \( z \). Inserting (71) in (70) results in an ODE with a variable coefficient: $$ \begin{equation} \frac{dp}{dz} = -\frac{Mg}{R^*T(z)} p \tag{72} \thinspace . \end{equation} $$
The atmosphere can be approximately modeled by seven layers. In each layer, (72) is applied with a linear temperature of the form $$ T(z) = \bar T_i + L_i(z-h_i),$$ where \( z=h_i \) denotes the bottom of layer number \( i \), having temperature \( \bar T_i \), and \( L_i \) is a constant in layer number \( i \). The table below lists \( h_i \) (m), \( \bar T_i \) (K), and \( L_i \) (K/m) for the layers \( i=0,\ldots,6 \).
| \( i \) | \( h_i \) | \( \bar T_i \) | \( L_i \) |
| 0 | 0 | 288 | -0.0065 |
| 1 | 11,000 | 216 | 0.0 |
| 2 | 20,000 | 216 | 0.001 |
| 3 | 32,000 | 228 | 0.0028 |
| 4 | 47,000 | 270 | 0.0 |
| 5 | 51,000 | 270 | -0.0028 |
| 6 | 71,000 | 214 | -0.002 |
For implementation it might be convenient to write (72) on the form $$ \begin{equation} \frac{dp}{dz} = -\frac{Mg}{R^*(\bar T(z) + L(z)(z-h(z)))} p, \end{equation} $$ where \( \bar T(z) \), \( L(z) \), and \( h(z) \) are piecewise constant functions with values given in the table. The value of the pressure at the sea level \( z=0 \), \( p_0=p(0) \), is \( 101325 \) Pa.
One commonly used simplification is to assume that the temperature is constant within each layer. This means that \( L=0 \).
Another commonly used approximation is to work with one layer instead of seven. This one-layer model is based on \( T(z)=T_0 - Lz \), with sea level standard temperature \( T_0=288 \) K and temperature lapse rate \( L=0.0065 \) K/m.
Sediments, originally made from materials like sand and mud, get compacted through geological time by the weight of new material that is deposited on the sea bottom. The porosity \( \phi \) of the sediments tells how much void (fluid) space there is between the sand and mud grains. The porosity reduces with depth because the weight of the sediments above and causes the void space to shrink and thereby increase the compaction.
A typical assumption is that the change in \( \phi \) at some depth \( z \) is negatively proportional to \( \phi \). This assumption leads to the differential equation problem $$ \begin{equation} \frac{d\phi}{dz} = -c\phi,\quad \phi(0)=\phi_0, \tag{73} \end{equation} $$ where the \( z \) axis points downwards, \( z=0 \) is the surface with known porosity, and \( c>0 \) is a constant.
The upper part of the Earth's crust consists of many geological layers stacked on top of each other, as indicated in Figure 17. The model (73) can be applied for each layer. In layer number \( i \), we have the unknown porosity function \( \phi_i(z) \) fulfilling \( \phi_i'(z)=-c_iz \), since the constant \( c \) in the model (73) depends on the type of sediment in the layer. From the figure we see that new layers of sediments are deposited on top of older ones as time progresses. The compaction, as measured by \( \phi \), is rapid in the beginning and then decreases (exponentially) with depth in each layer.
Figure 17: Illustration of the compaction of geological layers (with different colors) through time.
When we drill a well at present time through the right-most column of sediments in Figure 17, we can measure the thickness of the sediment in (say) the bottom layer. Let \( L_1 \) be this thickness. Assuming that the volume of sediment remains constant through time, we have that the initial volume, \( \int_0^{L_{1,0}} \phi_1 dz \), must equal the volume seen today, \( \int_{\ell-L_1}^{\ell}\phi_1 dz \), where \( \ell \) is the depth of the bottom of the sediment in the present day configuration. After having solved for \( \phi_1 \) as a function of \( z \), we can then find the original thickness \( L_{1,0} \) of the sediment from the equation $$ \int_0^{L_{1,0}} \phi_1 dz = \int_{\ell-L_1}^{\ell}\phi_1 dz \tp $$ In hydrocarbon exploration it is important to know \( L_{1,0} \) and the compaction history of the various layers of sediments.
A body moving vertically through a fluid (liquid or gas) is subject to three different types of forces: the gravity force, the drag force, and the buoyancy force.
The gravity force is \( F_g= -mg \), where \( m \) is the mass of the body and \( g \) is the acceleration of gravity. The uplift or buoyancy force ("Archimedes force") is \( F_b = \varrho gV \), where \( \varrho \) is the density of the fluid and \( V \) is the volume of the body. Forces and other quantities are taken as positive in the upward direction.
The drag force is of two types, depending on the Reynolds number $$ \begin{equation} \hbox{Re} = \frac{\varrho d|v|}{\mu}, \end{equation} $$ where \( d \) is the diameter of the body in the direction perpendicular to the flow, \( v \) is the velocity of the body, and \( \mu \) is the dynamic viscosity of the fluid. When \( \hbox{Re} < 1 \), the drag force is fairly well modeled by the so-called Stokes' drag, which for a spherical body of diameter \( d \) reads $$ \begin{equation} F_d^{(S)} = - 3\pi d\mu v \tp \end{equation} $$ For large Re, typically \( \hbox{Re} > 10^3 \), the drag force is quadratic in the velocity: $$ \begin{equation} F_d^{(q)} = -{1\over2}C_D\varrho A|v|v, \end{equation} $$ where \( C_D \) is a dimensionless drag coefficient depending on the body's shape, and \( A \) is the cross-sectional area as produced by a cut plane, perpendicular to the motion, through the thickest part of the body. The superscripts \( \,{}^q \) and \( \,{}^S \) in \( F_d^{(S)} \) and \( F_d^{(q)} \) indicate Stokes drag and quadratic drag, respectively.
All the mentioned forces act in the vertical direction. Newton's second law of motion applied to the body says that the sum of these forces must equal the mass of the body times its acceleration \( a \) in the vertical direction. $$ \begin{equation*} ma = F_g + F_d^{(S)} + F_b \tp\end{equation*} $$ Here we have chosen to model the fluid resistance by the Stokes drag. Inserting the expressions for the forces yields $$ ma = -mg - 3\pi d\mu v + \varrho gV \tp $$ The unknowns here are \( v \) and \( a \), i.e., we have two unknowns but only one equation. From kinematics in physics we know that the acceleration is the time derivative of the velocity: \( a = dv/dt \). This is our second equation. We can easily eliminate \( a \) and get a single differential equation for \( v \): $$ m{dv\over dt} = -mg - 3\pi d\mu v + \varrho gV \tp $$ A small rewrite of this equation is handy: We express \( m \) as \( \varrho_bV \), where \( \varrho_b \) is the density of the body, and we divide by the mass to get $$ \begin{equation} v'(t) = - \frac{3\pi d\mu}{\varrho_b V} v + g\left(\frac{\varrho}{\varrho_b} -1\right) \tag{74} \tp \end{equation} $$ We may introduce the constants $$ \begin{equation} a = \frac{3\pi d\mu}{\varrho_b V},\quad b = g\left(\frac{\varrho}{\varrho_b} -1\right), \end{equation} $$ so that the structure of the differential equation becomes obvious: $$ \begin{equation} v'(t) = -av(t) + b \tag{75} \tp \end{equation} $$ The corresponding initial condition is \( v(0)=v_0 \) for some prescribed starting velocity \( v_0 \).
This derivation can be repeated with the quadratic drag force \( F_d^{(q)} \), leading to the result $$ \begin{equation} v'(t) = -{1\over2}C_D{\varrho A\over\varrho_b V}|v|v + g\left({\varrho\over\varrho_b} - 1\right) \tp \tag{76} \end{equation} $$ Defining $$ \begin{equation} a = {1\over2}C_D{\varrho A\over\varrho_b V}, \end{equation} $$ and \( b \) as above, we can write (76) as $$ \begin{equation} v'(t) = -a|v|v + b \tp \tag{77} \end{equation} $$
An interesting aspect of (75) and (77) is whether \( v \) will approach a final constant value, the so-called terminal velocity \( v_T \), as \( t\rightarrow\infty \). A constant \( v \) means that \( v'(t)\rightarrow 0 \) as \( t\rightarrow\infty \) and therefore the terminal velocity \( v_T \) solves $$0 = -av_T + b$$ and $$ 0 = -a|v_T|v_T + b \tp $$ The former equation implies \( v_T = b/a \), while the latter has solutions \( v_T =-\sqrt{|b|/a} \) for a falling body (\( v_T < 0 \)) and \( v_T = \sqrt{b/a} \) for a rising body (\( v_T>0 \)).
Both governing equations, the Stokes' drag model (75) and the quadratic drag model (77), can be readily solved by the Forward Euler scheme. For higher accuracy one can use the Crank-Nicolson method, but a straightforward application this method results a nonlinear equation in the new unknown value \( v^{n+1} \) when applied to (77): $$ \begin{equation} \frac{v^{n+1}-v^n}{\Delta t} = -a\half(|v^{n+1}|v^{n+1} + |v^n|v^n) + b \tag{78} \tp \end{equation} $$ However, instead of approximating the term \( -|v|v \) by an arithmetic average, we can use a geometric mean: $$ \begin{equation} (|v|v)^{n+\half} \approx |v^n|v^{n+1} \tp \end{equation} $$ The error is of second order in \( \Delta t \), just as for the arithmetic average and the centered finite difference approximation in (78). With this approximation trick, the discrete equation $$ \frac{v^{n+1}-v^n}{\Delta t} = - a|v^{n}|v^{n+1} + b $$ becomes a linear equation in \( v^{n+1} \), and we can therefore easily solve for \( v^{n+1} \): $$ \begin{equation} v^{n+1} = \frac{v_n + \Delta t b^{n+\half}}{1 + \Delta t a^{n+\half}|v^{n}|}\tp \tag{79} \end{equation} $$
Suitable values of \( \mu \) are \( 1.8\cdot 10^{-5}\hbox{ Pa}\, \hbox{s} \) for air and \( 8.9\cdot 10^{-4}\hbox{ Pa}\, \hbox{s} \) for water. Densities can be taken as \( 1.2\hbox{ kg/m}^3 \) for air and as \( 1.0\cdot 10^3\hbox{ kg/m}^3 \) for water. For considerable vertical displacement in the atmosphere one should take into account that the density of air varies with the altitude, see the section Decay of atmospheric pressure with altitude. One possible density variation arises from the one-layer model in the mentioned section.
Any density variation makes \( b \) time dependent and we need \( b^{n+\half} \) in (79). To compute the density that enters \( b^{n+\half} \) we must also compute the vertical position \( z(t) \) of the body. Since \( v=dz/dt \), we can use a centered difference approximation: $$ \frac{z^{n+\half} - z^{n-\half}}{\Delta t} = v^n \quad\Rightarrow\quad z^{n+\half} = z^{n-\half}+\Delta t\, v^n\tp$$ This \( z^{n+\half} \) is used in the expression for \( b \) to compute \( \varrho(z^{n+\half}) \) and then \( b^{n+\half} \).
The drag coefficient \( C_D \) depends heavily on the shape of the body. Some values are: 0.45 for a sphere, 0.42 for a semi-sphere, 1.05 for a cube, 0.82 for a long cylinder (when the center axis is in the vertical direction), 0.75 for a rocket, 1.0-1.3 for a man in upright position, 1.3 for a flat plate perpendicular to the flow, and 0.04 for a streamlined, droplet-like body.
To verify the program, one may assume a heavy body in air such that the \( F_b \) force can be neglected, and further assume a small velocity such that the air resistance \( F_d \) can also be neglected. This can be obtained by setting \( \mu \) and \( \varrho \) to zero. The motion then leads to the velocity \( v(t)=v_0 - gt \), which is linear in \( t \) and therefore should be reproduced to machine precision (say tolerance \( 10^{-15} \)) by any implementation based on the Crank-Nicolson or Forward Euler schemes.
Another verification, but not as powerful as the one above, can be based on computing the terminal velocity and comparing with the exact expressions. The advantage of this verification is that we can also the test situation \( \varrho\neq 0 \).
As always, the method of manufactured solutions can be applied to test the implementation of all terms in the governing equation, but the solution then has no physical relevance in general.
Applying scaling, as described in the section Scaling, will for the linear case reduce the need to estimate values for seven parameters down to choosing one value of a single dimensionless parameter $$ \beta = \frac{\varrho_b gV\left(\frac{\varrho}{\varrho_b} -1\right)}{3\pi d\mu I},$$ provided \( I\neq 0 \). If the motion starts from rest, \( I=0 \), the scaled problem \( \bar u'=1-\bar u \), \( \bar u(0)=0 \), has no need for estimating physical parameters. This means that there is a single universal solution to the problem of a falling body starting from rest: \( \bar u(t) = 1 - e^{-\bar t} \). All real physical cases correspond to stretching the \( \bar t \) axis and the \( \bar u \) axis in this dimensionless solution. More precisely, the physical velocity \( u(t) \) is related to the dimensionless velocity \( \bar u(\bar t) \) through $$ u = \frac{\varrho_bgV\left(\frac{\varrho}{\varrho_b} -1\right)}{3\pi d\mu}\bar u(t/(g(\varrho/\varrho_b -1)))\tp$$
Suppose we have a partial differential equation $$ \frac{\partial u}{\partial t} = \alpha\frac{\partial^2u}{\partial x^2} + f(x,t), $$ with boundary conditions \( u(0,t)=u(L,t)=0 \) and initial condition \( u(x,0)=I(x) \). One may express the solution as $$ u(x,t) = \sum_{k=1}^m A_k(t)e^{ikx\pi/L},$$ for appropriate unknown functions \( A_k \), \( k=1,\ldots,m \). We use the complex exponential \( e^{ikx\pi/L} \) for easy algebra, but the physical \( u \) is taken as the real part of any complex expression. Note that the expansion in terms of \( e^{ikx\pi/L} \) is compatible with the boundary conditions: all functions \( e^{ikx\pi/L} \) vanish for \( x=0 \) and \( x=L \). Suppose we can express \( I(x) \) as $$ I(x) = \sum_{k=1}^m I_ke^{ikx\pi/L} \tp $$ Such an expansion can be computed by well-known Fourier expansion techniques, but the details are not important here. Also, suppose we can express the given \( f(x,t) \) as $$ f(x,t) = \sum_{k=1}^m b_k(t)e^{ikx\pi/L} \tp $$ Inserting the expansions for \( u \) and \( f \) in the differential equations demands that all terms corresponding to a given \( k \) must be equal. The calculations results in the follow system of ODEs: $$ A_k'(t) = -\alpha\frac{k^2\pi^2}{L^2} + b_k(t),\quad k=1,\ldots,m \tp $$ From the initial condition $$ u(x,0)=\sum_k A_k(0)e^{ikx\pi/L}=I(x)=\sum_k I_k e^{(ikx\pi/L)},$$ it follows that \( A_k(0)=I_k \), \( k=1,\ldots,m \). We then have \( m \) equations of the form \( A_k'=-a A_k +b \), \( A_k(0)=I_k \), for appropriate definitions of \( a \) and \( b \). These ODE problems independent each other such that we can solve one problem at a time. The outline technique is a quite common approach for solving partial differential equations.
Remark. Since \( a_k \) depends on \( k \) and the stability of the Forward Euler scheme demands \( a_k\Delta t \leq 1 \), we get that \( \Delta t \leq \alpha^{-1}L^2\pi^{-2} k^{-2} \). Usually, quite large \( k \) values are needed to accurately represent the given functions \( I \) and \( f \) and then \( \Delta t \) needs to be very small for these large values of \( k \). Therefore, the Crank-Nicolson and Backward Euler schemes, which allow larger \( \Delta t \) without any growth in the solutions, are more popular choices when creating time-stepping algorithms for partial differential equations of the type considered in this example.
Show in detail how we can apply the ideas of the Forward Euler,
Backward Euler, and Crank-Nicolson
discretizations to derive explicit
computational formulas for new temperature values in Newton's law of
cooling (see the section Newton's law of cooling):
$$
\begin{equation}
{dT\over dt} = -k(T-T_s(t)),\quad T(0)=T_0\tp
\tag{80}
\end{equation}
$$
Here, \( T \) is the temperature of the body, \( T_s(t) \) is the temperature
of the surroundings, \( t \) is time, \( k \) is the heat transfer
coefficient, and \( T_0 \) is the initial temperature of the body.
Summarize the discretizations in a \( \theta \)-rule
such that you can get the three
schemes from a single formula by varying the \( \theta \) parameter.
Filename: schemes_cooling.pdf.
The goal of this exercise is to implement the schemes from Exercise 10: Derive schemes for Newton's law of cooling and investigate several approaches for verifying the implementation.
a)
Implement the \( \theta \)-rule from
Exercise 10: Derive schemes for Newton's law of cooling
cooling(T0, k, T_s, t_end, dt, theta=0.5),
where T0 is the initial temperature, k is
the heat transfer coefficient, T_s is a function of t for
the temperature of the
surroundings, t_end is the end time of the simulation, dt is the
time step, and theta corresponds to \( \theta \). The cooling
function should return the temperature as an array T of values at
the mesh points and the time mesh t.
b) In the case \( \lim_{t\rightarrow\infty}T_s(t)=C=\mbox{const} \), explain why \( T(t)\rightarrow C \). Construct an example where you can illustrate this property in a plot. Implement a corresponding test function that checks the correctness of the asymptotic value of the solution.
c) A piecewise constant surrounding temperature, $$ T_s(t) = \left\lbrace\begin{array}{ll} C_0,& 0\leq t\leq t^*\\ C_1, & t>t^*,\end{array}\right. $$ corresponds to a sudden change in the environment at \( t=t^* \). Choose \( C_1=2T_0 \), \( C_2=\frac{1}{2}T_0 \), and \( t^*=4/k \). Plot the solution and explain why it seems physically reasonable.
d) Find a value of \( \Delta t \) such that you can illustrate stability problems of the Crank-Nicolson scheme in the test case from c).
e)
Find an expression for the exact solution of
(80). Construct a test case and compare the
numerical and exact solution in a plot. Find a value of the time step
\( \Delta t \) such that the two solution curves cannot (visually) be
distinguished from each other. Many scientists will claim that such a
plot provides evidence for a correct implementation, but point out why
there still may be errors in the code. Can you introduce bugs in the
cooling function and still achieve visually coinciding curves?
Hint. The exact solution can be derived by multiplying (80) by the integrating factor \( e^{kt} \).
f)
Implement a test function for checking that the solution returned by
the cooling function is identical to the exact numerical
solution of the problem (to machine precision) when \( T_s \) is constant.
Hint. The exact solution of the discrete equations in the case \( T_s \) is a constant can be found by introducing \( u=T-T_s \) to get a problem \( u'=-ku \), \( u(0)=T_0-T_s \). The solution of the discrete equations is then of the form \( u^{n}=(T_0-T_s)A^n \) for some amplification factor \( A \). This approach can be used to derive an expression for \( T^n \).
Filename: cooling.py.
A detective measures the temperature of a dead body to be 26.7 C at 2 pm. One hour later the temperature is 25.8 C. The question is when death occurred.
Assume that Newton's law of cooling (80) is an
appropriate mathematical model for the evolution of the temperature in
the body. First, determine \( k \) in (80) by
formulating a Forward Euler approximation with one time steep from
time 2 am to time 3 am, where knowing the two temperatures allows for
finding \( k \). Assume the temperature in the air to be 20 C. Thereafter,
simulate the temperature evolution from the time of murder, taken as
\( t=0 \), when \( T=37\hbox{ C} \), until the temperature reaches 25.8 C. The
corresponding time allows for answering when death occurred.
Filename: detective.py.
The surrounding temperature \( T_s \) in Newton's law of cooling
(80) may vary in time. Assume that the
variations are periodic with period \( P \) and amplitude \( a \) around
a constant mean temperature \( T_m \):
$$
\begin{equation}
T_s(t) = T_m + a\sin\left(\frac{2\pi}{P}t\right)
\tp
\end{equation}
$$
Simulate a process with the following data: \( k=20 \hbox{ min}^{-1} \),
\( T(0)=5 \) C, \( T_m=25 \) C, \( a=2.5 \) C, and \( P=1 \) h. Also experiment with
\( P=10 \) min and \( P=3 \) h. Plot \( T \) and \( T_s \) in the same plot.
Filename: osc_cooling.py.
The Carbon-14 isotope, whose radioactive decay is used extensively in dating organic material that is tens of thousands of years old, has a half-life of \( 5,730 \) years. Determine the age of an organic material that contains 8.4 percent of its initial amount of Carbon-14. Use a time unit of 1 year in the computations. The uncertainty in the half time of Carbon-14 is \( \pm 40 \) years. What is the corresponding uncertainty in the estimate of the age?
Hint. Use simulations with \( 5,730\pm 40 \) y as input and find the corresponding interval for the result.
Filename: carbon14.py.
The purpose of this exercise is to implement the stochastic model described in the section Radioactive Decay and show that its mean behavior approximates the solution of the corresponding ODE model.
The simulation goes on for a time interval \( [0,T] \) divided into \( N_t \) intervals of length \( \Delta t \). We start with \( N_0 \) atoms. In some time interval, we have \( N \) atoms that have survived. Simulate \( N \) Bernoulli trials with probability \( \lambda\Delta t \) in this interval by drawing \( N \) random numbers, each being 0 (survival) or 1 (decay), where the probability of getting 1 is \( \lambda\Delta t \). We are interested in the number of decays, \( d \), and the number of survived atoms in the next interval is then \( N-d \). The Bernoulli trials are simulated by drawing \( N \) uniformly distributed real numbers on \( [0,1] \) and saying that 1 corresponds to a value less than \( \lambda\Delta t \):
# Given lambda_, dt, N
import numpy as np
uniform = np.random.uniform(N)
Bernoulli_trials = np.asarray(uniform < lambda_*dt, dtype=np.int)
d = Bernoulli_trials.size
Observe that uniform < lambda_*dt is a boolean array whose true
and false values become 1 and 0, respectively, when converted to an
integer array.
Repeat the simulation over \( [0,T] \) a large number of times, compute the average
value of \( N \) in each interval, and compare with the solution of
the corresponding ODE model.
Filename: stochastic_decay.py.
Consider two radioactive substances A and B. The nuclei in substance A decay to form nuclei of type B with a half-life \( A_{1/2} \), while substance B decay to form type A nuclei with a half-life \( B_{1/2} \). Letting \( u_A \) and \( u_B \) be the fractions of the initial amount of material in substance A and B, respectively, the following system of ODEs governs the evolution of \( u_A(t) \) and \( u_B(t) \): $$ \begin{align} \frac{1}{\ln 2} u_A' &= u_B/B_{1/2} - u_A/A_{1/2},\\ \frac{1}{\ln 2} u_B' &= u_A/A_{1/2} - u_B/B_{1/2}, \end{align} $$ with \( u_A(0)=u_B(0)=1 \).
Make a simulation program that solves for \( u_A(t) \) and \( u_B(t) \). Verify the implementation by computing analytically the limiting values of \( u_A \) and \( u_B \) as \( t\rightarrow \infty \) (assume \( u_A',u_B'\rightarrow 0 \)) and comparing these with those obtained numerically.
Run the program for the case of \( A_{1/2}=10 \) minutes and \( B_{1/2}=50 \) minutes.
Use a time unit of 1 minute. Plot \( u_A \) and \( u_B \) versus time in the same
plot.
Filename: radioactive_decay_2subst.py.
We consider the models for atmospheric pressure in the section Decay of atmospheric pressure with altitude. Make a program with three functions,
atmospheric_pressure.py.
Implement the Stokes' drag model (74) and the quadratic drag model (76) from the section Vertical motion of a body in a viscous fluid, using the Crank-Nicolson scheme and a geometric mean for \( |v|v \) as explained, and assume constant fluid density. At each time level, compute the Reynolds number Re and choose the Stokes' drag model if \( \hbox{Re} < 1 \) and the quadratic drag model otherwise.
The computation of the numerical solution should take place either in a stand-alone function (as in the section Making a solver function) or in a solver class that looks up a problem class for physical data (as in the section ref{decay:prog:se:class}). Create a module (see the section ref{decay:prog:se:module}) and equip it with nose tests (see the section ref{decay:prog:se:nose}) for automatically verifying the code.
Verification tests can be based on
A solution that is linear in \( t \) will also be an exact solution of the discrete equations in many problems. Show that this is true for linear drag (by adding a source term that depends on \( t \)), but not for quadratic drag because of the geometric mean approximation. Use the method of manufactured solutions to add a source term in the discrete equations for quadratic drag such that a linear function of \( t \) is a solution. Add a nose test for checking that the linear function is reproduced to machine precision in the case of both linear and quadratic drag.
Apply the software to a case where a ball rises in water. The
buoyancy force is here the driving force, but the drag will be
significant and balance the other forces after a short time. A soccer
ball has radius 11 cm and mass 0.43 kg. Start the motion from rest, set
the density of water, \( \varrho \), to \( 1000\hbox{ kg/m}^3 \), set the
dynamic viscosity, \( \mu \), to \( 10^{-3}\hbox{ Pa s} \), and use a drag
coefficient for a sphere: 0.45. Plot the velocity of the rising ball.
Filename: vertical_motion.py.
The aim of this project is to develop a general solver for the vertical motion of a body with quadratic air drag, verify the solver, apply the solver to a skydiver in free fall, and finally apply the solver to a complete parachute jump.
All the pieces of software implemented in this project should be realized as Python functions and/or classes and collected in one module.
a) Set up the differential equation problem that governs the velocity of the motion. The parachute jumper is subject to the gravity force and a quadratic drag force. Assume constant density. Add an extra source term be used for program verification. Identify the input data to the problem.
b) Make a Python module for computing the velocity of the motion. Also equip the module with functionality for plotting the velocity.
Hint 1. Use the Crank-Nicolson scheme with a geometric mean of \( |v|v \) in time to linearize the equation of motion with quadratic drag.
Hint 2.
You can either use functions or classes for implementation.
If you choose functions, make a function
solver that takes all the input data in the problem as
arguments and that returns the velocity (as a mesh function) and
the time mesh. In case of a class-based implementation, introduce
a problem class with the physical data
and a solver class with the numerical data and a solve method
that stores the velocity and the mesh in the class.
Allow for a time-dependent area and drag coefficient in the formula for the drag force.
c) Show that a linear function of \( t \) does not fulfill the discrete equations because of the geometric mean approximation used for the quadratic drag term. Fit a source term, as in the method of manufactured solutions, such that a linear function of \( t \) is a solution of the discrete equations. Make a nose test to check that this solution is reproduced to machine precision.
d) The expected error in this problem goes like \( \Delta t^2 \) because we use a centered finite difference approximation with error \( \Oof{\Delta t^2} \) and a geometric mean approximation with error \( \Oof{\Delta t^2} \). Use the method of manufactured solutions combined with computing convergence rate to verify the code. Make a nose test for checking that the convergence rate is correct.
e) Compute the drag force, the gravity force, and the buoyancy force as a function of time. Create a plot with these three forces.
Hint.
You can either make a function forces(v, t, plot=None)
that returns the forces (as mesh functions) and t, and shows
a plot on the screen and also saves the plot to a file with name
stored in plot
if plot is not None, or you can extend the solver class with
computation of forces and include plotting of forces in the
visualization class.
f) Compute the velocity of a skydiver in free fall before the parachute opens.
Hint. Meade and Struthers [7] provide some data relevant to skydiving. The mass of the human body and equipment can be set to \( 100 \) kg. A skydiver in spread-eagle formation has a cross-section of 0.5 \( \hbox{m}^2 \) in the horizontal plane. The density of air decreases varies altitude, but can be taken as constant, 1 \( \hbox{kg/m}^3 \), for altitudes relevant to skydiving (0-4000 m). The drag coefficient for a man in upright position can be set to 1.2. Start with a zero velocity. A free fall typically has a terminating velocity of 45 m/s. (This value can be used to tune other parameters.)
g) The next task is to simulate a parachute jumper during free fall and after the parachute opens. At time \( t_p \), the parachute opens and the drag coefficient and the cross-sectional area change dramatically. Use the program to simulate a jump from \( z=3000 \) m to the ground \( z=0 \). What is the maximum acceleration, measured in units of \( g \), experienced by the jumper?
Hint. Following Meade and Struthers [7], one can set the cross-section area perpendicular to the motion to 44 \( \hbox{m}^2 \) when the parachute is open. Assume that it takes 8 s to increase the area linearly from the original to the final value. The drag coefficient for an open parachute can be taken as 1.8, but tuned using the known value of the typical terminating velocity reached before landing: 5.3 m/s. One can take the drag coefficient as a piecewise constant function with an abrupt change at \( t_p \). The parachute is typically released after \( t_p=60 \) s, but larger values of \( t_p \) can be used to make plots more illustrative.
Filename: skydiving.py.
Vertical motion of a body in the atmosphere needs to take into account a varying air density if the range of altitudes is many kilometers. In this case, \( \varrho \) varies with the altitude \( z \). The equation of motion for the body is given in the section Vertical motion of a body in a viscous fluid. Let us assume quadratic drag force (otherwise the body has to be very, very small). A differential equation problem for the air density, based on the information for the one-layer atmospheric model in the section Decay of atmospheric pressure with altitude, can be set up as $$ \begin{align} p'(z) &= -\frac{Mg}{R^*(T_0+Lz)} p,\\ \varrho &= p \frac{M}{R^*T} \tp \end{align} $$ To evaluate \( p(z) \) we need the altitude \( z \). From the principle that the velocity is the derivative of the position we have that $$ \begin{equation} z'(t) = v(t), \end{equation} $$ where \( v \) is the velocity of the body.
Explain in detail how the governing equations can be discretized
by the Forward Euler and the Crank-Nicolson methods.
Filename: falling_in_variable_density.pdf.
Implement the Forward Euler or the Crank-Nicolson scheme derived in Exercise 20: Formulate vertical motion in the atmosphere. Demonstrate the effect of air density variation on a falling human, e.g., the famous fall of Felix Baumgartner. The drag coefficient can be set to 1.2.
Remark.
In the Crank-Nicolson scheme one must solve a \( 3\times 3 \) system of
equations at each time level, since \( p \), \( \varrho \), and \( v \) are
coupled, while each equation can be stepped forward at a time with the
Forward Euler scheme.
Filename: falling_in_variable_density.py.
Consider the ODE problem
$$
y'(x) = \left\lbrace\begin{array}{ll}
-1, & x < 0,\\
1, & x \geq 0
\end{array}\right.\quad x\in (-1, 1],
\quad y(1-)=1,
$$
which has the solution \( y(x)=|x| \).
Using a mesh \( x_0=-1 \), \( x_1=0 \), and \( x_2=1 \), calculate by hand
\( y_1 \) and \( y_2 \) from the Forward Euler, Backward Euler, Crank-Nicolson,
and Leapfrog methods. Use all of the former three methods for computing
the \( y_1 \) value to be used in the Leapfrog calculation of \( y_2 \).
Thereafter, visualize how these schemes perform for a uniformly partitioned
mesh with \( N=10 \) and \( N=11 \) points.
Filename: signum.py.
The goal of this exercise is to compute the value of a fortune subject to inflation and a random interest rate. Suppose that the inflation is constant at \( i \) percent per year and that the annual interest rate, \( p \), changes randomly at each time step, starting at some value \( p_0 \) at \( t=0 \). The random change is from a value \( p^n \) at \( t=t_n \) to \( p_n +\Delta p \) with probability 0.25 and \( p_n -\Delta p \) with probability 0.25. No change occurs with probability 0.5. There is also no change if \( p^{n+1} \) exceeds 15 or becomes below 1. Use a time step of one month, \( p_0=i \), initial fortune scaled to 1, and simulate 1000 scenarios of length 20 years. Compute the mean evolution of one unit of money and the corresponding standard deviation. Plot the mean curve along with the mean plus one standard deviation and the mean minus one standard deviation. This will illustrate the uncertainty in the mean curve.
Hint 1. The following code snippet computes \( p^{n+1} \):
import random
def new_interest_rate(p_n, dp=0.5):
r = random.random() # uniformly distr. random number in [0,1)
if 0 <= r < 0.25:
p_np1 = p_n + dp
elif 0.25 <= r < 0.5:
p_np1 = p_n - dp
else:
p_np1 = p_n
return (p_np1 if 1 <= p_np1 <= 15 else p_n)
Hint 2.
If \( u_i(t) \) is the value of the fortune in experiment number \( i \),
\( i=0,\ldots,N-1 \),
the mean evolution of the fortune is
$$ \bar u(t)= \frac{1}{N}\sum_{i=0}^{N-1} u_i(t),
$$
and the standard deviation is
$$ s(t) = \sqrt{\frac{1}{N-1}\left(- (\bar u(t))^2 +
\sum_{i=0}^{N-1} (u_i(t))^2\right)}
\tp
$$
Suppose \( u_i(t) \) is stored in an array u.
The mean and the standard deviation of the fortune
is most efficiently computed by
using two accumulation arrays, sum_u and sum_u2, and
performing sum_u += u and sum_u2 += u**2 after every experiment.
This technique avoids storing all the \( u_i(t) \) time series for
computing the statistics.
Filename: random_interest.py.
We shall study a population modeled by (65) where the environment, represented by \( r \) and \( f \), undergoes changes with time.
a)
Assume that there is a sudden drop (increase) in the birth (death)
rate at time \( t=t_r \),
because of limited nutrition or food supply:
$$ a(t) =\left\lbrace\begin{array}{ll}
r_0, & t < t_r,\\
r_0 - A, & t\geq t_r,\end{array}\right.
$$
This drop in population growth is compensated by a sudden net immigration
at time \( t_f>t_r \):
$$ f(t) =\left\lbrace\begin{array}{ll}
0, & t < t_f,\\
f_0, & t\geq t_a,\end{array}\right.
$$
Start with \( r_0 \) and make \( A>r_0 \). Experiment with
these and other parameters to
illustrate the interplay of growth and decay in such a problem.
Filename: population_drop.py.
b)
Now we assume that the environmental conditions changes periodically with
time so that we may take
$$ r(t) = r_0 + A\sin\left(\frac{2\pi}{P}t\right)
\tp
$$
That is, the combined birth and death rate oscillates around \( r_0 \) with
a maximum change of \( \pm A \) repeating over a period of length \( P \) in time.
Set \( f=0 \) and experiment with the other parameters to illustrate typical
features of the solution.
Filename: population_osc.py.
Solve the logistic ODE
(66) using a Crank-Nicolson scheme where
\( (u^{n+\half})^2 \) is approximated by a geometric mean:
$$ (u^{n+\half})^2 \approx u^{n+1}u^n
\tp
$$
This trick makes the discrete equation linear in \( u^{n+1} \).
Filename: logistic_CN.py.
The ODE model (68) was derived under the assumption
that \( r \) was constant. Perform an alternative derivation without
this assumption: 1) start with (67);
2) introduce a time step \( \Delta t \) instead of \( m \): \( \Delta t = 1/m \) if
\( t \) is measured in years; 3) divide by \( \Delta t \) and take the
limit \( \Delta t\rightarrow 0 \). Simulate a case where the inflation is
at a constant level \( I \) percent per year and the interest rate oscillates:
\( r=-I/2 + r_0\sin(2\pi t) \).
Compare solutions for \( r_0=I, 3I/2, 2I \).
Filename: interest_modeling.py.