Preface


This section presents some examples to motivate the reader to study system of ordinary differential equations (ODE for short).

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 first course APMA0330
Return to the main page for the second course APMA0340
Return to Part II of the course APMA0340
Introduction to Linear Algebra with Mathematica

Motivation


This section presents motivating examples to study systems of ordinary differential equations (ODE for short),
\[ \begin{cases} \dot{y}_1 &= f_1 \left( y_1 , y_2 , \ldots , y_n , t \right) , \\ \dot{y}_2 &= f_2 \left( y_1 , y_2 , \ldots , y_n , t \right) , \\ \cdots & \qquad \cdots \\ \dot{y}_n &= f_n \left( y_1 , y_2 , \ldots , y_n , t \right) , \\ \end{cases} \]
or in vector form
\begin{equation} \label{EqMotiv.1} \dot{\bf y} = {\bf f} \left( {\bf y},t \right) , \end{equation}
where y(t) and f(y, t) are column-vectors:
\[ {\bf y} (t) = \begin{bmatrix} y_1 (t) \\ \vdots \\ y_n (t) \end{bmatrix} = \left[ y_1 (t), \ldots , y_n (t) \right]^{\mathrm{T}} , \qquad {\bf f}({\bf y}, t) = \begin{bmatrix} f_1 \left( y_1 , \ldots , y_n , t \right) \\ \vdots \\ f_n \left( y_1 , \ldots , y_n , t\right) \end{bmatrix} . \]
The system of differential equations \eqref{EqMotiv.1} occurs frequently in particular descriptions of physical phenomena from natural science. Autonomous systems have input functions that do not explicitly contain t: f = f(y). Using a natural extension of known statements for a single differential equation, we present the following existence and uniqueness theorem.
Theorem: Consider the initial value problem
\begin{equation} \label{EqMotiv.IVP} \dot{\bf y} = {\bf f} \left( {\bf y},t \right) , \qquad {\bf y}(0) = {\bf y}_0 , \end{equation}
where y0 ∈ ℝn. Suppose that f : ℝn → ℝn is a continuously differentiable function. Then there exists a unique solution of the given initial value problem \eqref{EqMotiv.IVP}.
Here overdot stands for the derivative with respect to time variable, \( \dot{y} = {\text d}y/{\text d}t . \)

Example 1: Consider the planar system of equations

\[ \begin{cases} \dot{x} = x^2 + y , \\ \dot{y} = y^2 + x . \end{cases} \]
     
VectorPlot[{y + x^2 , y^2 + x}, {x, -1, 3}, {y, -1, 3}, VectorMarkers -> "Drop"]
       Direction field for the given system of equations.            Mathematica code

Each solution x = x(t), y = y(t) of the equation above satisfying the initial condition \[ x(0) > 0 , \qquad y(0) > 0 , \] cannot exist on an interval 0 ≤ t < ∞

In fact, the trajectory of any solution initiating in the first quadrant is contained in the first quadrant. So its global solution does not exist.    ■

We start with motivating examples that you most likely saw in other courses.

