This section is presents some applications of the first order ordinary differential equations.

Return to computing page for the first course APMA0330
Return to computing page for the second course APMA0340
Return to Mathematica tutorial for the first course APMA0330
Return to Mathematica tutorial for the second course APMA0340
Return to the main page for the course APMA0330
Return to the main page for the course APMA0340
Return to Part II of the course APMA0330



We present some applications of differential equations of order one.
Example 1: A 500 liter container initially contains 10 kg of salt. A brine mixture of 100 grams of salt per liter is entering the container at 6 liter per minute. The well-mixed contents
are being discharged from the tank at the rate of 6 liters per minute.

Express the amount of salt in the container as a function of time.

Salt is coming at the rate: 6*(0.1)=0.6 kg/min
d yin /dt =0.6 ; d yout/dt = 6x/500

dx/dt = 0.6 -6x/500 ; x(0)=10

DSolve[{x'[t]==6/10-6*x[t]/500, x[0]==10},x[t],t]
Out[15]= {{x[t] -> 10 E^(-3 t/250) (-4 + 5 E^(3 t/250))}}

Suppose that the rate of discharge is reduced to 5 liters per minute.

DSolve[{x'[t]==6/10-5*x[t]/(500+t), x[0]==10},x[t],t]
SolRule[t_] = Apart[x[t] /. First[%]]
Out[1]= {{x[t] -> (1/( 10 (500 + t)^5))(3125000000000000 + 187500000000000 t + 937500000000 t^2 + 2500000000 t^3 + 3750000 t^4 + 3000 t^5 + t^6)}}
Out[2]= 50 + t/10 - 1250000000000000/(500 + t)^5
To check the answer:
Simplify[SolRule'[t] == 6/10 - 5*SolRule[t]/(500 + t)]
Together[SolRule'[t] == Together[6/10 - 5*SolRule[t]/(500 + t)]]
Out[3]= True
Out[4]= True
Out[5]= 10

Example 2: The Solow–Swan model is an economic model of long-run economic growth set within the framework of neoclassical economics. It attempts to explain long-run economic growth by looking at capital accumulation, labor or population growth, and increases in productivity, commonly referred to as technological progress. At its core is a neoclassical (aggregate) production function, often specified to be of Cobb–Douglas type, which enables the model "to make contact with microeconomics". The model was developed independently by Robert Solow and Trevor Swan in 1956, and superseded the Keynesian Harrod–Domar model.

The neo-classical model was an extension to the 1946 Harrod–Domar model that included a new term: productivity growth. Important contributions to the model came from the work done by Solow and by Swan in 1956, who independently developed relatively simple growth models. The world is too complex to describe it inevery detail. We always need a simple but effective model in predicting how the world works. Solow's model fitted available data on US economic growth with some success. In 1987 Solow was awarded the Nobel Prize in Economics for his work. Today, economists use Solow's sources-of-growth accounting to estimate the separate effects on economic growth of technological change, capital, and labor.

The Solow model is based on the following assumptions:

  • Time variable is continuous.
  • Single good produced with a constant technology.
  • No government or international trade.
  • All factors of production are fully employed.
  • Labor force grows at constant rate n.
  • nitial values for capital, K0 and labor, L0 are given.
The neoclassical (Cobb--Douglas) aggregate production function function is assumed to be
\[ Y(t) = F \left( K(t), L(t) \right) = K(t)^{\alpha} \left( A(t)\,L(t) \right)^{1-\alpha} , \]
where t denotes time, 0 < α < 1 is the elasticity of output with respect to capital, and Y( t ) represents total production. A refers to labor-augmenting technology or “knowledge”, thus A L represents effective labor. All factors of production are fully employed, and initial values A( 0 ), K( 0 ), and L( 0 ) are given. The number of workers, i.e. labor, as well as the level of technology grow exogenously at rates n n and g, respectively:
  1. Solow, Robert M. (February 1956). "A contribution to the theory of economic growth". Quarterly Journal of Economics. 70 (1): 65–94. doi:10.2307/1884513.
  2. Solow, Robert M. (1957). "Technical change and the aggregate production function". Review of Economics and Statistics. 39 (3): 312–320. doi:10.2307/1926047.
  3. Swan, Trevor W. (November 1956). "Economic growth and capital accumulation". Economic Record. 32 (2): 334–361. doi:10.1111/j.1475-4932.1956.tb00434.x.


Vertical Motion

It is known from elementary physics that, in the absence of air friction, a ball thrown up from the ground into earth's atmosphere with initial speed v0 would attain a maximum altitude of \( v_0^2 /(2g) .\) In this case the return time is \( 2\, v_0 /g, \) independent of the ball's mass. Here g is the acceleration due to gravity. If the ball is thrown up from altitude y0 (which we later assume to be zero), then the time T0 spent traveling is given by

\[ \left\{ \begin{array}{c} \mbox{travel time with no air resistance} \\ \mbox{when thrown from height $y_0$} \end{array} \right\} \ = \ T_0 \ =\ \frac{v_0 + \sqrt{v_0^2 +2\, y_0\, g}}{g} \]

The presence of air influences the ball's motion: it experiences two forces acting on it---the force of gravity and the air resistance force. Let's define the symbol T as follows:

\[ \left\{ \begin{array}{c} \mbox{travel time {\it with} air resistance} \\ \mbox{when thrown from height $y_0$} \end{array} \right\} \ = \ T \]

Without air resistance, the object travels farther up than with air resistance. On the way down, without air resistance the object travels a larger distance, but it also gathers more speed. A natural question is, which travel time (with air resistance vs. no air resistance) is larger? Also, it is of interest to find the maximum altitude \( y_{\max} \) of the ball, the time Tmax to reach maximum altitude, and the time \( T_{\text{down}} \) to return back from ymax. Therefore, \( T_{\max} + T_{\text{down}} = T_{\text{total}} \) (the total time the ball spent in the air). The landing velocity is denoted by \( v_{\ell} . \)

Example 3: Air resistance is the force that acts in the direction opposite to the motion of an object through air. Air resistance depends on the shape, material, and orientation of the object, the density of the air, and the object's relative speed.

We would like to think that there is a nice formula for the air resistance in terms of speed and other variables. Such a formula would help in making calculations and predicting various quantities. A starting point for obtaining such a formula is our everyday experience. Based on our experience, a reasonable assumption to make\footnote{Our intuition based on everyday experience is limited to a small range of conditions. This may lead to erroneous assumptions.

It has been observed that, under suitable conditions, the magnitude of the air resistance is proportional to a power of the speed s=|v|:

\begin{equation} \tag{M} F \propto s^p, \qquad \mbox{which may be written as}\qquad F(v) = k \, |v|^p. \end{equation}
Here v is velocity, and both k and p are positive constants. For very small objects, such as a speck of dust (about 1 micrometer or 0.001mm), p=1 seems to give a reasonable formula for the air resistance. For larger, human scale objects moving at relatively large speed, p=2 works better. Therefore, the magnitude of the air resistance F as a function of velocity v is assumed be given by formula (M).

The air resistance force depends on the velocity (v) of the object at time t, so let us denote this force with the symbol F(v). Note that the air resistance, force F(v), always acts in the direction opposite to the motion. Therefore, F(v) acts in the down (negative) direction when the ball is moving up, and it acts in the up (positive) direction when the ball is moving down. If we measure the displacement y = y(t) vertically upwards from the ground, then \( v= {\text d}y/{\text d}t = \dot{y} \) is the velocity of the object. Newton's law of motion for the ball on the way up gives the differential equation

\begin{equation} \label{vert.1} m\, \dot{v} = -m\,g -F(v) \ , \qquad\mbox{or}\qquad \frac{{\text d}v}{{\text d}t} = - g - \frac{1}{m}\, F(v), \end{equation}
and on the way down,
\begin{equation} \label{vert.2} m\, \dot{v} = -m\,g +F(v)\ , \qquad\mbox{or}\qquad \frac{{\text d}v}{{\text d}t} = - g + \frac{1}{m}\, F(v). \end{equation}
Since we assume \( F(v) = k \, |v|^p , \) the equation of motion on way up becomes
\begin{equation} \label{vert.3} m\, \dot{v} = -m\,g -k\,|v|^p \qquad\mbox{or}\qquad \frac{{\text d}v}{{\text d}t} = -g-\frac{k}{m}\,|v|^p , \end{equation}
and on the way down,
\begin{equation} \label{vert.4} m\, \dot{v} = -m\, g +k\, |v|^p \qquad\mbox{or}\qquad \frac{{\text d}v}{{\text d}t} = -g+\frac{k}{m}\,|v|^p . \end{equation}

To find an equation for $v_{\ell}$, the landing velocity, we rewrite \refeq{vert.3} as \( {\text d}\,t = - {\text d}v/(g +F(v)/m) \) and integrate both sides from t=0 and v = v0 to \( t=T_{\max} \) and v=0. Here \( T_{\max} \) is the time to reach the maximum altitude \( y_{max} , \) which is also the time to have velocity v=0. We obtain,

\begin{equation} \label{Tmax.1} \int_{0}^{T_{max}}\,{\text d}t = - \int_{v_0}^{0} \frac{{\text d}\,v}{g+\frac{k}{m}\,|v|^p} \, , \quad \mbox{that is,} \quad T_{max} = - \int_{v_0}^{0} \frac{{\text d}\,v}{g+\frac{k}{m}\,|v|^p} \end{equation}
A similar formula is valid for time T. From \refeq{vert.4}, we get
\begin{equation} \label{Tmax.2} \int_{T_{max}}^T \, {\text d}t = \int_0^{v_\ell} \frac{ \,{\text d}v}{-g+\frac{k}{m}\,|v|^p} \quad\Longrightarrow \quad T - T_{\max} = \int_0^{v_\ell} \frac{ \,{\text d}v}{-g+\frac{k}{m}\,|v|^p} \end{equation}
Equating the time \( T_{max} \) in the above equations, we have
\begin{equation} \label{T landing} T \ = \ - \int_{v_0}^{0} \frac{{\text d}\,v}{g+\frac{k}{m}\,|v|^p} + \int_0^{v_\ell} \frac{{\text d}\,v}{-g+\frac{k}{m}\,|v|^p} \end{equation}
This equation gives T as a function of \( v_\ell . \) The next step is to find an equation for \( v_\ell . \)

To find an equation for $v_\ell$, we rewrite \refeq{vert.3} as \( v\,{\text d}t = -v\, {\text d}v/(g +F(v)/m) \) and integrate both sides from t=0 and v = v0 to \( t=T_{\max} \) and v=0. Here \( T_{\max} \) is the time to reach the maximum altitude \( y_{\max} , \) which is also the time to have velocity v=0. Using the fact that the integral of the velocity is the displacement, we obtain,

\begin{equation} \label{ymax.1} \int_{0}^{T_{max}} v\,{\text d}t = - \int_{v_0}^{0} \frac{v\,{\text d}v}{g+\frac{k}{m}\,|v|^p} \, , \quad \mbox{that is,} \quad y_{max} - y_0 = - \int_{v_0}^{0} \frac{v\,{\text d}v}{g+\frac{k}{m}\,|v|^p} \end{equation}
A similar formula is valid for the distance traveled down. From \refeq{vert.4} we get
\begin{equation} \label{ymax.2} \int_{T_{max}}^T v\, {\text d}t = \int_0^{v_\ell} \frac{v\,{\text d}v}{-g+\frac{k}{m}\,|v|^p} \quad\Longrightarrow \quad - y_{\max} = \int_0^{v_\ell} \frac{v\,{\text d}v}{-g+\frac{k}{m}\,|v|^p} \end{equation}
Equating the maximum distance \( y_{\max} \) traveled in both ways, we get an equation involving the landing velocity \( v_{\ell} : \)
\begin{equation} \label{v landing} \int_0^{v_\ell} \frac{v\,{\text d}v}{-g+\frac{k}{m}\,|v|^p} \ = \ - y_0 + \int_{v_0}^{0} \frac{v\,{\text d}v}{g+\frac{k}{m}\,|v|^p} \end{equation}
This equation is an equation where the unknown is \( v_{\ell} , \) which does not appear explicitly solved for.

g := 9.806
Solve[Integrate[v/(-g + 0.01*v^2), {v, 0, vl}] == Integrate[v/(g + 0.01*v^2), {v, 50, 0}], vl]
Out[2]= {{vl -> -26.5393}, {vl -> 26.5393}}

We would like to determine the ratio: \[ \gamma = T/T_0 = Tg/ \left( v_0 + \sqrt{v_0^2 +2 x_0 g}\right), \] where T time in air with air resistance and T0 is the time in air without air resistance.

For a tennis ball thrown upward with the initial velocity v0 =10, it is possible to find x0 that γ > 1 when p=0.9. In general, it is unknown for what values of p< 1, x0, and v0 we can achieve γ > 1.

Example 4: We consider a model of falling object, say a tennis ball, to a flat surface that moves up and down periodically. Using vertical axis directed upward, we denote v(t) as the velocity of the ball and y(t) as its position/height at time t. It has been observed that, under suitable conditions, the magnitude of the air resistance is proportional to the power of speeds=|v|:
\begin{equation} \tag{M} F \propto s^p, \qquad \mbox{which may be written as}\qquad F(v) = k \, |v|^p. \end{equation}
Note that the air resistance force F(v) always acts in the direction opposite to the motion. Therefore, F(v) acts in the down (negative) direction when the ball is moving up, and it acts in the up (positive) direction when the ball is moving down.

Suppose that initially at t=0, the ball of mass m is dropped from the altitude \( y=h > 1 \) without initial velocity. At the same time, it is assumed that the floor starts moving according to the formula \( z= \sin \omega t. \) When elastic ball hits a hard flat surface, it bounces back with the same velocity. It is assumed that the collision is totally elastic, so the ball loses no kinetic energy in the collision, and its speed after collision is the same as before the collision. At this point, ignore the time needed for the ball to be deformed during collision before fully rebounded and has lifted off from the surface instantly. Hence the ball can be treated as a rigid body with negligible deformation during impact.

After collision, the ball climbs up until its velocity becomes zero, and then the ball falls vertically downward under the influence of gravity, hits the the moving floor, and bounces back.

Derivation of a differential equation

Newton's law of motion for the ball on the way down is

\begin{equation*} %\label{vert.1} m\, \dot{v} = -m\,g +F(v)\ , \qquad\mbox{or}\qquad \frac{{\text d}v}{{\text d}t} = - g + \frac{1}{m}\, F(v), \end{equation*}
and on the way up
\begin{equation*} %\label{vert.2} m\, \dot{v} = -m\,g -F(v) \ , \qquad\mbox{or}\qquad \frac{{\text d}v}{{\text d}t} = - g - \frac{1}{m}\, F(v), \end{equation*}
where g is the acceleration due to gravity. Since we assume \( F(v) = k \, |v|^p ,\) the equation of motion on the way down becomes
\begin{equation*} %\label{vert.3} m\, \dot{v} = -m\, g +k\, |v|^p \qquad\mbox{or}\qquad \frac{{\text d}v}{{\text d}t} = -g+\frac{k}{m}\,|v|^p . \end{equation*}
and on way up
\begin{equation} %\label{vert.4} m\, \dot{v} = -m\,g -k\,|v|^p \qquad\mbox{or}\qquad \frac{{\text d}v}{{\text d}t} = -g-\frac{k}{m}\,|v|^p , \end{equation}
Assuming that p=2 for a tennis ball, the above differential equations can be integrated using separation of variables:
\begin{equation} %\label{vert.1} \int \frac{{\text d}v}{-g+\frac{k}{m}\,|v|^2} = \int {\text d}t \qquad\Longrightarrow \qquad \sqrt{\frac{m}{g\,k}}\, \mbox{Arctanh} \left( v\,\sqrt{\frac{k}{g\,m}} \right) = t , \end{equation}
and on way up
\begin{equation} %\label{vert.2} \int \frac{{\text d}v}{g+\frac{k}{m}\,|v|^2} = -\int {\text d}t \qquad\Longrightarrow \qquad \sqrt{\frac{m}{g\,k}}\, \mbox{Arctan} \left( v\,\sqrt{\frac{k}{g\,m}} \right) = -t , \end{equation}
\[ \mbox{Arctanh} \, (x) = \frac{1}{2} \, \ln \frac{1+x}{1-x} \quad \mbox{for} \quad |x| < 1. \]
Since the velocity v(t) is the derivative of ball's position y(t), you may need to integrate
\[ \int \arctan (\alpha v)\, {\text d}v = v\,\arctan (\alpha v) - \frac{1}{2\alpha} \, \ln \left( 1 + \alpha^2 v^2 \right) , \quad \int \mbox{arctanh} (\alpha v)\, {\text d}v = v\,\mbox{arctanh} (\alpha v) - \frac{1}{2\alpha} \, \ln \left\vert 1 - \alpha^2 v^2 \right\vert . \]
For example, at the initial stage, you need to solve two IVPs:
\[ \frac{{\text d}v}{{\text d}t} = -g+\frac{k}{m}\,|v|^2, \quad v(0) =0; \qquad \frac{{\text d}y}{{\text d}t} = v, \quad y(0) =h . \]

Input parameters

 \( g \approx 9.806 \) m/sec2  the acceleration due to gravity near sea level at 45 deg. latitude;
 m=0.08  mass of the object, in kg because a tennis ball is about 80 grams;
 k=0.02   drag coefficient, positive;
 p=2   power of the speed term in the resistance force;
  ω = π   frequency of the oscillating floor;
 \( y_0 =h > 1 \)  initial altitude, positive, in meters.

Derivation of solution

A ball that is dropped from height h> 1 can be described by its velocity v(t) and position y(t):

\[ \frac{{\text d}v}{{\text d}t} = -g+\frac{k}{m}\,|v|^2, \quad v(0) =0; \qquad \frac{{\text d}y}{{\text d}t} = v, \quad y(0) =h . \]
Separation of variables yields
\[ \int_0^v \frac{{\text d}v}{-g+\frac{k}{m}\,|v|^2} = \int_0^t {\text d}t \qquad \Longrightarrow \qquad \mbox{Arctanh} \left( \sqrt{\frac{k}{gm}} \, v \right) = - \sqrt{\frac{gk}{m}} \,t . \]
Therefore, we find the velocity v(t) on first stage of falling from hight h> 1 explicitly:
\[ v(t) = - \sqrt{\frac{gm}{k}} \,\tanh \left( \sqrt{\frac{gk}{m}} \,t \right) , \qquad k\,v^2 < g\,m. \]
Then we find ball's position by integrating v(t):
\[ y(t) = h- \frac{m}{k}\,\ln \cosh \left( \sqrt{\frac{gk}{m}} \,t \right) . \]
This equation is valid until the ball meets the surface \( z= \sin \pi t. \) Therefore, we need to solve the transcendent equation:
\[ \sin \pi t = h- \frac{m}{k}\,\ln \cosh \left( \sqrt{\frac{gk}{m}} \,t \right) . \]

Mathematica confirms:
Assuming[v > 0 && k > 0 && m > 0 && g > 0, Integrate[1/(-g + k/m*v1^2), {v1, 0, v}]]
Integrate[A*Tanh[B*t], t]
k := 0.02; g := 9.806; m := 0.08; h0 := 2; (* h0 is the initial height *)
FindRoot[Sin[\[Pi] t] == h0 - 4*Log[Cosh[1.5657266683556232*t]], {t, 0.4}]
vv[t_] = omega*Cos[omega*t]
omega = Pi
T1 = 0.4716548296910227
y1[t_] = h0 - m/k*Log[Cosh[t*Sqrt[g*k/m]]]
V1 = kk*vv[T1] + Sqrt[g*m/k]*Tanh[T1*Sqrt[g*k/m]]
Y1 = Sin[omega*T1]
Since Mathematica provides the approximate value of T1 to be 0.4716548296910227 for h0 = 2, we denote it by T1. The velocity \( v(t) = \omega\,\cos (\omega \,t) \) of the floor/racket at T1 is \( vv(T1) = V1 \approx 0.279386 \) positive, which indicate that the racquet is moving up. So we add 80% of its velocity to the ball.

On the second stage, the ball bounced up with the initial velocity \( V1 = kk*vv(T1) + \sqrt{\frac{gm}{k}} \tanh \left( \sqrt{\frac{gk}{m}}\,T1 \right) \approx 4.158 \) and from the position \( Y1 = \sin \left( \omega\,T1 \right) \approx 0.996038. \) Here kk is the coefficients of elastic damping of the racket, which we set to be 0.8. Therefore, we need to solve two initial value problems:

\[ \frac{{\text d}v}{{\text d}t} = -g-\frac{k}{m}\,|v|^2 , \quad v(T1) = V1, \qquad \frac{{\text d}y}{{\text d}t} = v(t), \quad y(T1) = Y1, \quad \mbox{for }t\ge T1. \]
Separating variables and integrating, we obtain
\[ v(t) = V1 + \sqrt{\frac{gm}{k}} \, \tan \left[ \sqrt{\frac{gk}{m}} \left( T1 -t \right) \right] , \quad T1 \le t \le T2 , \]
where T2 is the value of time when v(T2) =0. Then we find ball's position:
\[ y(t) = Y1+ V1\left( t- T1 \right) +\frac{m}{k} \,\ln \cos \sqrt{\frac{gk}{m}} \left( T1 -t \right) , \quad T1 \le t \le T2 . \]
Mathematica find its value to be \( T2 \approx 0.846 : \)
v2[t_] := V1 + Sqrt[g*m/k]*Tan[Sqrt[g*k/m]*(T1 - t)]
T2 = t /. FindRoot[v2[t] == 0, {t, 0.8}]
The altitude of the ball with zero speed will be
y2[t_] = Y1 + V1*(t - T1) + m/k*Log[Cos[Sqrt[g*k/m]*(T1 - t)]]
v3[t_] = -Sqrt[g*m/k]*Tanh[Sqrt[g*k/m]*(t - T2)]
Y2 = y2[T2]
Then we find time (which we denote as T3) when the ball meets the racket:
y3[t_] = Y2 - m/k*Log[Cosh[Sqrt[g*k/m]*(t - T2)]]
T3 = t /. FindRoot[y3[t] == Sin[omega*t], {t, 2}]
The velocity of the racket at time contact is positive and we transfer 80% of its speed to the ball.
The position of the ball at t=T2 will be \( y(T2) \approx 1.822 . \) After T2, the ball starts falling down and we need to solve the initial value problems:
\[ \frac{{\text d}v}{{\text d}t} = -g+\frac{k}{m}\,|v|^2, \quad v(T2) =0; \qquad \frac{{\text d}y}{{\text d}t} = v(t), \quad y(T2) =Y2 \approx 1.822 . \]
We again separate variables and obtain
\[ v(t) = -\sqrt{\frac{gm}{k}} \, \tanh \left[ \sqrt{\frac{gk}{m}} \left( t- T2 \right) \right] , \quad T2 \le t\le T3 , \]
\[ y(t) = Y2 - \frac{m}{k} \, \ln \cosh \left[ \sqrt{\frac{gk}{m}} \left( t- T2 \right) \right] , \quad T2 \le t\le T3 . \]
To find T3, we need to solve the transcendent equation:
\[ \sin \left( \omega t \right) = Y2 - \frac{m}{k} \, \ln \cosh \left[ \sqrt{\frac{gk}{m}} \left( t- T2 \right) \right] . \]
We plot with Mathematica to be sure that the racket does not hit the ball on this stage:
y[t_] = Piecewise[{{h0 - m/k*Log[Cosh[t*Sqrt[g*k/m]]],
0 < t < T1}, {Y1 + V1*(t - T1) + m/k*Log[Cos[Sqrt[g*k/m]*(T1 - t)]],
T1 < t <= T2}, {Y2 - 4*Log[Cosh[Sqrt[g*k/m]*(t - T2)]], T2 < t <= T3}}]
Plot[{y[t], Sin[omega*t]}, {t, 0, T3}]
At \( t= T3 \approx 1.67 , \) the ball meets the racket and starts climbing up. So its velocity should be the solution of the following initial value problem:
\[ \frac{{\text d}v}{{\text d}t} = -g-\frac{k}{m}\,|v|^2 , \quad v(T3) = V3, \qquad \frac{{\text d}y}{{\text d}t} = v(t), \quad y(T3) = Y3, \quad \mbox{for } T3 \le t \le T4, \]
where T4 is time where the velocity is zero, and \( V3 = kk*\sin (\omega\,T3) + \sqrt{\frac{gm}{k}} \tanh \left( \sqrt{\frac{gk}{m}}\,T3 \right) \approx 4.69381 . \) Mathematica provides the numerical values:
V3 = kk*Sin[omega*T3] + Sqrt[g*m/k]*Tanh[Sqrt[g*k/m]*(T3 - T2)]
Then we define ball's velocity and its position:
Y3 = Sin[omega*T3]
v4[t_] = V3 + Sqrt[g*m/k]*Tan[(T3 - t)*Sqrt[g*k/m]]
y4[t_] = Y3 + V3*(t - T3) + m/k*Log[Cos[(T3 - t)*Sqrt[g*k/m]]]
T4 = t /. FindRoot[v4[t] == 0, {t, 2}]
Next, we check whether the racket hits the ball from below
TT4 = t /. FindRoot[y4[t] == Sin[omega*t], {t, 2}]
Since this number TT4 is less than previously found T4, we conclude that the racket hits the ball on its way up. At this moment TT4, we define the position and the velocity of the ball:
VV3 = vv[TT4] + v4[TT4]
YY4 = Sin[omega*TT4]
and then define the velocity and new position
vv4[t_] = VV3 + Sqrt[g*m/k]*Tan[(TT4 - t)*Sqrt[g*k/m]]
yy4[t_] = YY4 + VV3*(t - TT4) + m/k*Log[Cos[(TT4 - t)*Sqrt[g*k/m]]]
y[t_] = Piecewise[{{h0 - m/k*Log[Cosh[t*Sqrt[g*k/m]]],
0 < t < T1}, {Y1 + V1*(t - T1) + m/k*Log[Cos[Sqrt[g*k/m]*(T1 - t)]],
T1 < t <= T2}, {Y2 - 4*Log[Cosh[Sqrt[g*k/m]*(t - T2)]],
T2 < t <= T3}, {Y3 + V3*(t - T3) +
m/k*Log[Cos[(T3 - t)*Sqrt[g*k/m]]], T3 < t < TT4}, {yy4[t], TT4 < t < 2.081}}]
Plot[{y[t], Sin[omega*t]}, {t, 0, 2.081}, PlotStyle -> {Thick, Thick}]
The time interval where the ball oscillates over the racket is not enough to show on the picture because time intervals are too small.
Graphs of sin (π t) and position of y(t).
 Claude Perrault  Sir Isaac Newton  Christiaan Huygens

A tractrix (from the Latin verb trahere "pull, drag"; plural: tractrices) is the curve along which an object moves, under the influence of friction, when pulled on a horizontal plane by a line segment attached to a tractor (pulling) point that moves at a right angle to the initial line between the object and the puller at an infinitesimal speed. It is therefore a curve of pursuit. It was first introduced by Claude Perrault in 1670, and later studied by Isaac Newton (1676) and Christiaan Huygens (1692).

Claude Perrault (1613 -- 1688) was a French architect, best known for his participation in the design of the east façade of the Louvre in Paris. He also achieved success as a physician and anatomist, and as an author, who wrote treatises on physics and natural history. Perrault was born and died in Paris. Aside from his influential architecture, he became well known for his translation of the ten books of Vitruvius, the only surviving Roman work on architecture, into French, written at the instigation of Colbert, and published, with Perrault's annotations, in 1673. His treatise on the five classical orders of architecture followed in 1683. As physician and natural philosopher with a medical degree from the University of Paris, Perrault became one of the first members of the French Academy of Sciences when it was founded in 1666.

Sir Isaac Newton (1642 -- 1726/27) was an English mathematician, astronomer, theologian and physicist (described in his own day as a "natural philosopher") who is widely recognized as one of the most influential scientists of all time and a key figure in the scientific revolution. His book Philosophiae Naturalis Principia Mathematica ("Mathematical Principles of Natural Philosophy"), first published in 1687, laid the foundations of classical mechanics. Newton also made pathbreaking contributions to optics, and he shares credit with Gottfried Wilhelm Leibniz for developing the infinitesimal calculus.

Christiaan Huygens (Latin: Hugenius; 1629 -- 1695) was a prominent Dutch mathematician and scientist. He is known particularly as an astronomer, physicist, probabilist and horologist. Huygens was a leading scientist of his time. His work included early telescopic studies of the rings of Saturn and the discovery of its moon Titan, the invention of the pendulum clock and other investigations in timekeeping. He published major studies of mechanics and optics (having been one of the most influential proponents of the wave theory of light), and pioneered work on games of chance.

Example 5: In order to plot tractrix curve, use the following code:

tractrix[a_][t_] := a*{Sin[t], Cos[t] + Log[Tan[t/2]]};
ParametricPlot[tractrix[a][t] // Evaluate, {t, 0, .99*\[Pi]},
PlotRange -> {0, 7}], {a, 1, 6}]
Plot[y'[x] = -Sqrt[a^2 - x^2]/x, {x, 0, 20},
PlotRange -> {-10, 10}], {a, 0, 20}]
sol = DSolve[{y'[x] == -Sqrt[a^2 - x^2]/x, y[a] == 0}, y, x]
Plot[-Sqrt[a^2 - x^2] + a Log[a] - a Log[a^2] - a Log[x] +
a Log[a^2 + a Sqrt[a^2 - x^2]], {x, 0, 20}, PlotRange -> All], {a,
1, 20}]

Example 6:: If a person becomes susceptible to a desease after recovering from it (like gonorrhea, meningitis, and streptococcal sore throat), then the percent of persons susceptible to becoming infected with the desease, S(t), and the percent of people in the population infected with the desease, I(t), can be modeled by the system of differential equations

\[ \begin{cases} \dot{S} &= -\lambda\,I(t)\,S(t) + \gamma\, I(t) + \mu - \mu\,S(t) , \\ \dot{I} &= \lambda\,I(t)\,S(t) - \gamma\, I(t) - \mu\, I(t) , \end{cases} \qquad\mbox{subject to} \qquad \begin{cases} S(t) + I(t) =1, \\ S(0) = S_0, \quad I(0) = I_0 . \end{cases} \]
As usual, we use dot to represent the derivative with respect to time variable t. The above model is called an SIS model (susceptible-infected-susceptible model) with vital dynamics (where μ and γ represent the birth and death rates) because once an individual has recovered from the desease, the individual again becomes susceptible to the desease. The λSI term is understood as follows: An average infected individuals makes contact sufficient to infect λN (here N is total population) others per unit time. Also, the probability that a given individual that each infected individual comes in contact with is susceptible is S/N. Thus, each infected individual causes \( \left( \lambda\,N \right) \left( S/N \right) = \lambda\,S \) infections per unit time. Therefore, I infected individuals cause a total number of infections per unit time of λSI. The γI term is even simpler to understand: γ is the fraction of infected individuals who recover (and re-enter the susceptible class) per unit time.

Using the relation \( S+I =1 , \) we can eliminate one dependent variable and reduce the problem to an ordinary differential equation

\[ \frac{{\text d}I}{{\text d}t} = \left[ \lambda - (\gamma + \mu ) \right] I(t) - \lambda\, I^2 , \qquad I(0) = I_0 . \]
For disease to spread, the derivative \( \dot{I} >0 \) must be positive. Since the above equation is of Bernoulli type, Mathematica can find its general solution without a problem.
DSolve[{y'[t] == a*y[t] - b*(y[t])^2, y[0] == i0}, y[t], t]
{{y[t] -> (a E^(a t) i0)/(a - b i0 + b E^(a t) i0)}}
Correspondingly, the SIS model without vital dynamics has the simpler form:
\[ \begin{cases} \dot{S} &= -\lambda\,I(t)\,S(t) + \gamma\, I(t) , \\ \dot{I} &= \lambda\,I(t)\,S(t) - \gamma\, I(t) , \end{cases} \qquad\mbox{subject to} \qquad \begin{cases} S(t) + I(t) =1, \\ S(0) = S_0, \quad I(0) = I_0 . \end{cases} \]
Now we plot solutions depending on the initial conditions.
sol = DSolve[{y'[t] == 1.1*y[t] - 1.4*(y[t])^2, y[0] == i0}, y[t], t];
toplot = Table[sol[[1, 1, 2]], {i0, 0.3, 1.1, 0.1}];
Plot[Evaluate[toplot], {t, 0, 6}, PlotStyle -> Thick]
The same graph could be obtained with the aid of Table to define the list of initial conditions, which we call inits, and then use Map to apply the function graph to the list of numbers inits. The result is a list of nine graphics objects that we name toshow.
graph[i0_] :=
numsol = NDSolve[{y'[t] == 1.1*y[t] - 1.4*(y[t])^2, y[0] == i0}, y[t], {t, 0, 10}];
Plot[y[t] /. numsol, {t, 0, 10}, DisplayFunction -> Identity]]
inits = Table[i, {i, 0.3, 1.1, 0.1}];
toshow = Map[graph, inits];
Show[toshow, DisplayFunction -> $DisplayFunction, PlotRange -> {0, 1}]
Example 7: Forms of a given element with different numbers of protons, neutrons, and its nuclear energy are called nuclides. Some nuclides are not stable. For example, potassium-40 (40K) naturally very clowly decays to reach argon-40 (40Ar). This decay was first observed, but not understood, by the Nobel laureate Henri Becquerel (1852--1908) in 1896. However, Marie Skłodowska Curie (1867--1934) began studying this decay in 1898, cointed out it as radioactivity, and discovered the radioactive substances polonium and radium.

Given a sample of 40K of sufficient size, after \( 1.2 \times 10^{9} \) years approximately half of the sample will have decayed to 40Ar. The half-life of a nuclide is the time for half the nuclei in a given sample to decay.

Suppose you know that the half-life of polonium-209 is about 125 years. What percent of the original amount of 209Po remains after 75 years?

Let y0 represent te original amount of polonium-209 that is present and y(t) denote the amount of 209Po after t years. Then y(t) is the solution to the following initial value problem:

\[ \dot{y} = -r\,y(t), \qquad y(0) = y_0 , \]
given that \( y(125) = y_0 /2 . \) Using separation of variables, we find the solution of the given IVP:
\[ \frac{{\text d}y}{y} = -r\,{\text d}t \qquad \Longrightarrow \qquad y(t) = y_0 e^{-rt} . \]
Using half-life value, we determine r to be
\[ y(125) = y_0 e^{-125\,r} = y_0 /2 \qquad \Longrightarrow \qquad e^{125\,r} =2 \qquad \Longrightarrow \qquad r = \ln 2 / 125 . \]
Substituting this value of r into the solution formula, we get
\[ y(t) = y_0 e^{- \ln 2\,t/125} = y_0 \left( e^{- \ln 2\,t} \right)^{1/125} = y_0 \left( e^{\ln 2^{-1/2}} \right)^{t/125} = y_0 \left( \frac{1}{2} \right)^{t/125} . \]
After 75 years, there will remain
\[ 2^{-75/125} = 2^{-3/5} \approx 0.659754 \]
or about 66\% of the original amount.


Return to Mathematica page
Return to the main page (APMA0330)
Return to the Part 1 (Plotting)
Return to the Part 2 (First Order ODEs)
Return to the Part 3 (Numerical Methods)
Return to the Part 4 (Second and Higher Order ODEs)
Return to the Part 5 (Series and Recurrences)
Return to the Part 6 (Laplace Transform)
Return to the Part 7 (Boundary Value Problems)