Example 2: A projectile in vacuum is an object upon which the only force acting is gravity, so the influence of air resistance is negligible. More realistic motion under drag force will be considered in the next example. Projectile motion is a form of motion experienced by a launched object. Ballistics (Greek: βάλλειν, romanized: ba'llein, lit. 'to throw') is the science of dynamics that deals with the flight, behavior and effects of projectiles, especially bullets, unguided bombs, rockets, or the like; the science or art of designing and accelerating projectiles so as to achieve a desired performance. In this example, we consider only two dimensional projectile motion in a vacuum. An object dropped from rest is a projectile; an object that is thrown vertically upward is also a projectile. So we first consider a projectile that is thrown upward at an angle α to the horizontal line.
     
       Vertical trajectories.            Projectile thrown upward at angle α to abscissa.

Consider an object of mass m which is launched from the origin in a uniform gravitational field g with an initial launch speed of v0 at an angle of α to the horizontal (abscissa). For simplicity, it will be assumed that the motion of the projectile is under the influence of gravity only and that its flight is over level ground so that the end of the flight is at the same level as the beginning. The position vector of the projectile at time t after firing is defined by the vector

\[ {\bf r} (t) = \left[ x(t), y(t) \right] = x(t)\,{\bf i} + y(t)\,{\bf j} \]
and was set into motion with the initial velocity
\[ \left( v_0 \cos\alpha, v_0 \sin\alpha \right) = v_0 \cos\alpha\,{\bf i} + v_0 \sin\alpha\,{\bf j} . \]
The equations of motion are
\[ \begin{split} \ddot{x} &= \phantom{-}0 , \qquad x(0) = 0, \quad \dot{x}(0) = v_0 \cos\alpha \\ \ddot{y} &= -g , \qquad y(0) =h, \quad \dot{y}(0) = v_0 \sin\alpha , \end{split} \tag{2.1} \]
where g is the acceleration due to gravity. Solving these equations, we obtain
\[ \begin{split} x(t) &= v_0 t\,\cos\alpha , \\ y(t) &= v_0 t\,\sin \alpha - \frac{g\,t^2}{2} , \end{split} \tag{2.2} \]
      One can eliminate t from the equations above and get \begin{equation*} y(x) = \left( \tan\alpha \right) x - \frac{g\,\sec^2 \alpha}{2\,v_0^2} \, x^2 . \tag{2.3} \end{equation*} The equation above represents a quadratic function: thus, its graph is a parabola. Therefore, in the absence of air resistance, the path of the projectile is a parabola.
g[x_, a_] = x*(Tan[a]) - (9.81/2)*x^2/Cos[a]^2;
Plot[Table[g[x, c], {c, 0, Pi/2, 0.3}], {x, 0, 0.1}, PlotStyle -> Thick, ImageSize -> 250, PlotRange -> {0, 0.056}]
       Trajectories of a projectile.            Mathematica code

The time T taken by the projectile to reach the maximum height H, and the range R of the projectile are given by
\[ \begin{split} T &= \frac{v_0}{g}\,\sin\alpha , \\ H &= \frac{v_0^2}{2g}\,\sin^2 \alpha , \\ R &= \frac{v_0^2}{g}\,\sin (2 \alpha ) = v_0^2 \cos \alpha \, \frac{\sin \alpha + \left( \sin^2 \alpha + 2gh/v_0^2 \right)^{1/2}}{g} , \end{split} \tag{2.4} \]
where h is a height above ground at which the object was initially released, so y(0) = h. Usually, the initial speed v0 is not known and may be inconvenient to measure. The angle α which will give the maximum range is determined from the equation
\[ \sin^2 \alpha = \left( 2 + 2gh/ v_0^2 \right)^{-1} \]
      Let (xm, ym) be coordinates of the vertex corresponding to the parabola solution curve, so ym is the maximum height attained and xm is the corresponding horizontal coordinate. \[ x_m = \frac{v_0^2}{2g} \,\sin \left( 2\alpha \right) , \qquad y_m = \frac{v_0^2}{2g} \,\sin^2 \alpha \]
ParametricPlot[{Sin[2*a], (Sin[a])^2}, {a, 0, Pi/2}, {r, 0, 1}, Mesh -> Automatic, MeshStyle -> Directive[Thick, Red]]
       Trajectory of the vertex.            Mathematica code

      Now we plot trajectories and their corresponding vertices together.
g[x_, a_] = x*(Tan[a]) - (9.81/2)*x^2/Cos[a]^2;
plot = Plot[Table[g[x, c], {c, 0, Pi/2, 0.3}], {x, 0, 0.1}, PlotStyle -> Thick, ImageSize -> 250, PlotRange -> {0, 0.056}];
ver = ParametricPlot[{Sin[2*a]/2/9.81, (Sin[a])^2/2/9.81}, {a, 0, Pi/2}, {r, 0, 1}, Mesh -> Automatic, MeshStyle -> Directive[Thick, Red]];
Show[plot, ver]
       Projectile trajectories along with their vertices.            Mathematica code

 

  1. de Alwise, T., Projectile motion with Mathematica, International Journal of Mathematical Education in Science and Technology, 2000, Vol. 31, No. 5, pp.749--755. doe: 10.1080/002073900434413
  2. Benacka, I., Stubna, I., Ball launched against an inclined plane---an example of using recurrent sequences in school physics, International Journal of Mathematical Education in Science and Technology, 2009, Vol.40, No. 5, pp. 696--705.
  3. Bose, S.K., Thoughts on projectile motion, American Journal of Physics, 1985, Vol. 53, No. 2, pp. 175
  4. de Mestre, N., The Mathematics of Projectiles in Sport, Cambridge University Press, 2012, https://doi.org/10.1017/CBO9780511624032
  5. Donnelly, D., The parabolic envelope of constant initial speed trajectories, American Journal of Physics, 1992, Vol. 60, No. 12, pp. 1149--1150.
  6. Fernández-Chapou, J.L., Salas-Brito, A.L. and Vargas, C.A., An elliptic property of parabolic trajectories, American Journal of Physics, 2004, Vol. 72, pp. 1109--1110.
  7. Hu, H., Yu, J., Another look at Projectile motion, The Physics Teacher, 2000, Vol. 38, No. 10, page 423.
  8. Padyala, R., An alternative view of the elliptic property of a family of projectile paths, The Physics Teacher, 2019, Vol. 57, No. 9, pp376--377.
  9. Sarafian, H., On projectile motion, The Physics Teacher, 1999, Vol. 37, No. 2, [[. 86--88.
   ■

Example 3: A particular practical interest is a situation when the projectile moves under gravity in a medium whose resistance is not negligible. The drag force acting on the projectile depends on many factors, including material of the projectile, its shape, speed, and some other parameters. Most of these factors are incorporated in the Reynolds number, which we denote by Re to avoid confusion with "Re" (it is a custom to denote as ℜ), which serves to identify a "real part" of a complex number (a standard notation for the Reynolds number is Re). The Reynolds number is the ratio of inertial forces to viscous forces within a fluid that is subjected to relative internal movement due to different fluid velocities. This number was introduced by Sir George Stokes (1819--1903), and popularized by Osborne Reynolds (1842--1912).

Besides gravitational and resistance force, there is another force acting on a ball traveling in air or liquid---the aerodynamic lift due to the backspin of the ball. The air pressure is decreased above the ball and increased below the ball, providing an aerodynamic lift force (Kutta–Joukowski theorem, also Zhukovsky). Therefore, analysis of ball movement in air requires utilization of three forces: gravity, drag, and lift. For a golf ball traveling with speed in range 42 -- 70 m/sec (Reynolds numbers from about 1.25×105 to 1.9×105), the lift force varies approximately linearly with the backspin angular speed.

The air resistance force depends on many things: the speed and direction of motion of the object (called projectile), its shape, size, whether or not it is spinning, the properties of air, whether or not there is any wind, etc.. Any attempt to construct a model of air resistance to be used in projectile motion treating all of these variables rigorously is a tremendously hard job. Since an accurate determination of the resistive force is very complicated, we concentrate our attention on a projectile having shape of a sphere. Fortunately, the dependents of the air resistive force acting on a ball is known and it depends on projectile's speed. For a sphere of radius r moving a fluid of density ρ and viscosity η, the drag force and the Reynolds number are

\[ \begin{split} F = 6\pi \eta\rho v \qquad \mbox{for} \quad R_e < 1, \\ R_e = 2r\rho v/\eta . \end{split} \]
For a ball with a density of 1 g/cm³, and using the value of ratio η/ρ = 0.15 cm²/sec for air, one finds utilization of a linear grad force appropriate when
\[ r \le 4 \times 10^{-3} \mbox{ cm}. \]
Spheres with this small of a radius would seem to be of little interest. For example, a golf ball (42.67 mm in diameter) has a resistive force with is proportional to the speed raised to power 1.3, so Fv1.3. For spheres with radii and speeds of practical interest a reasonable approximation to the drag force is
\[ F = \pi r^2 v^2 /4 \qquad \mbox{for} \quad 1 \ll R_e \le 10^5 . \]
      The dependence of the drag force on velocity of a ball moving in air is well documented. Experiments measuring drag force dependence on the ball velocity were conducted by the faculty at MIT, which shows a complicated behavior. For small velocities, it is approximately linear, but then drag force starts closer to the quadratic function with speed increases. At large Reynolds numbers, we observe turbulence behind the ball. So when speed increases (and the Reynolds number exceeds 105), the drag force declines and oscillates as shown in the figure. Such phenomena is observed in baseball (73–75 mm in diameter with a mass of 142 to 149 g) and volleyball (20.7--21.3 cm in diameter with a mass of 260-280 g).
f[t_, y_] = Sin[2.35*t*y]*Cos[t + 1.0*y];
s = NDSolve[{y'[t] == f[t, y[t]], y[0] == 1}, y, {t, 0, 30}];
Plot[Evaluate[y[t] /. s], {t, 3.9, 25}, PlotRange -> {0.2, 1}, PlotTheme -> "Web"]
       Figure 6: Dependence of the drag force on ball's speed.            Mathematica code

Usually the drag force is modeled by a power function depending on projectile velocity v:
\[ F({\bf v}) \sim \left( a + b\,\| {\bf v} \| \right) {\bf v}^p \tag{3.1} \]
with some coefficients 𝑎 and b. When the power p of the speed is an integer, there is an analytic solution involving the speed and the direction of the projectile implicitly in terms of logarithmic and trigonometric functions. When the power is not an integer, the relation between the projectile's speed and direction at any instant must be obtained by numerical methods. In both cases the time, range and height of the projectile must be obtained numerically.

 

Quadratic Resistance and Speed


The equations of motion when drag force is proportional to the square of speed are
\[ \begin{split} m\,\frac{{\text d} v_x}{{\text d}t} &= -b\,v_x \left( v_x^2 + v_y^2 \right)^{1/2} , \\ m\,\frac{{\text d} v_y}{{\text d}t} &= -mg -b\,v_y \left( v_x^2 + v_y^2 \right)^{1/2} , \end{split} \]
where the vertical y-direction is taken upward and the constant b (depending on the radius of the ball) represents the coefficient in the resistance force.

Under the conditions of no drift, no wind, an inertial coordinate frame, the density and velocity of sound as functions of the height only, plus the assumption that the only forces with any appreciable effect are gravity and the drag, it is possible to write the equations for a trajectory in the form

\begin{equation} m\,\ddot{x} = -\rho\,d^2 v \,K_D \left( \frac{v}{a} \right) \dot{x} , \qquad m\,\ddot{y} = -\rho\,d^2 v \,K_D \left( \frac{v}{a} \right) \dot{y} - mg . \end{equation}
with
\[ v = \sqrt{\dot{x}^2 + \dot{y}^2} , \qquad \rho = \rho (y) , \qquad a= a(y), \]
where (x, y) are the coordinates in the horizontal and vertical directions, m is the projectile's mass, d its diameter, \( K_D (v/a) \) is the dimensionless drag coefficient, ρ is the density of air and 𝑎 is the speed of sound in air, which depends on the type of medium and the temperature of the medium. Upon performing numerical experiments, we can determine the angle you need for maximum projectile range.

     
a=0.3;
s3 = NDSolve[{x''[t] + x'[t]*Sqrt[(x'[t])^2 + (y'[t])^2] == 0, y''[t] + 1 + y'[t]*Sqrt[(x'[t])^2 + (y'[t])^2] == 0, x[0] == 0, x'[0] == Cos[a], y[0] == 0, y'[0] == Sin[a]}, {x[t], y[t]}, {t, 0, 1}];
p3 = ParametricPlot[{x[t], y[t]} /. s3, {t, 0, 0.54}, PlotStyle -> Thickness[0.01], PlotRange -> {{0, 0.55}, {0, 0.33}}];
a=0.5;
s5 = NDSolve[{x''[t] + x'[t]*Sqrt[(x'[t])^2 + (y'[t])^2] == 0, y''[t] + 1 + y'[t]*Sqrt[(x'[t])^2 + (y'[t])^2] == 0, x[0] == 0, x'[0] == Cos[a], y[0] == 0, y'[0] == Sin[a]}, {x[t], y[t]}, {t, 0, 1}];
p5 = ParametricPlot[{x[t], y[t]} /. s5, {t, 0, 0.84}, PlotStyle -> {Purple, Thickness[0.01]}, PlotRange -> {{0, 0.55}, {0, 0.33}}];
a=0.7;
s7 = NDSolve[{x''[t] + x'[t]*Sqrt[(x'[t])^2 + (y'[t])^2] == 0, y''[t] + 1 + y'[t]*Sqrt[(x'[t])^2 + (y'[t])^2] == 0, x[0] == 0, x'[0] == Cos[a], y[0] == 0, y'[0] == Sin[a]}, {x[t], y[t]}, {t, 0, 2}];
p7 = ParametricPlot[{x[t], y[t]} /. s7, {t, 0, 1.1}, PlotStyle -> {Green, Thickness[0.01]}, PlotRange -> {{0, 0.6}, {0, 0.33}}]];
a=1,0;
s10 = NDSolve[{x''[t] + x'[t]*Sqrt[(x'[t])^2 + (y'[t])^2] == 0, y''[t] + 1 + y'[t]*Sqrt[(x'[t])^2 + (y'[t])^2] == 0, x[0] == 0, x'[0] == Cos[a], y[0] == 0, y'[0] == Sin[a]}, {x[t], y[t]}, {t, 0, 2}];
p10 = ParametricPlot[{x[t], y[t]} /. s10, {t, 0, 1.41}, PlotStyle -> {Red, Thickness[0.01]}, PlotRange -> {{0, 0.6}, {0, 0.33}}]
a=1,2;
s12 = NDSolve[{x''[t] + x'[t]*Sqrt[(x'[t])^2 + (y'[t])^2] == 0, y''[t] + 1 + y'[t]*Sqrt[(x'[t])^2 + (y'[t])^2] == 0, x[0] == 0, x'[0] == Cos[a], y[0] == 0, y'[0] == Sin[a]}, {x[t], y[t]}, {t, 0, 10}];
p12 = ParametricPlot[{x[t], y[t]} /. s12, {t, 0, 1.555}, PlotStyle -> {Black, Thickness[0.01]}, PlotRange -> {{0, 0.6}, {0, 0.33}}];
Show[p12, p10, p7, p5, p3]
       Figure 7: Trajectories of a projectile depending on the initial angle.            Mathematica code

Now we compare all models when a ball is thrown under the same angle, 45°, and the same speed.
     
a = Pi/4;
ss = NDSolve[{x''[t] + x'[t]*(1 + 0.4*Sqrt[(x'[t])^2 + (y'[t])^2]) == 0, y''[t] + 1 + y'[t]*(1 + 0.4*Sqrt[(x'[t])^2 + (y'[t])^2]) == 0, x[0] == 0, x'[0] == Cos[a], y[0] == 0, y'[0] == Sin[a]}, {x[t], y[t]}, {t, 0, 1}];
pp = ParametricPlot[{x[t], y[t]} /. ss, {t, 0, 2.5}, PlotStyle -> {Orange, Thickness[0.01]}, PlotRange -> {{0, 0.6}, {0, 0.25}}];
s = NDSolve[{x''[t] + x'[t] == 0, y''[t] + 1 + y'[t] == 0, x[0] == 0, x'[0] == Cos[a], y[0] == 0, y'[0] == Sin[a]}, {x[t], y[t]}, {t, 0, 1}];
p = ParametricPlot[{x[t], y[t]} /. s, {t, 0, 2.5}, PlotStyle -> {Black, Thickness[0.01]}, PlotRange -> {{0, 0.6}, {0, 0.25}}];
s3 = NDSolve[{x''[t] + Sqrt[(x'[t])^2 + (y'[t])^2]*Abs[x'[t]]^(0.3) == 0, y''[t] + 1 + Sqrt[(x'[t])^2 + (y'[t])^2]*Abs[y'[t]]^(0.3) == 0, x[0] == 0, x'[0] == Cos[a], y[0] == 0, y'[0] == Sin[a]}, {x[t], y[t]}, {t, 0, 1.5}];
p3 = ParametricPlot[{x[t], y[t]} /. s3, {t, 0, 105}, PlotStyle -> {Blue, Thickness[0.01]}, PlotRange -> {{0, 0.6}, {0, 0.25}}];
s2 = NDSolve[{x''[t] + x'[t]*Sqrt[(x'[t])^2 + (y'[t])^2] == 0, y''[t] + 1 + y'[t]*Sqrt[(x'[t])^2 + (y'[t])^2] == 0, x[0] == 0, x'[0] == Cos[a], y[0] == 0, y'[0] == Sin[a]}, {x[t], y[t]}, {t, 0, 1}];
p2 = ParametricPlot[{x[t], y[t]} /. s2, {t, 0, 2.5}, PlotStyle -> {Purple, Thickness[0.01]}, PlotRange -> {{0, 0.6}, {0, 0.25}}];
s0 = NDSolve[{x''[t] == 0, y''[t] + 1 == 0, x[0] == 0, x'[0] == Cos[a], y[0] == 0, y'[0] == Sin[a]}, {x[t], y[t]}, {t, 0, 1}];
p0 = ParametricPlot[{x[t], y[t]} /. s0, {t, 0, 2.5}, PlotStyle -> {Green, Thickness[0.01]}, PlotRange -> {{0, 0.6}, {0, 0.25}}];
Show[pp, p, p3, p2, p0, PlotRange -> {{0, 0.6}, {0, 0.25}}]

The graphs of trajectories correspond to the following models:
\( \begin{cases} \ddot{x} + \dot{x} \left( 1 + 0.4\,\sqrt{\dot{x}^2 + \dot{y}^2} \right) &=0 , \\ \ddot{y} +1+ + \dot{y} \left( 1 + 0.4\,\sqrt{\dot{x}^2 + \dot{y}^2} \right) &=0 \end{cases} \)     (orange),
\( \begin{split} \ddot{x} + \dot{x} =0 , \\ \ddot{y} +1+ + \dot{y} =0 \end{split} \)     (black),
\( \begin{cases} \ddot{x} + \sqrt{\dot{x}^2 + \dot{y}^2}\, \dot{x}^{0.3} &=0 , \\ \ddot{y} +1+ + \sqrt{\dot{x}^2 + \dot{y}^2}\, \dot{y}^{0.3} &=0 \end{cases} \)     (blue),
\( \begin{cases} \ddot{x} + \dot{x} \,\sqrt{\dot{x}^2 + \dot{y}^2} &=0 , \\ \ddot{y} +1+ + \dot{y} \,\sqrt{\dot{x}^2 + \dot{y}^2} &=0 \end{cases} \)     (purple),
and finally, the trajectory of projectile in vacuum (green).

 

Linear Resistance and Speed


In simple linear resisting medium, the drag force at time t, which acts on the projectile, is taken to be directly proportional to the instantaneous velocity of the projectile at that time and directed in the opposite direction. From Newton’s second law, we have for the horizontal and vertical components of the motion to be

\[ \begin{split} \ddot{x} &= -\gamma\,\dot{x} , \\ \ddot{y} &= -\gamma\,\dot{y} -g , \end{split} \]
where γ is a positive drag coefficient per unit mass and is assumed to remain constant during the motion of the projectile. This system of differential equations is dec oupled and each equation can be solved independently of another one.
DSolve[{x''[t] + gamma*x'[t] == 0, x[0] == 0, x'[0] == v*Cos[alpha]}, x[t], t]
{{x[t] -> (E^(-gamma t) (-1 + E^(gamma t)) v Cos[alpha])/gamma}}
DSolve[{y''[t] + gamma*y'[t] + g == 0, y[0] == 0, y'[0] == v*Sin[alpha]}, y[t], t]
{{y[t] -> -(( E^(-gamma t) (g - E^(gamma t) g + E^(gamma t) g gamma t + gamma v Sin[alpha] - E^(gamma t) gamma v Sin[alpha]))/ gamma^2)}}
Since coefficient γ is per unit mass, inclusion of the mass here leads to a slight simplification in the resulting equations of motion. Solving these two differential equations in the usual way, on applying the initial conditions yields the following well-known equations of motion
\[ \begin{split} x(t) &= \frac{v_0 \cos\alpha}{\gamma} \left( 1 - e^{-\gamma\,t} \right) \\ y(t) &= \left( \frac{g}{\gamma^2} + \frac{v_0 \sin\alpha}{\gamma} \right) \left( 1 - e^{-\gamma\,t} \right) - \frac{gt}{\gamma} . \end{split} \]
The trajectory equation for a projectile in a linear resisting medium can be found by eliminating the time t:
\[ y = \left( \frac{g\,\sec \alpha}{v_0 \gamma} + \tan \alpha \right) x + \frac{g}{\gamma^2} \,\ln \left( 1- \frac{\gamma\, \sec \alpha}{v_0}\, x \right) . \]
      Now we plot trajectories in the presence of linear drag force.
gamma = 0.5;
g[x_, a_] = x*(Tan[a] + 9.81/Cos[a]/gamma) + 9.81/gamma^2*Log[1 - gamma*x/Cos[a]];
Plot[Table[g[x, c], {c, 0, Pi/2, 0.3}], {x, 0, 0.1}, PlotStyle -> Thick, ImageSize -> 250, PlotRange -> {0, 0.055}]
       Figure 7: Trajetories in the presence of linear resistance.            Mathematica code

Considering the value at the maximum, via dy/dx = 0, we obtain
\[ x_m = \rho \,\frac{\sin\alpha\,\cos\alpha}{1+\varepsilon\,\sin\alpha} , \qquad y_m = \frac{\rho}{\varepsilon^2} \left[ \varepsilon\,\sin\alpha - \ln \left( 1+\varepsilon\,\sin\alpha \right) \right] , \]
where the dimensionless parameters are
\[ \varepsilon = v_0 \gamma /g, \qquad \rho = v_0^2 /g. \]
In the absence of air resistance, the trajectory is parabolic in shape and determined by just two parameters, the initial launch speed v0 and the angle of projection α. The path of the trajectory in this case is therefore symmetrical about its peak. On the other hand, in the linearly resistive model the trajectory of the projectile becomes asymmetrical as its shape is modified by the presence of an additional parameter introduced through the damping coefficient γ. Here the peak in the trajectory occurs closer to its final impact point than its initial point of projection. This forward blunting behavior in the trajectory paths is characteristic of projectiles in a linear resisting medium and is one of a number of results which will be confirmed analytically from the closed-form solutions we obtain for the time of flight and range.

It is important to know where the trajectories peak is in a linear resisting medium. Equating the rime derivative of y(t) to zero and solving for t, we determine the time tm taken to reach its maximum height:
\[ t_m = \frac{1}{\gamma}\,\ln \left( 1 + \frac{\gamma v_0}{g}\,\sin \alpha \right) . \]
If we introduce the dimensionless quantity c = γv0/g and write ζ(α) = 1 + csin(α), then the equation for tm can be rewritten more compactly as
\[ t_m = \frac{1}{\gamma}\,\ln \zeta\left( \alpha\right) = \frac{v_0 \sin\alpha}{g} \cdot \frac{\ln \zeta (\alpha )}{\zeta (\alpha ) -1} . \]
Then the coordinates of the vertex become
\[ \begin{split} x_m (\alpha ) &= \frac{v_0^2 \sin (2 \alpha )}{2g} \,\frac{1}{\zeta (\alpha )} \\ y_m (\alpha ) &= \frac{v_0^2 \sin^2 \alpha}{g} \,\frac{\zeta (\alpha ) -1 - \ln \zeta (\alpha )}{\left( \zeta (\alpha ) -1 \right)^2} . \end{split} \]
So the maximum height of the projectile is ym(α) in the linearly resistive case. The curve connecting the peak heights in the trajectories of a projectile in a linear resisting medium turn out to be elliptic-like,

When a projectile is launched vertically upwards (i.e. α = π/2), the horizontal range is xm(π/2) = 0 and its peak heights becomes

\[ y_m (\pi /2) = \frac{v_0}{\gamma} - \frac{g}{\gamma^2} \,\ln \left( 1 + \frac{\gamma v_0}{g} \right) . \]
An explicit formula for the range of the projectile is available only in terms of a special function---the Lambert W function, also called the omega function or product logarithm. Its origins date back to the mid-eighteenth century in the initial work of Johann Lambert in 1758 and subsequent work of Leonhard Euler in 1783. Its modern day definition was given in 1980s by a group of mathematicians. More information can be found in a special section on Lambert function.

The Lambert W function is defined to be the inverse of the function \( w\,e^w = z, \) and is therefore the solution of the defining Lambert equation

\[ W(z) \,e^{W(z)} = z . \]
In general, for a given complex z, the defining equation will have an infinite number of solutions and is therefore a multivalued analytic function. Its principal branch is denoted by W0(z); other branches are labeled with Wk(z), k ∈ ℤ. However, only two of them are real-valued functions that correspond k = 0 and k = −1. So one real solution of the defining Lambert equation is W0(x) or simply W(x) for −1/ex < ∞, while the secondary real branch satisfying W(x) ≤ −1 is denoted by W-1(x), −1/ex < 0.

The time of flight T is found by setting y = 0 in the formula for y(t), solution of the system of differential equations. Then solving for the time. When this is done, after a little algebraic rearranging of terms

\[ u = \zeta \left( 1 - e^{-u} \right) , \]
where we have set u = γT. The time of flight T is therefore obtained as the solution to the transcendental equation above. Rewriting this equation as
\[ -\zeta\,e^{-\zeta} = \left( u - \zeta \right) e^{u-\zeta} \]
we see that its solution is expressed in terms of Lambert W function
\[ T = \frac{v_0 \sin\alpha}{g} \left( \frac{\zeta (\alpha ) + W_0 \left( -\zeta (\alpha )\,e^{-\zeta (\alpha )}\right)}{\zeta (\alpha ) -1} \right) . \]
upon substituting for u and recasting into the zeta function form. The relationship between the time of flight of a projectile in a linear resisting medium and the Lambert W function is apparent.

An exact analytic expression for the range R of a projectile in a linear resisting medium can also be found in terms of the Lambert W function.

\[ R= \frac{v_0^2 \sin 2\alpha}{2g} \left( \frac{1 - \exp \left\{ -\zeta (\alpha ) - W_0 \left( -\zeta (\alpha )\,e^{-\zeta (\alpha )}\right) \right\}}{\zeta (\alpha ) -1} \right) . \]
It is important to recognize that the form of expressions involving the Lambert W function are not unique. The expression for the range given above can be re-written in the equivalent form of
\[ R= \frac{v_0^2 \sin 2\alpha}{2g} \left( \frac{\zeta (\alpha ) + W_0 \left( -\zeta (\alpha )\,e^{-\zeta (\alpha )}\right\}}{\zeta (\alpha ) \left( \zeta (\alpha ) -1 \right)} \right) . \]

The angle αmax that maximizes the range of the projectile is found from the angle that satisfies dR/dt = 0. With the availability of an analytical expression for the range as a function of the launch angle α in terms of the Lambert W function, a function whose derivative is known, it is tempting to differentiate one of the equivalent expressions for the range, set the result equal to zero and solve for αmax. This yields

\[ \alpha_{\max} = \arcsin \left[ \frac{c\,W_0 \left( c^2 - 1/e \right)}{c^2 -1 - W_0 \left( (c^2 -1 )/e \right)} \right] , \qquad c = \frac{\gamma v_0}{g} \ne 1 \]
The angle which maximizes the range of a projectile in a linear resisting medium can therefore be expressed explicitly as a function of the dimensionless drag coefficient c = γv0/g in terms of the Lambert W function. Note that for the special case of c = 1, the angle which maximizes the range of the projectile is given by
\[ \alpha_{\max} = \arcsin \left[ \frac{1}{e-1} \right] . \]

The peak trajectories curve is

\[ y_m = \frac{\gamma^2}{g} \left[ c\,\sin \alpha - \ln \left( 1+ c\,\sin\alpha \right) \right] , \qquad c = \frac{\gamma v_0}{g} \ne 1 \]
It is another transcendental equation which can be solved for α in terms of the Lambert W function.

 

  1. Brancazio, P.J., Trajectory of a fly ball, The Physics Teacher, 1985, Vol. 23, No. 1, pp. 20-23. https://doi.org/10.1119/1.234170
  2. Coutis, P., Modelling the projectile motion of a cricket ball, International Journal of Mathematical Education in Sciences and Technology, 1998, Vol. 29, No. 6, pp. 789--798.https://doi.org/10.1080/0020739980290601
  3. Erlichson, H., Maximum projectile range with drag and lift, with particular application to golf, American Journal of Physics, 1983, Vol. 51, No. 4, pp. 357--362. https://doi.org/10.1119/1.13248
  4. Groetsch, C., Cipra, B., Halley's comment—projectiles with linear resistance, Mathematics Magazine, 1997, Vol. 70, No. 4, pp. 273--280. doi: 10.1080/0025570X.1997.11996553
  5. Hackborn, W.W., Motion through air: what a drag, Canadian Applied Mathematics Quarterly, 2006, vol. 14, pp. 285–298.
  6. Hayen, J.C>, Projectile motion in a resistant medium (part I: exact solution and properties), International Journal of Non-Linear Mechanics, 2003, Vol. 38, No. 3, pp.357--369.
  7. Hernández-Saldaña, H., On the locus formed by the maximum heights of projectile motion with air resistance, European Journal of Physics, 2010, Vol. 31, No. 8, 1319--1329. doi:10.1088/0143-0807/31/6/00
  8. Kantrowitz, R., The English Galileo and his vision of projectile motion under air resistance, International Journal of Mathematics and Mathematical Sciences, 2020, Article ID 9695053 | https://doi.org/10.1155/2020/9695053
  9. de Mestre, N., The Mathematics of Projectiles in Sport, Cambridge University Press, 2012, https://doi.org/10.1017/CBO9780511624032
  10. Mohazzabi, P., When does air resistance become significant in projectile motion? The Physics teacher, 2018, Vol. 68, pp.168--169.
  11. Packel E and Yuen D 2004 Projectile motion with resistance and the Lambert function, Coll. Math. J., Vol. 35, pp. 337
  12. Parker, G.W., Projectile motion with air resistance quadratic in the speed, American Journal of Physics, 1977, Vol. 45, No. 7, pp. 606--610.
  13. Stewart, S.M., An analytic approach to projectile motion in a linear resisting medium, International Journal of Mathematical Education in Science and Technology, 2006, Vol. 37, No. 4, pp. 411-431. doi: 10.1080/00207390600594911
  14. Warburton, R.D.H. and Wang, J., Analysis of asymptotic motion with air resistance using the Lambert W function, American Journal of Physics, 2004, Vol. 72, pp. 1404--1407.
  15. Yabushita, K., Yamashita, M., and Tsuboi, K., An analytic solution of projectile motion with the quadratic resistance law using the homotopy analysis method, Journal of Physics A Mathematical and Theoretical 2007, Vol. 40, No 29:8403--8416. doi: 10.1088/1751-8113/40/29/015
   ■
Our next example is from mechanics and it requires the application of Euler--Lagrange equations for a conservative case:
\begin{equation} \label{EqMotiv.2} \frac{\text d}{{\text d}t}\,\frac{\partial {\cal L}}{\partial \dot{q}_i} \left( t, {\bf q}(t), \dot{\bf q}(t) \right) = \frac{\partial {\cal L}}{\partial q_i} \left( t, {\bf q}(t), \dot{\bf q}(t) \right) , \qquad i=1,2,\ldots n, \end{equation}
where q = (q1, q2, … , qn) are generalized (canonical) coordinates and the Lagrangian is the difference between the kinetic and potential energy:
\begin{equation} \label{EqMotiv.3} {\cal L} = {\text K} - \Pi . \end{equation}
For a simple autonomous system \( \dot{\bf y} + {\bf f}({\bf y}) = 0, \) with K = \( \dot{\bf y}^2 = \dot{y}^2_1 + \cdots + \dot{y}^2_n \) and \( \Pi = {\bf f}^2 ({\bf y}) , \) we have
\[ \pi_y = \frac{\partial {\cal L}}{\partial \dot{\bf y}} = \left[ \frac{\partial {\cal L}}{\partial \dot{y}_1} , \ldots , \frac{\partial {\cal L}}{\partial \dot{y}_n} \right] = 2\,\dot{\bf y}, \qquad \frac{\partial {\cal L}}{\partial {\bf y}} = 2{\bf f}({\bf y}) \cdot {\bf f}' ({\bf y}) , \]
where πy is the moment canonically conjugated to the variable y. The corresponding Hamiltonian is
\[ {\cal H} = {\mbox K} + \Pi = \pi_y \,\dot{\bf y} - {\cal L} = \dot{\bf y}^2 + {\bf f}^2 ({\bf y}) . \]

Example 4: We start with modeling a simple pendulum, which consists of a body (bob) suspended from a fixed point (pivot) so that it can swing back and forth under the influence of gravity. We consider the simple pendulum that is characterized by the following assumptions:

  • The bob is free to move within a plane (so we consider only two-dimensional oscillations).
  • The system is conservative, so the pendulum rotates in a vacuum and friction in the pivot is negligible.
  • The bob of mass m is attached to one end of a rigid, but weightless rod of length ℓ, which is assumed to be constant during the pendulum motion. The other end of the rod (or rigid spring) is supported at the pivot.
First, we illustrate with Mathematica the motion of the simple pendulum. Introducing the function f for circle arrow
start = -\[Pi]/2; \[Theta] = start + Pi/6;
f[s_Circle] := s /. Circle[a_, r_, {start_, end_}] :> ({s, Arrow[{# - r/10^6 {-Sin@end, Cos@end}, #}]} &[ a + r {Cos@end, Sin@end}])
we proceed
arc = Graphics@f[Circle[{0, 0}, 1, {-Pi/2, -1.0}]];
polygon = Polygon[{{-1, 0}, {1, 0}, {1, 1/4}, {-1, 1/4}}];
top = Graphics[{Gray, polygon}];
circle1 = Graphics[Circle[{0, -1/20}, 1/20]];
circle2 = Graphics[Circle[{2, -3}, 1/5]];
circle3 = Graphics[{Dashed, Circle[{0, 0}, 3.6, {-1.57, -1.0}]}];
line1 = Graphics[{Dashed, Line[{{0, -1/20}, {0, -3.8}}]}];
line2 = Graphics[{Thick, Line[{{0.03, -0.07}, {1.9, -2.83}}]}];
arrow = Graphics[{Red, Arrowheads[0.08], Arrow[{{1.986, -3.0}, {1.986, -4.6}}]}]; line3 = Graphics[{Line[{{2.2, -3.8}, {2.9, -3.8}}]}];
arrow2a = Graphics[{Arrowheads[0.04], Arrow[{{0.7, -3.8}, {0, -3.8}}]}];
arrow2b = Graphics[{Arrowheads[0.04], Arrow[{{1.3, -3.8}, {1.98, -3.8}}]}];
text1 = Graphics[ Text[Style["L sin\[Theta]", FontSize -> 14, Purple], {1.0, -3.8}]];
text2 = Graphics[ Text[Style["mg", FontSize -> 14, Purple], {2.1, -4.8}]];
line3 = Graphics[{Line[{{1.0, -3.6}, {2.9, -3.6}}]}];
line4 = Graphics[{Line[{{2.2, -3.0}, {2.9, -3.0}}]}];
text3 = Graphics[ Text[Style["L(1 - cos\[Theta])", FontSize -> 14, Purple], {2.7, -3.3}]];
arrow3a = Graphics[{Arrowheads[0.04], Arrow[{{2.6, -3.35}, {2.6, -3.6}}]}];
arrow3b = Graphics[{Arrowheads[0.04], Arrow[{{2.6, -3.27}, {2.6, -3.0}}]}];
text4 = Graphics[ Text[Style["L", FontSize -> 14, Purple], {1.2, -1.5}]];
text5 = Graphics[ Text[Style["\[Theta]", FontSize -> 16, Purple], {0.2, -0.7}]];
text6 = Graphics[ Text[Style["O", FontSize -> 14, Purple], {-0.16, -0.2}]];
textT = Graphics[ Text[Style["T", FontSize -> 14, Black], {1.15, -2.1}]];
arrowT = Graphics[{Blue, Arrowheads[0.08], Arrow[{{1.88, -2.8}, {1.2, -1.8}}]}];
Show[top, circle1, circle2, circle3, line1, line2, line3, arrow, arrow2a, arrow2b, text1, text2, line4, text3, arrow3a, arrow3b, text4, text5, text6, arrowT, textT, arc]
      The modeling of the motion is greatly simplified when the given body (bob) is considered essentially as a point particle. The position of the bob is described by the angle θ between the rod and the downward equilibrium vertical position, with the counterclockwise direction taken as positive. The only force acting on the pendulum is the gravitational force mg, acting downward, where g denotes the acceleration due to gravity. The position of the bob can be determined in Cartesian coordinates as
\[ x = \ell \,\sin \theta , \qquad y = -\ell\,\cos\theta , \]
where the origin is taken at the pivot and the positive vertical direction is upward. Using \( {\cal L} = \mbox{K} - \Pi , \) the Lagrangian, which is the difference of the kinetic energy K and the potential energy Π of the system, we have
\[ \frac{\text d}{{\text d}t} \,\frac{\partial {\cal L}}{\partial \dot{\theta}} = \frac{\partial {\cal L}}{\partial \theta} . \]

With the kinetic energy expressed via the angular displacement θ
\[ \mbox{K} = \frac{m}{2} \left( \dot{x}^2 + \dot{y}^2 \right) = \frac{m}{2}\, \ell^2 \left[ \left( \dot{\theta} \,\cos\theta \right)^2 + \left( \dot{\theta} \,\sin \theta \right)^2 \right] = \frac{m}{2}\,\ell^2 \dot{\theta}^2 . \]
and the potential energy
\[ \Pi = mg \left( \ell + y \right) = mg\ell \left( 1- \cos \theta \right) , \]
we derive the corresponding derivatives
\[ \frac{\partial \mbox{K}}{\partial \dot{\theta}} = I_0 \dot{\theta} = m\ell^2 \dot{\theta} , \qquad \frac{\partial \Pi}{\partial \theta} = mgy = mg\ell \, \sin \theta . \]
Using these expressions, we obtain from the Euler--Lagrange equation \( \displaystyle \frac{\text d}{{\text d} t} \left( \frac{\partial {\cal L}}{\partial \dot{\theta}} \right) = \frac{\partial {\cal L}}{\partial \theta} , \) for the Lagrangian \( {\cal L} = \mbox{K} - \Pi , \) the pendulum equation in a vacuum:
\[ \ddot{\theta} + \left( g/\ell \right) \sin \theta =0 \qquad \mbox{or} \qquad \ddot{\theta} + \omega_0^2 \sin \theta =0 \qquad \left( \omega_0^2 = g/\ell \right) , \]
where \( \ddot{\theta} = {\text d}^2 \theta / {\text d}t^2 \) , \( \omega_0 = \sqrt{g/\ell} >0 ,\) and g is gravitational acceleration.

 

      Now we consider a double pendulum, that consists of one pendulum attached to another. Consider a double bob pendulum with masses m1 and m2 attached by rigid massless wires of lengths ℓ1 and ℓ2. Further, let the angles the two wires make with the vertical be denoted θ1 and θ2, as illustrated in the figure to the left. Finally, let gravity acceleration be given by g. Then the positions of the bobs are given by
\begin{align*} x_1 &= \ell_1 \sin\theta_1 , \\ y_1 &= - \ell_1 \cos\theta_1 , \\ x_2 &= \ell_1 \sin\theta_1 + \ell_2 \sin\theta_2 , \\ y_2 &= - \ell_1 \cos\theta_1 - \ell_2 \cos\theta_2 . \end{align*}

The potential energy of the system is then given by
\begin{align*} \Pi &= m_1 g \left( \ell_1 + y_1 \right) + m_2 g \left( \ell_2 + y_2 \right) \\ &= m_1 g \ell_1 \left( 1 - \cos\theta_1 \right) + m_2 g \ell_2 \left( 1 - \cos\theta_2 \right) , \end{align*}
and the kinetic energy by
\begin{align*} \mbox{K} &= \frac{1}{2}\,m_1 v_1^2 + \frac{1}{2}\,m_2 v_2^2 \\ &= \frac{1}{2}\,m_1 \ell_1^2 \dot{\theta}^2_1 + \frac{1}{2}\,m_2 \left[ \ell_1^2 \dot{\theta}^2_1 + \ell_2^2 \dot{\theta}^2_2 + 2\ell_1 \ell_2 \dot{\theta}_1 \dot{\theta}_2 \cos \left( \theta_1 - \theta_2 \right) \right] . \end{align*}
The Lagrangian is then
\begin{align*} {\cal L} &= \mbox{K} - \Pi \\ &= \frac{1}{2} \left( m_1 + m_2 \right) \ell_1^2 \dot{\theta}_1^2 + \frac{1}{2} \,m_2 \ell_2^2 \dot{\theta}_2^2 + m_2 \ell_1 \ell_2 \dot{\theta}_1 \dot{\theta}_2 \cos \left( \theta_1 - \theta_2 \right) . \end{align*}
Performing calculations for one variable θ1, we obtain
\begin{align*} \frac{\partial {\cal L}}{\partial \dot{\theta}_1} &= m_1 \ell_1^2 \dot{\theta}_1 + m_2 \ell_2^2 \dot{\theta}_2 + m_2 \ell_1 \ell_2 \cos \left( \theta_1 - \theta_2 \right) , \\ \frac{\text d}{{\text d}t} \left( \frac{\partial {\cal L}}{\partial \dot{\theta}_1} \right) &= \left( m_1 + m_2 \right) \ell_1^2 \ddot{\theta}_1 + m_2 \ell_1 \ell_2 \ddot{\theta}_2 \cos \left( \theta_1 - \theta_2 \right) - m_2 \ell_1 \ell_2 \dot{\theta}_2 \sin \left( \theta_1 - \theta_2 \right) \left( \dot{\theta}_1 - \dot{\theta}_2 \right) , \\ \frac{\partial {\cal L}}{\partial \theta_1} &= - \ell_1 g \left( m_1 + m_2 \right) \sin\theta_1 - m_2 \ell_1 \ell_2 \dot{\theta}_1 \dot{\theta}_2 \sin \left( \theta_1 - \theta_2 \right) , \end{align*}
so the Euler-Lagrange differential equation \eqref{EqMotiv.2} for variable θ1 becomes
\[ \left( m_1 + m_2 \right) \ell_1^2 \ddot{\theta}_1 + m_2 \ell_1 \ell_2 \ddot{\theta}_2 \cos \left( \theta_1 - \theta_2 \right) + m_2 \ell_1 \ell_2 \dot{\theta}_2^2 \sin \left( \theta_1 - \theta_2 \right)+ \ell_1 g \left( m_1 + m_2 \right) \sin\theta_1 =0 . \]
Dividing through by ℓ1, this simplifies to
\begin{equation} \label{EqMotiv.4} \left( m_1 + m_2 \right) \ell_1 \ddot{\theta}_1 + m_2 \ell_2 \ddot{\theta}_2 \cos \left( \theta_1 - \theta_2 \right) + m_2 \ell_2 \dot{\theta}_2^2 \sin \left( \theta_1 - \theta_2 \right)+ g \left( m_1 + m_2 \right) \sin\theta_1 =0 . \end{equation}
      Similarly, for θ2:
\begin{align*} \frac{\partial {\cal L}}{\partial \dot{\theta}_2} &= m_2 \ell_2^2 \dot{\theta}_2 + m_2 \ell_1 \ell_2 \cos \left( \theta_1 - \theta_2 \right) , \\ \frac{\text d}{{\text d}t} \left( \frac{\partial {\cal L}}{\partial \dot{\theta}_2} \right) &= m_2 \ell_2^2 \ddot{\theta}_2 + m_2 \ell_1 \ell_2 \ddot{\theta}_2 \cos \left( \theta_1 - \theta_2 \right) - m_2 \ell_1 \ell_2 \dot{\theta}_1 \sin \left( \theta_1 - \theta_2 \right) \left( \dot{\theta}_1 - \dot{\theta}_2 \right) , \\ \frac{\partial {\cal L}}{\partial \theta_2} &= m_2 \ell_1 \ell_2 \dot{\theta}_1 \dot{\theta}_2 \sin \left( \theta_1 - \theta_2 \right) - \ell_2 m_2 g \sin\theta_2 , \end{align*}
so the Euler-Lagrange differential equation \eqref{EqMotiv.2} for variable θ2becomes

\[ m_2 \ell_2^2 \ddot{\theta}_2 + m_2 \ell_1 \ell_2 \ddot{\theta}_1 \cos \left( \theta_1 - \theta_2 \right) - m_2 \ell_1 \ell_2 \dot{\theta}_1^2 \sin \left( \theta_1 - \theta_2 \right) + \ell_2 m_2 g\,\sin \theta_2 =0 . \]
Dividing through by ℓ2, this simplifies to
\begin{equation} \label{EqMotiv.5} m_2 \ell_2 \ddot{\theta}_2 + m_2 \ell_1 \ddot{\theta}_1 \cos \left( \theta_1 - \theta_2 \right) - m_2 \ell_1 \dot{\theta}_1^2 \sin \left( \theta_1 - \theta_2 \right) + m_2 g\,\sin \theta_2 =0 . \end{equation}
Therefore, the double pendulum is described by the system of two ordinary differential equations of the second order \eqref{EqMotiv.4}, \eqref{EqMotiv.5}.    ■

Example 5: Consider two pendulums that are coupled by a Hookian spring with spring constant k (see figure below). Suppose that the spring is attached to each rod at a distance 𝑎 from their pivots, and that the pendula (plural from pendulum) are far apart so that the spring can be assumed to be horizontal during their oscillations. Let θ1, θ2 and ℓ1, ℓ2 be the angle of inclination of the shafts with respect to the downward vertical lines and lengths of the rods for the pendula, respectively. The kinetic energy is the sum of kinetic energies of two individual pendula:

\[ {\mbox K} = \frac{m_1}{2} \left( \ell_1 \dot{\theta}_1 \right)^2 + \frac{m_2}{2} \left( \ell_2 \dot{\theta}_2 \right)^2 . \]
The potential energy is accumulated by the spring, which amounts to \( \displaystyle \frac{k}{2} \left( a\,\sin\theta_1 - a\,\sin\theta_2 \right)^2 = \frac{a^2 k}{2} \left( \sin\theta_1 - \sin\theta_2 \right)^2 , \) and by lifting both bobs:
\[ \Pi = m_1 g \ell_1 \left( 1 - \cos\theta_1 \right) + m_2 g \ell_2 \left( 1 - \cos\theta_2 \right) + \frac{a^2 k}{2} \left( \sin\theta_1 - \sin\theta_2 \right)^2 . \]
     
polygon = Graphics[Polygon[{{-3, 0}, {3, 0}, {3, 1/4}, {-3, 1/4}}]];
line1 = Graphics[{Dashed, Line[{{2.5, 0}, {2.5, -5}}]}];
line2 = Graphics[{Dashed, Line[{{-2.5, 0}, {-2.5, -5}}]}];
line3 = Graphics[{Thickness[0.01], Line[{{2.5, 0}, {4.5, -3.6}}]}];
line4 = Graphics[{Thickness[0.01], Line[{{-2.5, 0}, {-3.9, -3.6}}]}];
circle1 = Graphics[Circle[{4.66, -4}, 0.38]];
circle2 = Graphics[Circle[{-4.0, -4}, 0.38]];
txt1 = Graphics[ Text[Style[Subscript[m, 1], FontSize -> 14, Purple], {-3.95, -4.0}]];
txt2 = Graphics[ Text[Style[Subscript[m, 2], FontSize -> 14, Purple], {4.68, -4.0}]];
spring = Graphics[{Thickness[0.01], Line[{{-3.35, -2.3}, {-1.7, -2.3}, {-1.5, -2.0}, {-1.5, -2.6}, \ {-1.0, -2.0}, {-1.0, -2.6}, {-0.5, -2}, {-0.5, -2.6}, {0, -2}, {0, \ -2.6}, {0.5, -2.0}, {0.5, -2.6}, {1, -2}, {1, -2.6}, {1.5, -2}, {1.5, \ -2.6}, {1.7, -2.3}, {3.69, -2.3}}]}];
txtk = Graphics[Text[Style["k", FontSize -> 18, Blue], {0.25, -1.4}]];
txt11 = Graphics[ Text[Style[Subscript[\[Theta], 1], FontSize -> 14, Purple], {-2.77, -1.8}]];
txt22 = Graphics[ Text[Style[Subscript[\[Theta], 2], FontSize -> 14, Purple], {2.98, -1.8}]];
Show[polygon, line1, line2, line3, line4, circle1, circle2, txt1, \ txt2, spring, txtk, txt11, txt22]
       Two pendula.            Mathematica code

Substituting these expressions into the Euler--Lagrange equations \eqref{EqMotiv.2}, we obtain the system of motion:
\[ \begin{split} m_1 \ell_1^2 \ddot{\theta}_1 + m_1 g\ell_1 \sin \theta_1 + a^2 k \left( \sin\theta_1 - \sin \theta_2 \right) \cos \theta_1 &= 0 , \\ m_2 \ell_2^2 \ddot{\theta}_2 + m_2 g\ell_2 \sin \theta_2 - a^2 k \left( \sin\theta_1 - \sin \theta_2 \right) \cos \theta_2 &= 0 \end{split} \]
   ■

Example 6: When a system with masses connected by springs is in motion, the springs are subject to both elongation and compression. It is clear from experience that there is some force, called the restoring force, which tends to return the attached mass to its equilibrium position. According to Hooke's law, the restoring force has direction opposite to the elongation (or compression) and is proportional to the distance of the mass from its equilibrium position. The corresponding constant of proportionality is called a spring constant.

We consider a situation when two bodies of masses m1 and m2 that are connected to three springs of negligible mass having spring constants k1 k2, and k3, respectively. We denote by x1 and x2 displacement of each body from its equilibrium position.
     
mass1 = Graphics[{Thickness[0.01], Line[{{-12.5, 0.5}, {-2.5, 0.5}, {-2.5, 4.5}, {-12.5, 4.5}, {-12.5, 0.5}}]}, Ticks -> None, Axes -> False];
mass2 = Graphics[{Thickness[0.01], Line[{{2.5, 0.5}, {2.5, 4.5}, {10.5, 4.5}, {10.5, 0.5}, {2.5, 0.5}}]}, Ticks -> None, Axes -> False];
spring1 = Graphics[{Thickness[0.005], Line[{{-12.5, 2.5}, {-14, 2.5}, {-14.4, 1.5}, {-14.8, 3.5}, {-15.2, 1.5}, {-15.6, 3.5}, {-16, 1.5}, {-16.4, 3.5}, {-16.6, 2.5}, {-17.5, 2.5}}]}, Ticks -> None, Axes -> False];
spring2 = Graphics[{Thickness[0.005], Line[{{-2.5, 2.5}, {-1.5, 2.5}, {-1, 1.5}, {-0.5, 3.5}, {0, 1.5}, {0.5, 3.5}, {1, 1.5}, {1.5, 3.5}, {2, 2.5}, {2.5, 2.5}}]}, Ticks -> None, Axes -> False];
spring3 = Graphics[{Thickness[0.005], Line[{{10.5, 2.5}, {11.5, 2.5}, {12, 1.5}, {12.5, 3.5}, {13, 1.5}, {13.5, 3.5}, {14, 1.5}, {14.5, 3.5}, {15, 2.5}, {16.5, 2.5}}]}, Ticks -> None, Axes -> False];
c1 = Graphics[{Thickness[0.01], Circle[{-10, 0}, 0.5]}];
c2 = Graphics[{Thickness[0.01], Circle[{-5, 0}, 0.5]}];
c3 = Graphics[{Thickness[0.01], Circle[{5, 0}, 0.5]}];
c4 = Graphics[{Thickness[0.01], Circle[{8, 0}, 0.5]}];
polygon = Graphics[{LightGray, Polygon[{{-17.5, -0.5}, {16.5, -0.5}, {16.5, -2.5}, {-17.5, -2.5}}]}];
poly1 = Graphics[{LightGray, Polygon[{{-17.5, -2.5}, {-18, -2.5}, {-18, 5}, {-17.5, 5}}]}];;
poly2 = Graphics[{LightGray, Polygon[{{16.5, -2.5}, {18, -2.5}, {18, 5}, {16.5, 5}}]}];
l1 = Graphics[{Dashed, Line[{{-12.5, 4.5}, {-12.5, 6.5}}]}];
l2 = Graphics[{Dashed, Line[{{2.5, 4.5}, {2.5, 6.5}}]}];
a1 = Graphics[Arrow[{{-12.5, 6}, {-8, 6}}]];
a2 = Graphics[Arrow[{{2.5, 6}, {7, 6}}]];
txt1 = Graphics[ Text[Style[Subscript[x, 1], FontSize -> 14, Purple], {-7, 6}]];
txt2 = Graphics[ Text[Style[Subscript[x, 2], FontSize -> 14, Purple], {8, 6}]];
t1 = Graphics[ Text[Style[Subscript[m, 1], FontSize -> 18, Purple], {-7, 2.5}]];
t2 = Graphics[ Text[Style[Subscript[m, 2], FontSize -> 18, Purple], {7, 2.5}]];
k1 = Graphics[ Text[Style[Subscript[k, 1], FontSize -> 14, Purple], {-15, 5}]];
k2 = Graphics[ Text[Style[Subscript[k, 2], FontSize -> 14, Purple], {0, 5}]];
k3 = Graphics[ Text[Style[Subscript[k, 3], FontSize -> 14, Purple], {13.5, 5}]];
Show[mass1, mass2, spring1, spring2, spring3, c1, c2, c3, c4, polygon, poly1, poly2, l1, l2, a1, a2, txt1, txt2,t1,t2,k1,k2,k3]
       Spring-mass system.            Mathematica code

This example can be used to model several mechanical systems because there are many situations when we observe masses connected with each other. For example, a multi-story building can be considered as masses (corresponding to each floor) connecting with springs. You can model an entire set of automobiles with several hundred masses connected to each other by several undred springs and can analyze how each part of the whole car vibrates when you drive this caravan along a bumpy road. Or you can model an entire set of rail cars with several masses, one for each car and one for the locomotive, all connected to each other by several springs and then analyze how each part of the whole train vibrates when the train chugs along climbing the mountain rail. You may think this kind of simple two mass and three spring model is not adequate for complicated models from real life, but in reality the logic and process of modeling is exactly the same. You would just have several dozens of differential equations instead of two or three equations, which is very similar to what you see here.

Let x1 and x2 be displacements of masses m1 and m2, respectively, from their equilibrium positions. So we consider these variables as canonical coordinates. Then the kinetic energy of these two objects will be

\[ \mbox{K} = \frac{1}{2}\,m_1 \dot{x}_1^2 + \frac{1}{2}\,m_2 \dot{x}_2^2 , \]
where \( \dot{x} = {\text d}x/{\text d}t \) is the velocity corresponding to the displacement x. The potential energy is concentrated in three springs and it is equal to the sum of potential energies of each spring. Correspondingly, the potential energy of each spring is proportional to the square of the amount the spring is stretched, so
\[ \Pi = \frac{1}{2}\,k_1 x_1^2 + \frac{1}{2}\,k_2 \left( x_1 - x_2 \right)^2 + \frac{1}{2}\,k_3 x_2^2 . \]
The Euler--Lagrange equations of the first kind for a system with two degrees of freedom are
\[ \frac{\text d}{{\text d}t} \,\frac{\partial {\cal L}}{\partial \dot{x}_i} - \frac{\partial {\cal L}}{\partial x_i} = 0, \qquad i=1,2. \]
When the kinetic energy does not depend on canonical coordinates (displacements in our case), and the potential energy does not depends on moments (velocities), the Euler--Lagrange equations become
\[ \frac{\text d}{{\text d}t} \,\frac{\partial {\text K}}{\partial \dot{x}_i} + \frac{\partial \Pi}{\partial x_i} = 0, \qquad i=1,2. \]
Since we have
\[ \frac{\partial {\cal L}}{\partial \dot{x}_i} = \frac{\partial {\text K}}{\partial \dot{x}_i} \qquad \Longrightarrow \qquad \frac{\text d}{{\text d}t} \,\frac{\partial {\text K}}{\partial \dot{x}_i} = m_i \ddot{x}_i , \]
and
\[ \frac{\partial \Pi}{\partial x_1} = k_1 x_1 + k_2 \left( x_1 - x_2 \right) , \qquad \frac{\partial \Pi}{\partial x_2} = - k_2 \left( x_1 - x_2 \right) + k_3 x_2 , \]
the Euler--Lagrange equations read as
\[ \begin{split} m_1 \ddot{x}_1 + k_1 x_1 + k_2 \left( x_1 - x_2 \right) &=0, \\ m_2 \ddot{x}_1 + k_3 x_2 - k_2 \left( x_1 - x_2 \right) &=0 . \end{split} \]
   ■

Example 7: We start with two loop circuit:
     
coil = ParametricPlot[{1*Cos[t*3 + Pi] + 1.5*t - 9, 2*Sin[t*3] - 8}, {t, 0, 5*Pi}, PlotLabel -> "L", Ticks -> None, Axes -> False, ImageSize -> Tiny, PlotStyle -> {Thickness[0.01]}];
resistor2 = Graphics[{Thickness[0.01], Line[{{-15, -8}, {-15, -2}, {-17, 0}, {-13, 0}, {-17, 2}, {-13, 2}, {-17, 4}, {-13, 4}, {-17, 6}, {-13, 6}, {-17, 8}, {-13, 8}, {-15, 9}, {-15, 11}}]}, PlotLabel -> Subscript[R, 1], Ticks -> None, Axes -> False];
resistor1 = Graphics[{Thickness[0.01], Line[{{15.5619, -8}, {18, -8}, {18, -2}, {20, 0}, {16, 0}, {20, 2}, {16, 2}, {20, 4}, {16, 4}, {20, 6}, {16, 6}, {20, 8}, {16, 8}, {18, 9}, {18, 11}}]}, PlotLabel -> Subscript[R, 1], Ticks -> None, Axes -> False];
c1 = Graphics[{Thickness[0.01], Circle[{-30, 2}, 2]}];
l1 = Graphics[{Thickness[0.008], Line[{{-25.4, 9.5}, {-25.4, 12.5}}]}];
l2 = Graphics[{Thickness[0.008], Line[{{-24.6, 9.5}, {-24.6, 12.5}}]}];
l3 = Graphics[{Thickness[0.01], Line[{{-10, -8}, {-30, -8}, {-30, 0}}]}];
l4 = Graphics[{Thickness[0.01], Line[{{-30, 4}, {-30, 11}, {-25.4, 11}}]}];
l5 = Graphics[{Thickness[0.01], Line[{{-24.6, 11}, {18, 11}}]}];
txtC = Graphics[ Text[Style["C", Bold, FontSize -> 14, Blue], {-25.5, 8.0}]];
txtV = Graphics[ Text[Style["V(t)", Bold, FontSize -> 14, Blue], {-25.5, 2.0}]];
txtL = Graphics[ Text[Style["L", Bold, FontSize -> 14, Blue], {2.0, -4.5}]];
txtR1 = Graphics[ Text[Style[Subscript[R, 1], Bold, FontSize -> 14, Blue], {-10.5, 2.5}]];
txtR2 = Graphics[ Text[Style[Subscript[R, 2], Bold, FontSize -> 14, Blue], {12.5, 2.5}]];
Show[l1, l2, l3, l4, l5, coil, resistor1, resistor2, c1, txtC, txtV, \ txtL, txtR1, txtR2]
       Two-loop circuit.            Mathematica code

We denote the current flowing through the left-hand loop by I1 and the current in the right-hand loop by I2 (both in the clockwise direction). The current through the resistor R1 is I1I2. The voltage drops on each element are

\begin{align*} \Delta V_{R_1} &= R_1 \left( I_1 - I_2 \right) \\ \Delta V_{R_2} &= R_2 I_2 , \\ \Delta V_{C} &= \frac{q_1}{C} , \\ \Delta V_{L} &= L\,\frac{{\text d}I_2}{{\text d}t} \end{align*}
where q1 is the electric charge in the left-hand side conductor C. The voltage analysis using Kirchhoff's law leads to the following equations:
\begin{align*} V(t) &= \frac{q_1}{C} + R_1 I_1 - R_1 I_2 , \\ 0&= R_2 I_2 + R_1 I_2 - R_1 I_1 + L\,\frac{{\text d}I_2}{{\text d}t} , \\ I_1 &= \frac{{\text d}q_1}{{\text d}t} . \end{align*}
So we obtain a system of differential equations with three unknows: q1, I1, and I2.

 

Now we consider another electric circuit with three loops.
     
connect[pointList_] := {Line[pointList], Map[Text[Style["", FontSize -> 18]] &, pointList[[{1, -1}]]]};
resistor[l_ : 1, n_ : 3] := Line[Table[{i l/(4 n), 1/3 Sin[i Pi/2]}, {i, 0, 4 n}]];
gap[l_ : 1] := Line[l {{{0, 0}, {1/3, 0}}, {{2/3, 0}, {1, 0}}}]
coil[l_ : 1, n_ : 3] := Module[{scale = l/(5/16 n + 1/2), pts = {{0, 0}, {0, 1}, {1/2, 1}, {1/2, 0}, {1/2, -1}, {5/16, -1}, {5/16, 0}}}, Append[Table[ BezierCurve[scale Map[{d 5/16, 0} + # &, pts]], {d, 0, n - 1}], BezierCurve[scale Map[{5/16 n, 0} + # &, pts[[1 ;; 4]]]]]];
capacitor[l_ : 1] := {gap[l], Line[l {{{1/3, -.5}, {1/3, .5}}, {{2/3, -.5}, {2/3, .5}}}]}
battery[l_ : 1] := {gap[ l], {Rectangle[l {1/3, -(2/3)}, l {1/3 + 1/9, 2/3}], Line[l {{2/3, -1}, {2/3, 1}}]}}
Options[display] = {Frame -> True, FrameTicks -> None, PlotRange -> All, GridLines -> Automatic, GridLinesStyle -> Directive[Orange, Dashed], AspectRatio -> Automatic};
display[d_, opts : OptionsPattern[]] := Graphics[Style[d, Thick], Join[FilterRules[{opts}, Options[Graphics]], Options[display]]]
at[position_, angle_ : 0][obj_] := GeometricTransformation[obj, Composition[TranslationTransform[position], RotationTransform[angle]]]
label[s_String, color_ : RGBColor[.3, .5, .8]] := Text@Style[s, FontColor -> color, FontFamily -> "Geneva", FontSize -> Large];
circuit = display[{connect[{{11.5, 0}, {0, 0}, {2.5, 2.5}}], coil[] // at[{2.5, 2.5}, Pi/4], connect[{{3.2, 3.2}, {5.75, 5.75}}], connect[{{5.75, 5.75}, {5.75, 4}}], resistor[] // at[{5.75, 3}, Pi/2], connect[{{5.75, 3}, {5.75, 2}}], connect[{{0, 0}, {2.5, .85}}], capacitor[] // at[{7.5, 1.4}, -Pi/10], connect[{{3.45, 1.2}, {5.75, 2}}], connect[{{5.75, 2}, {7.8, 1.3}}], resistor[] // at[{2.5, .85}, Pi/10], connect[{{8.15, 1.15}, {11.5, 0}}], connect[{{5.75, 5.75}, {8.5, 3}}], battery[] // at[{8.25, 3.25}, -Pi/4], connect[{{8.75, 2.75}, {11.5, 0}}] }]
r1 = Graphics[ Text[Style[Subscript[R, 1], Bold, FontSize -> 16, Black], {4.2, 0.7}]];
r2 = Graphics[ Text[Style[Subscript[R, 2], Bold, FontSize -> 16, Black], {4.8, 3.5}]];
txtC = Graphics[ Text[Style["C", Bold, FontSize -> 16, Black], {7.2, 0.8}]];
txtL = Graphics[ Text[Style["L", Bold, FontSize -> 16, Black], {1.5, 3.0}]];
Show[circuit, r1, r2, txtC, txtL]
       Three-loop circuit.            Mathematica code

We assume that the voltage source satisfies to US standards (Kurtus) with an amplitude of 127 volts and a frequency of 60 Hz to give the frequency-domain function: v(t) = 127 cos(120πt). For the resistors, capacitors, and inductors, I choose to keep the values in simple SI units of ohms, farads, and henries. For each resistor ​R₁ and ​R₂​, I used a magnitude of 10 ohms. For the capacitor and inductor, I used magnitudes of 5 farads and 5 henries, respectively.

To analyze the AC circuit, I utilized mesh analysis applied to the three loops of the circuit. Using the voltage equations for each of the three types of elements, I derived a system of three equations, but 4 unknowns.

\begin{align*} \left( I_1 - I_2 \right) R_1 + \left( I_1 - I_2 \right) R_2 + L\, \frac{{\text d}I_1}{{\text d}t} &= 0, \\ V_0 + V_c + \left( I_1 - I_2 \right) R_1 &= 0, \\ V_c + \left( I3 - I_1 \right) R_2 &= 0 , \\ I_3 - I_2 &= L\, \frac{{\text d}V_c}{{\text d}t} \end{align*}
The system contains four unknowns (Vc, I₁, I₂, I₃) and two differentials in the system and one for each capacitor and inductor.    ■

Example 8: According to statistical data, about 65% of marriages in the United States end in divorce within a 40-year period, resulting in huge social and economic consequences. The divorce rate for second marriages is even higher, about 75%. Professor John Gottman of the University of Washington identified five types of married couples based on their interactions.

To mathematically model the interplay of a married couple, let x(t) denote a measure of the husband's positivity (e.g., happiness) and let y(t) be the corresponding measurement for the wife's positivity (particular numerical values for such measurements can be found in John's work). In the absence of marital interaction, a single person tends to her/his own ``uninfluenced steady state:''     x0 and y0. This process is modeled by \[ \frac{{\text d}x}{{\text d}t} = h (x-x_0 ) , \qquad \frac{{\text d}y}{{\text d}t} = w (y-y_0 ) . \] After marriage, their interaction can be modeled as \begin{equation} \label{E511.m1} \frac{{\text d}x}{{\text d}t} = h (x-x_0 ) + I_1 (y) , \qquad \frac{{\text d}y}{{\text d}t} = w (y-y_0 ) + I_2 (x), \end{equation} where I1(y) is the influence exerted by the wife on the husband and I2(x) is the influence exerted by the husband on the wife. In a validating style of interaction, these functions can be modeled by \begin{equation} \label{E511.m2} I_k (z) = \begin{cases} a_k z , & \ \mbox{if } z>0 , \\ b_k z , & \ \mbox{if } z<0 \end{cases} \qquad (k=1,2). \end{equation} In a conflict avoiding style of interaction, the spouse who adopts this style avoids interacting with the other spouse through the negative range of his/her emotions. The corresponding function I(z) can be chosen as I(z) = H(z), where H(t) is the Heaviside function.    ■

Example 9: To attract your interest in the systems of differential equations, we present a simple model for the dynamics of love affairs (Strogatz 1988). A story written by William Shakespeare about 1594--96, describes a love affear between Romeo and Juliet.

Romeo is in love with Juliet, but in our version of this story, Juliet is a fickle lover. The more Romeo loves her, the more Juliet wants to run away and hide. But when Romeo gets discouraged and backs off, Juliet begins to find him strangely attractive. Romeo, on the other hand, tends to echo her: he warms up when she loves him, and grows cold when she disfavors him.

Let us introduce variables

\[ \begin{split} R(t) &= \mbox{ Romeo's love (or disfavor) for Juliet at time } t, \\ J(t) &= \mbox{ Juliet's love (or disfavor) for Romeo at time }t . \end{split} \]
Positive values of R, J signify love, negative values signify disfavor. Then a model for their star-crossed romance is
\[ \begin{split} \dot{R} &= a\,J , \\ \dot{J} &= -b\,R , \end{split} \]
where the parameters 𝑎 and b are positive, to be consistent with the story.
      Upon plotting the phase portrait, we see that the governing system has a center at (R, J) = (0,0).
VectorPlot[{2*y, -x}, {x, -1, 1}, {y, -1, 1}, ClippingStyle -> {Blue, Red}]
       Figure 1: Romeo and Juliet affair.            Mathematica code

   ■

Example 10: This example presents an important application of systems of differential equations to so-called mixed problems. Such a problem usually involves two or more tanks interconnected with pipes. The tanks contain a chemical dissolved in a fluid, which we take as salt (NaCl) and water, respectively. It is often happened that external substance enters one of the tanks and solution leaves the system at some rate. One simplifying assumption that we might make is that the concentration of the substance (salt) is kept uniform in te mixture. Then we can calculate the concentration of the substance in the mixture by dividing the instateneous amount of salt by the volume of the mixture in the tank at any time t. Multiplying this concentration by the exit rate of the mixture then gives the desired output rate of the substance. At a given time t, we are interested in the amount of salt in each tank.

Two tanks have 200 liters (L) and 300 liters of salt water, respectively. Tank A starts with 10 kg of salt. Tank B starts with 0 kg of salt, so it is full of pure water. 10% salt water enters tank A at 2 L/min. Solution exits the well-mixed tank A at 3 L/min and flows into tank B. Solution exits tank B and enters tank A at 1 L/min. Solution exits tank B and goes into the environment at 2.5 L/min. Therefore, the volume of A is unchanged, and the volume of B changes at a rate of -0.5 L/min.

We assume that the amount of water in the pipe between the tanks is negligible. Let x(t) be the amount of salt in kilograms (kg) in tank A at time t. Let y(t) be the amount of salt in kg in tank B at time t. The concentration in tank A is x/200, and the concentration of salt in tank B is y/(300 −0.5t). The equation describing the rate of change of the amount of salt in tank A is

\[ \frac{{\text d}x}{{\text d}t} = \begin{bmatrix} \mbox{inflow rate} \\ \mbox{of salt in A} \end{bmatrix} - \begin{bmatrix} \mbox{outflow rate} \\ \mbox{of salt from A} \end{bmatrix} = \begin{bmatrix} \mbox{inflow rate} \\ \mbox{of solution in A} \end{bmatrix} \times \begin{bmatrix} \mbox{concentration} \\ \mbox{of salt in water} \end{bmatrix} - \begin{bmatrix} \mbox{outflow rate} \\ \mbox{of solution in A} \end{bmatrix} \times \begin{bmatrix} \mbox{concentration} \\ \mbox{of salt in water} \end{bmatrix} . \]
Substituting the given data, we obtain
\[ \frac{{\text d}x}{{\text d}t} = 0.2 - 3\cdot \frac{x}{200} + 1\cdot \frac{y(t)}{300-0.5\,t} . \]
For tank B, we have
\[ \frac{{\text d}y}{{\text d}t} = \begin{bmatrix} \mbox{inflow rate} \\ \mbox{of salt from A} \end{bmatrix} - \begin{bmatrix} \mbox{outflow rate} \\ \mbox{of salt from B} \end{bmatrix} = 3 \left( \frac{x}{200} \right) - 3.5 \left( \frac{y}{300 -0.5\,t} \right) . \]
Combining these equations and using the given initial conditions, we obtain the initial value problem:
\[ \begin{split} \frac{{\text d}x}{{\text d}t} = 0.2 - 3\cdot \frac{x}{200} + 1\cdot \frac{y(t)}{300-0.5\,t} , \qquad x(0) = 10, \\ \frac{{\text d}y}{{\text d}t} = 3 \left( \frac{x}{200} \right) - 3.5 \left( \frac{y}{300 -0.5\,t} \right) , \qquad y(0) = 0. \end{split} \]
     
Graphics[{Red, Rectangle[{0, 0.275}, {0.2, 0.3}], Blue, Rectangle[{0.2, 0}, {0.7, 0.5}], Red, Rectangle[{0.7, 0.1}, {0.8, 0.125}], Red, Rectangle[{0.7, 0.15}, {0.8, 0.175}], Blue, Rectangle[{0.8, 0}, {1.3, 0.5}], Red, Rectangle[{1.3, 0.1}, {1.5, 0.125}], Green, Ellipsoid[{0.45, 0.52}, {0.25, 0.1}], Green, Ellipsoid[{1.05, 0.52}, {0.25, 0.1}], Blue, Ellipsoid[{1.05, 0}, {0.25, 0.1}], Blue, Ellipsoid[{0.45, 0}, {0.25, 0.1}], Yellow, Arrow[{{0, 0.29}, {0.3, 0.29}}], Yellow, Arrow[{{0.9, 0.11}, {0.6, 0.11}}], Yellow, Arrow[{{0.7, 0.16}, {0.9, 0.16}}], Yellow, Arrow[{{1.1, 0.11}, {1.3, 0.11}}]}]
       Figure 10: Two tanks connected.            Mathematica code

   ■

Example 11: Economists and financial analysts often utilize systems of ordinary differential equations (ODE’s) in order to model complex economic phenomena such as resource allocation and optimal growth. These systems often have no explicit analytical solution, but can be visualized and analyzed numerically in order to interpret long-term economic growth in terms of capital accumulation and consumption growth. In 1928, a British philosopher, mathematician, and economist Frank P. Ramsey (1903--1930) discovered what is regarded as the Ramsey growth model, explicating the optimal saving for society using differential equations. The model was expounded upon in 1965 by David Cass (1937--2008) and Tjalling Koopmans (1910--1985), forming the Ramsey--Cass--Koopmans (RCK) model that is currently used widely in macroeconomic theory in order to model optimal growth for an economy possessing labor augmenting technological progress.

We need to remind the main terminology.

  • Resource Allocation is the distribution of resources---usually financial---among competing groups of people or programs.
  • Optimal Growth: a firm's attempt to balance loss of current utility, as consumption is reduced to finance investment, against the future gain in utility, as the benefit from investment is realized.
  • Economic Growth: an increase in the production of economic goods and services, compared from one period of time to another.
  • Consumption Growth: an increase in the value of goods and services bought and consumed by people.
  • Capital Accumulation: an increase in assets from investments or profits.
  • RCK model is the canonical model of optimal growth for an economy with exogenous ‘labor-augmenting’ technological progress.
  • Labor Augmenting Technological Progress: technical progress that increases the effective labor input.
  • Labor Productivity: the hourly output of a country's economy.
  • Labor Supply: the amount of labor workers are willing to offer at various prices.
  • Depreciation: the cost of a tangible or physical asset over its useful life or life expectancy.
  • Elasticity; a measure of a variable's sensitivityto a change in another variable.
The core Ramsey-Cass-Koopmans (RCK for short) model is two-dimensional, comprising two coupled ODEs for per-capita wealth (k) and per-capita consumption (c).
\[ \frac{{\text d}k}{{\text d}t} = f(k) - c - k \left( \phi + \xi + \delta \right) , \qquad \frac{{\text d}c}{{\text d}t} = \frac{c}{\rho} \left( f' (k) - \theta -\xi \delta - \rho\phi \right) . \tag{11.1} \]
The terms in these equations are as follows:
  • f(k) = kα is a production function measuring the relative economic output in terms of k and a capital elasticity parameter α (the responsiveness of the output production to changes in the input capital).
  • φ is the growth rate of labor productivity (for example, due to technological innovation or efficiency improvements).
  • ξ is the growth rate of labor supply (for example, due to migration or population increase).
  • δ is the depreciation rate of capital (for example, due to inflation).
  • f'(k) = αkα−1 is the rate of change (derivative) of the production function f(k).
  • θ is an elasticity parameter indicating the tendency of consumers to smooth out their consumption over time.
  • ρ is the rate at which consumers discount their future consumption (for example, by indicating a preference for immediate consumption or attempting to preserve their long-term average consumption).

 

  1. Gourinchas, P.-O., The Ramsey-Cass-Koopmans Model, UC Berkeley, 2014.
  2. Simulating the Ramsey-Cass-Koopmans Model Using MATLAB and Simulink
  3. The Ramsey-Cass-Koopmans Model
  4. The Ramsey-Cass-Koopmans Model
  5. Sargent, Thomas J. & Stachurski, John, Cass-Koopmans Competitive Equilibrium with Python.
   ■

 

 

Return to Mathematica page
Return to the main page (APMA0340)
Return to the Part 1 Matrix Algebra
Return to the Part 2 Linear Systems of Ordinary Differential Equations
Return to the Part 3 Non-linear Systems of Ordinary Differential Equations
Return to the Part 4 Numerical Methods
Return to the Part 5 Fourier Series
Return to the Part 6 Partial Differential Equations
Return to the Part 7 Special Functions