A Promise kept

You might recall from my blog Expressions, Equations, and Formulae that my friend Solus “Sol” Simkin had asked me to write on differential equations as well, when he first requested a blog. I refused outright then because the subject—which has a history spanning more than 300 years—is both too vast and too complex for a blog at this level. But I later relented, and many were the occasions when I regretted my commitment to Sol.

As I gnawed on the problem over time, I decided on a modest presentation of differential equations with the following properties:

1. Simple definitions with links to more complete explanations;
2. A threefold approach to solving simple differential equations:
   1. Qualitatively using slope fields;
   2. Analytically (where possible); and
   3. Numerically using Euler’s method.
3. References will be given to well-known texts on differential equations to fill the large gaps left untouched;
4. After one anchor problem, dealt with thoroughly, I would bid my subject goodbye, hopefully leaving the reader with enough references for him or her not to feel stranded.

This blog is my promise kept to Sol.

A Gentle Entry

When you make a cup of tea, you must allow it to infuse. Likewise, with a cup of coffee, you need to let it percolate. With new knowledge, it is no different. You need time to allow the new knowledge to steep into your mind and find its place alongside old knowledge. It is only then that you become comfortable with what is new. Be gentle on yourself if things do not pop out of the page at once. You are making friends with mathematical ideas that took centuries to evolve and mature.

Differential equations model mechanical, electrical, chemical, biological, financial and other systems as they evolve in time or space. By solving them we can predict the future state of a system from its current state, provided that the model is adequate.

Variables, Functions, and Derivatives

Before we launch into what a differential equation is, we need to establish familiarity with certain basic concepts. In order not to leave anyone behind, I am proceeding slowly but surely below. This review is necessarily brief; if you need more details, please refer to the hyperlinks for each definition.

Variables are used to study relationships. An independent variable is a quantity that can take on different values, without restriction from other variables. Independent variables are often denoted by the symbols \(t\) for time and \(x\) for distance. A dependent variable is a quantity whose value is constrained by the independent variable through some function, and may be denoted by symbols such as \(y\), \(v\), \(F\).

Mathematically, a function \(f\) is a binary relation that uniquely assigns values in one set \(X\) to values in another set \(Y\), denoted by \(f: X \rightarrow Y\). Explicitly, this is written as \[ f: x \mapsto y \text{ where } x \in X \text{ and } y \in Y. \]

In the context of calculus and differential equations, we also write \(y(x)\) instead of \(y = f(x)\) to denote the functional dependence of \(y\) on \(x\). This notation is slightly less cumbersome when denoting the derivatives of \(y\) with respect to \(x\): \[ \begin{aligned} \frac{\mathrm{d}f(x)}{\mathrm{d}x} = \frac{\mathrm{d}y(x)}{\mathrm{d}x} = \frac{\mathrm{d}y}{\mathrm{d}x} &= y'(x)\quad\text{ for the first derivative}\\[1em] \frac{\mathrm{d^2}f(x)}{\mathrm{d}x^2} = \frac{\mathrm{d^2}y(x)}{\mathrm{d}x^2} = \frac{\mathrm{d^2}y}{\mathrm{d}x^2} &= y''(x)\quad\text{ for the second derivative}. \end{aligned} \] The first and second derivatives are often written using the convenient abbreviations \(y'(x)\) and \(y''(x)\) respectively.

What is a Differential Equation?

A differential equation (DE) connects a function with its derivative(s). In calculus, the function \(y(x)\) is known, and its derivative, \(y'(x)\), needs to be found. In differential equations, the derivative, \(y'(x)\), is known, and the function, \(y(x)\), needs to be identified. As with indefinite integration in calculus, the solution will throw up one or more arbitrary constants. The values of these constants must be determined by the initial conditions provided in the problem in order to obtain a unique solution.

A trivial differential equation

The simplest, most trivial differential equation is one which equates the first derivative of a function to zero: \[ y'(x) = 0. \qquad{(1)}\] If you think about the graph of the function, it will have to be parallel to the \(x\)-axis to exhibit zero gradient. But there are countless lines that may be drawn that are parallel to the \(x\)-axis. Which one is our solution? See Figure 1 which shows the graphs of four lines that each satisfy Equation 1, and there are countless more.

This is where we need an initial condition to force a unique solution. Will this condition relate to \(y'\)? Surely not, for the gradient has been defined in Equation 1 to be uniformly zero everywhere.

So, the constraint must be imposed by a value of \(y\) at some particular \(x\). In this case, the value of the particular \(x\) does not matter because \(y\) will be the same for all \(x\). It is conventional to give the constraint as the initial condition \(y(0) = C\) where \(C\) is some constant.

If we are told that \(y(0) = 5\), our solution to Equation 1 is simply \(y = 5\). The point to note here is that there is an infinite family of curves¹ which satisfy the differential equation, but it is the initial condition that selects one from this multiplicity.

The last step in solving any differential equation is to check the solution by substituting it back into the original differential equation. In this case, we have \(y = 5\), and we need to evaluate \(y'\) and verify it is zero. \[ y(x) = 5 \implies \frac{\mathrm{d}y}{\mathrm{d}x} = \frac{\mathrm{d}}{\mathrm{d}x}(5) = 0. \] And we are done.

A very simple second example

The second simple example is the DE \[ y' = 5. \qquad{(2)}\] Note Equation 2 is not very different from Equation 1 in that the right hand side (RHS) is still a constant, but non-zero this time. Well, what family of solutions does it lead to? The possible solutions are all straight lines with a slope of \(5\), as illustrated in Figure 2.

Again we depend upon an initial condition to compel a unique solution. Suppose we are told that \(y(0) = -5\). Then, out of the myriad straight lines with a slope of \(5\), we should select the one which satisfies the equation \(y = 5x - 5\), which is the lowest straight line in Figure 2.

As always, we check the solution \(y(x) = 5x - 5\). Then, \(y'(x) = \frac{\mathrm{d}y}{\mathrm{d}x} = \frac{\mathrm{d}}{\mathrm{d}x}(5x - 5) = 5\) which agrees with Equation 2. Are we done? Not yet. We still have to evaluate the initial condition which is \(y(0)\). Substituting for \(x = 0\) in \(y(x) = 5x - 5\), we get \(y(0) = 5(0) - 5 = -5\) which agrees with the initial condition \(y(0) = -5\). The solution has now been verified.

A not so trivial example

So far, we have got away with inspecting or guessing solutions to the DEs because they were so easy. Let us now tackle a DE which appears deceptively simple, but requires us to do some hard yards. This time, we will attach an initial condition and also explicitly show time \(t\) as the independent variable. The DE is: \[ \begin{aligned} y'(t) &= y(t);\\ y(0) &= 25. \end{aligned} \qquad{(3)}\] The first line of Equation 3 is interesting. It says \(\dfrac{\mathrm{d}y}{\mathrm{d}t} = y(t)\). The derivative equals the function itself. And we have encountered this very situation in the blog Eigenvalues and Eigenvectors—Why are they important?, where we identified the exponential function \(e^t\) as the eigenfunction of the differential operator \(\dfrac{\mathrm{d}}{\mathrm{d}t}\).

So, we may take the easy route and claim \(y(t) = e^t\), and throw in a constant \(C\) after that, to satisfy the initial condition. But where does that constant sit? Is it added to \(y(t)\) or does it multiply \(y(t)\)?² Let us look at both possibilities, though:

1. \(y(t) = e^t + C\): Together with \(y(0) = e^0 + C = 1 + C = 25 \implies C = 24\). Now, \(\dfrac{\mathrm{d}y}{\mathrm{d}t} = e^t\), but this does not equal \(e^t + 24\). So, adding \(C\) to \(e^t\) does not give us a solution.

2. \(y(t) = Ce^t\): Together with \(y(0) = Ce^0 = C(1) = C = 25 \implies C = 25\). Again, substituting into the original DE, we have \(y' = \dfrac{\mathrm{d}y}{\mathrm{d}t}(Ce^t) = C\dfrac{\mathrm{d}y}{\mathrm{d}t}(e^t) = Ce^t = y\). This solution satisfies the original DE when back-substituted, and the value \(C=25\) also satisfies the initial condition. So, \(y(t) = 25e^t\) is the correct solution.

But surely, mathematics should not be this slippery! And that observation is spot on. The established mathematical technique for solving such DEs is called separation of variables. If we have two variables, say, \(u\) and \(v\), we re-arrange terms by transposing them so that all terms involving \(u\) are on one side, and all terms involving \(v\) are on the other side. We then apply all available resources to compute what is being asked.

The notation used for calculus, while allowing for easy manipulation of symbols, is still a little difficult to justify from a philosophical viewpoint; some books have called it “informal” algebra, although each step may be rigorously justified. Each interpretation has its drawbacks. So, just learn the process, and apply it, if digging deeper troubles your heart or head. To know more, consult the listed references.

Separating variables, we proceed so: \[ \begin{aligned} \frac{\mathrm{d}y}{\mathrm{d}t} &= y; \text{ separate the variables}\\ \frac{\mathrm{d}y}{y} &= dt\\ \int\frac{1}{y}\mathrm{d}y &= \int\mathrm{d}t\\ \ln{\lvert y\rvert} &= t + k; \text{ $y \neq 0$ and $k$ constant; apply $\exp$ function to both sides}\\ \lvert y\rvert &= e^{t + k}\\ &= (e^t)(e^k); \text{ noting that $e^k = C$ is constant and positive}\\ y &= \pm Ce^t; \text{ allowing for $\lvert y \rvert$ with $C > 0$ to give new constant $C'$}\\ &= C'e^t. \end{aligned} \qquad{(4)}\] We can then apply the initial condition \(y(0) = 25 = C'e^0 = C'\) as before to get the solution \(y(t) = 25e^t\).

As before, the family of possible solutions is shown in Figure 3, where the solution we are after is represented by the topmost curve.

The Falling Skydiver Problem

The fourth and final problem we look at in this blog is also a simple one. It is called the “falling skydiver problem”. It has gained prominence in this century and is now a fixture in most textbooks on differential equations [2,3]. It has been treated in greater or lesser detail in most of the recommended texts. This is also my promised anchor problem and will therefore be looked at more thoroughly.

This problem forces us to integrate the real world of physics with the abstract world of mathematics to obtain a solution. It is simple enough to be rooted within our daily experience. Gravity and air resistance are no strangers to us. Who has not seen an object fall? And who has not experienced the rush of wind brush past us as travellers in a vehicle?

As we peel away layer by layer, onion-like, the factors impacting on the problem, we will also get a taste of what it means to solve a differential equation applied to the real world. Much ancillary material from related subjects must be absorbed and given its rightful place in the “drama of the solution” like actors in a stage play.

We start with a brief peek at the basic nomenclature on differential equations. The wise reader should consult the recommended texts. An algorithm to solve DE problems is presented next, followed by some definitions of variables and their SI units.

Nomenclature: ODE, PDE, Order, etc.

A differential equation involves derivatives and quantifies how certain dependent variables change with respect to independent variables. Differential equations may be ordinary (ODEs) where there is only one independent variable, or partial (PDEs) where there are two or more independent variables. The order of a differential equation equals the highest order of derivative appearing in the equation. Differential equations might also be linear, like \(y' + 2y +x = 0\), or nonlinear like \(y'=y^2\) or \(yy' = 2\). If the dependent and independent variables may be separated, the differential equation is called separable, and non-separable otherwise.

This is a blog: not a watered-down textbook. So, the definitions here are perforce brief. To get a firmer grounding in these terms, you should refer to one or more of the recommended texts [4–12]. The basic concepts and terminology are explained clearly and concisely in the text by Howell [12].

Algorithm to solve DE problems

1. When a quantitative problem is presented by Nature, the very first step to solve it is to define the independent variables—generally time \(t\) and space \(x, y, z\)—and the dependent variables.

2. The second step is to model the physical situation first using words to get a qualitative grasp, and then equations to quantify what we are seeking. This is abstraction. The model could be linear, nonlinear, dynamic, stochastic, etc., depending on the nature of the problem.

3. In the case of differential equations, the third step is to solve the equations graphically, analytically, numerically, or otherwise, to get a meaningful solution that can be used in practice.

Units in the real world

Differential equations arise naturally whenever mathematics is applied to problems that arise in nature. Therefore, one needs to become aware of the technical nomenclature of fields like physics, chemistry, biology, etc., and their units where applicable.

Because differential equations arise naturally in the solution of real world problems [13], they are usually entwined with the names and units for different quantities that must be defined and explained. In the context of mechanics let us address them sequentially, using SI units. This section is heavy on Physics.

1. Space: This is the three dimensional Euclidean space in which we live and through which matter moves. The length between two points in space is their distance. The set of points in a one-dimensional line is usually denoted by \(\mathbb{R}\). For two-dimensions, the set is the Cartesian product \(\mathbb{R} \times \mathbb{R} = \mathbb{R}^2\), and the points are denoted \((x, y)\) where each \(x \in \mathbb{R}\) and each \(y \in \mathbb{R}\). Likewise for three dimensional space where \((x, y, z) \in \mathbb{R}^3\).³

2. Time: This is the time we experience—24 hours or 86,400 seconds in a day—usually denoted by the letter \(t\).

3. Scalar: A quantity that only has magnitude. Examples include time, temperature, distance, mass, energy, etc.

4. Vector: A quantity that has both magnitude and direction. While distance is a scalar that defines the length between two points, displacement is a vector and encodes both the length between two points and their direction. Vectors are usually denoted in boldface, like \(\mathbf{a}\), or when written by hand, with a squiggle atop or underneath the symbol.

5. Distance: The difference in length between two points in space, without concern for their orientation.

6. Displacement: Displacement is the vector analogue of distance, and is often denoted by \(\mathbf{s}\). The displacement from \(A\) to \(B\) is the negative of the displacement from \(B\) to \(A\).

7. Speed: When a body moves a distance \(d\) metres in time \(t\) seconds, its average speed is \(\frac{d}{t}\) metres per second.

8. Velocity: When a body moves, velocity captures both its speed and the direction in which the change occurs. Velocity is the vector analogue to speed, and is usually denoted by \(\mathbf{v}\). A body whose speed is constant but whose direction is changing constantly, like a ball revolving on a string, experiences a changing velocity, and is hence undergoing acceleration.

9. Mass: Mass represents the resistance matter presents to motion. Roughly, on Earth, it is proportional to weight, which is a vector. Mass is generally denoted by \(m\). So, the heavier an object, the greater its mass.⁴ Mass multiplied by velocity equals momentum.

10. Momentum: Momentum, often denoted by \(\mathbf{p}\), is a vector and is the product of mass and velocity, \(\mathbf{p} = m\mathbf{v}\).

11. Acceleration: Acceleration, \(\mathbf{a}\) is defined as the time rate of change of the velocity vector. Using the symbology of calculus, \[ \mathbf{a} \triangleq \frac{\mathrm{d}\mathbf{v}}{\mathrm{d}t}. \] The symbol \(\triangleq\) is often used to denote a definition.

12. Force: Newton’s Second Law states that force is proportional to the rate of change of momentum it induces. Stated mathematically, bearing in mind that force is a vector, \[ \mathbf{F} \propto \frac{\mathrm{d}\mathbf{p}}{\mathrm{d}t} \triangleq k\frac{\mathrm{d}\mathbf{p}}{\mathrm{d}t}. \] The symbol \(\propto\) stands for is proportional to. To convert it to an equality, we must introduce a constant \(k\). By judiciously choosing units, it is possible to make the value of \(k\) equal to one. We then have the famous vector equation \[ \begin{aligned} \mathbf{F} &= \frac{\mathrm{d}\mathbf{p}}{\mathrm{d}t}\\ &= \frac{\mathrm{d}\left(m\mathbf{v}\right)}{\mathrm{d}t}; \text{ extracting the $m$ out as a constant}\\ &= m\frac{\mathrm{d}\mathbf{v}}{\mathrm{d}t}; \text{ acceleration is the derivative of velocity}\\ &= m\mathbf{a}. \end{aligned} \qquad{(5)}\] And there you have it: Newton’s second law stated mathematically in vector form as a differential equation in the second last line of Equation 5. Force is measured in newtons, mass in kilograms, and acceleration in metres per second per second, written \(\mathrm{m \; s^{-2}}\) or \(\mathrm{m/s^{2}}\). Note that the mass \(m\) increases with velocity at relativistic speeds, and is then not a constant but a variable \(m(v)\), and Equation 5 does not hold.

13. Newton’s Second Law: Perhaps the very first encounter with differential equations occurs for most people with the statement of Newton’s Second Law of Motion:

The net force on a body is proportional to the rate of change of the body’s momentum.

This is, of course, Equation 5 as we have just seen. Recall that \[ \mathbf{F} = m\frac{\mathrm{d}\mathbf{v}}{\mathrm{d}t} = m\mathbf{a} \] is a vector differential equation. The net force, \(\mathbf{F}\) and the acceleration \(\mathbf{a}\) are positive in the same direction and collinear. If you ever get confused, remember that we are dealing with vectors where both magnitudes and directions matter.

Gravity and the acceleration it induces

An object in free fall experiences the gravitational force pulling it toward the earth. To better understand this, we need a short detour on gravitation as enunciated by Newton.

The classical example of how objects fall from a height to the earth is the (apocryphal?) experiment of Galileo throwing rocks from atop the leaning tower of Pisa. If the two rocks have different weights but are otherwise similar, meaning air resistance may be neglected, they will strike the ground together.

We already know that the earth exerts a gravitational force on all earthbound objects. Newton proposed and concluded that, for an object on the surface of the Earth, the attraction of the earth on a body of mass \(m\) is \[ F = G\frac{mM}{R^2} \qquad{(6)}\] where \(G\) is the universal gravitational constant, \(m\) is the mass of the object of interest, \(M\) is the mass of the earth, and \(R\) is the radius of the earth, assuming a perfect sphere. The acceleration, \(g\), experienced by the mass \(m\) in free-fall is, from Equations 5, 6, the force of gravitational attraction divided by the mass of the object of interest. We may now write \[ \begin{aligned} \frac{F}{m} &= G\frac{M}{R^2}; \text{ noting that the LHS is $g$}\\ g &= G\frac{M}{R^2}. \end{aligned} \qquad{(7)}\] Let us substitute the values of \(G\), \(M\), and \(R\) in SI units to get a value for \(g\) the acceleration due to gravity felt by all objects on earth. \[ \begin{aligned} G &= 6.6743 × 10^{-11}\; \mathrm{m^{3} kg^{−1} s^{−2}}\\ M &= 5.9722 × 10^{24}\; \mathrm{kg}\\ R &= 6.3781 × 10^{6}\; \mathrm{m}; \text{ giving}\\ g &= G\frac{M}{R^2} = 9.7985 \approx 9.8\; \mathrm{m\,s^{-2}}. \end{aligned} \qquad{(8)}\] The value of \(g\) we have calculated is in the right ball park.

The force exerted by gravity on an object of mass \(m\), whether stationary or falling vertically, is its weight, \(W = mg\; \mathrm{N}\), where the \(\mathrm{N}\) stands for newtons, the SI unit of force.

Air resistance

Air resistance is usually ignored in problems meant for elementary school. However, in the case of the parachuting skydiver, it is air resistance and the drag it produces that protects the skydiver from a serious accident or worse.

Air resistance is a deceptively complex phenomenon that is usually discounted in elementary treatments of free fall in Physics. One way to appreciate it is to think of swimming in water. Water is denser than air and we feel its resistance whether we are wading through it, swimming in it, or diving through it. Since air is a fluid like water, it too presents resistance just like water, and that resistance cannot be ignored at high speeds.

Three points make for tractable analysis of air resistance [12]:

1. It is assumed that there are no air currents or draughts; think of water currents at sea to understand why.

2. For an object that does not move, the air resistance is zero. It increases with the speed of the moving body. This is not an assumption but a fact.

3. Like friction, air resistance always opposes motion. In one dimension, the sign of the force of air resistance and the sign of the velocity are opposite. This is a fact and a convention. A skyward air resistance and an earthward velocity have opposite signs.

The magnitude of the drag is usually expressed as [14] \[ R = \frac{1}{2}C\rho\ Av^{2}. \qquad{(9)}\] where \(C\) is the drag coefficient, \(\rho\) is the density of air, \(A\) is the area presented by the falling body to the air, and \(v\) is the magnitude of the velocity of the falling object.

What are the equations governing air resistance? The answers are not simple.

1. At low velocities, air resistance is proportional to the velocity of the object [6].
2. At higher velocities, the drag is proportional to the square of the velocity [5].

For this problem, we will use [14]: \[ \begin{aligned} R &= \frac{1}{2}C\rho\ Av^{2}\\ &= kv^{2}\; \text{ where}\\ k &= \frac{1}{2}C\rho\ A \end{aligned} \qquad{(10)}\]

Let us now compute the value of \(k\). We have \[ \begin{aligned} \rho &= 1.21\, \mathrm{kg}\,\mathrm{m}^{-3}\; & \text{ density of air}\\ C &= 0.70\; & \text{ drag coefficient}\\ A &= 0.18\, \mathrm{m}^2\; & \text{ area presented by the body to air}\\ k &= \frac{1}{2}C\rho\ A\; \text{ gives}\\ &= 0.076 23\, \mathrm{N\, m^{-2}\, s^{2}}\\ &= 0.076 23\, \mathrm{kg\,m^{-1}}. \end{aligned} \qquad{(11)}\] We will revert to Equation 11 when we solve the DE stepwise later in this blog, and will need these numbers then.

Modelling the fall

Let us backtrack a bit here. The skydiver starts from a helicopter or aeroplane with zero downward velocity. At that point, the chute is unopened, and the skydiver falls more or less like a rock hurtling to the earth, opposed by air resistance. By assuming a horizontal position, the skydiver presents more surface area to the air, and thereby experiences increased drag, which slows his or her fall.

Let us assume that positive velocity \(\mathbf{v}\) is toward the earth and denote its magnitude by \(v\). The air resistance varies with \(v\), but is negative because it is in the direction opposite to the downward vector \(\mathbf{v}\).

As the skydiver’s earthward speed increases, the equations governing the drag force change from being proportional to \(v\) to being proportional to \(v^2\) [14].

The parachute does not open the moment the skydiver falls from the aircraft. Only after some time, when the rip cord is pulled, does the parachute unfurl and engage. Once the chute opens, the skydiver experiences an abrupt upward pull. See this YouTube video [15] where this is apparent. After this initial skyward yank, there are two forces: gravity and air resistance. As the skydiver’s speed increases, drag is assumed to be proportional to \(v^2\).

The four stages of the fall

To recapitulate, there are four stages in the fall:

1. Start of the jump: The skydiver starts with zero downward velocity and zero upward air resistance. Both of these values increase with time.

2. Parachute just opens: The downward falling skydiver is pulled up abruptly as the chute opens [15], but begins descent again shortly afterward. The drag is better modelled by \(R = kv\) [14].

3. Parachute has fully opened: As the velocity increases, the retarding force due to air-resistance is more accurately modelled by \(R = kv^2\).

4. Terminal descent: There will be a velocity at which the weight of the skydiver equals the drag. At this point of force balance, we move from Newton’s second law to the first law. The skydiver experiences uniform motion in a straight line at a constant speed called the terminal velocity [14].

Numerical values

When modelling the skydiver’s fall and estimating speeds and forces, symbols alone will not give us a feel for what is happening. To this end. I have lifted some numbers from a Physics text [14] to inject realism into our analysis.

1. The acceleration due to gravity \(g\) is \(9.8\; \mathrm{m}\mathrm{s}^{-2}\), as derived above Equation 8;

2. The mass \(m\) of the skydiver is assumed to be \(70\; \mathrm{kg}\);

3. In the air-resistance equation \(kv^2\), the coefficient \(k\) is, from Equation 11, taken to be: \[ k = 0.07\; \mathrm{N} \mathrm{s}{^2 }\mathrm{m}^{-2}. \]

Falling skydiver: problem statement

The problem stated here has been adapted from Blanchard, Devaney, and Hall [5].

The velocity \(v\) of a freefalling skydiver is well modeled by the differential equation \[ m\frac{\mathrm{d}v}{\mathrm{d}t} = mg - kv^2 \qquad{(12)}\] where \(m\) is the mass of the skydiver, \(g\) is the gravitational constant, and \(k\) is the drag coefficient determined by the position of the diver during the dive. The constants \(m\), \(g\), and \(k\) are positive.

Assume that \(v = 0\) at \(t = 0\), i.e., \(v(0) = 0\).

1. Perform a qualitative analysis of this model.

2. Calculate the terminal velocity of the skydiver. Express your answer in terms of \(m\), \(g\), and \(k\).

3. Plot a slope field curve for Equation 12. What can you say about the value computed for (b) above?

4. Find the general solution of the differential equation Equation 12.

5. Confirm your answer to (b) above by calculating the limit of \(v(t) \text{ as } t\to\infty\).

6. Solve Equation 12 numerically using software of your choice.

Falling skydiver: solution

For convenience, each part of the question is answered under a separate heading below.

We have been told to assume that when \(t\) is zero, \(v\) is also zero. Moreover, the differential equation stated in the problem obviously applies throughout. Therefore, while we have identified four stages in the skydiver’s fall, we need to deceive ourselves into believing that there are only two stages:

1. fall with air resistance proportional to \(v^2\), and

2. fall at constant, terminal velocity.

This is a convenient fiction, but is necessary to make problems mathematically solvable. After the solution, errors may be analyzed to determine whether the simplification introduced was justified or whether the problem needs to be re-worked more accurately. If our focus is on the terminal velocity alone, the simplified assumptions we have introduced are justified and sufficient, because they are valid for the last stage of the fall.

Qualitative analysis

Qualitative analysis means different things to different people. I will answer it below as I understand it.

1. In a free body diagram, we abstract the object of interest and replace all other objects by the forces they exert on the abstracted body.

I wanted to inject some realism, eye candy, and levity into the problem, and so have called it an “almost free body diagram” because other objects are also shown. But the two forces acting on the skydiver are clearly depicted.

Let the downward direction denote positive velocity \(v\). Then, the net downward force \(F\) is the difference between the downward acting weight \(W = mg\) and the upward acting drag \(R = kv^2\): \[ F = m\frac{\mathrm{d}v}{\mathrm{d}t} = W - R = mg - kv^2. \qquad{(13)}\]

The qualitative analysis goes thus. The skydiver starts with zero velocity in the vertical direction. Therefore, the initial drag is also zero.

Initially, the skydiver accelerates like any downward falling object. Because \(g\) is constant, the downward force is constant, but the velocity increases linearly with time. But the drag kicks in as the skydiver falls and his/her velocity becomes non-zero. The drag increases quadratically with time. So, at some point, the drag will equal the gravitational pull of the earth. At that point in time the two forces—weight and drag—will be equal and Newton’s first Law will apply. The skydiver will fall to the earth at constant speed. This we will call the terminal velocity.

Terminal velocity

2. The downward force of gravity is constant at \(mg\). The upward force of air resistance is quadratic with velocity. Therefore, a plot of these two forces with time will look like the graphs shown in Figure 5. Where the two curves intersect, the net or residual force on the skydiver will be zero. Newton’s First Law will kick in, and the skydiver will fall to earth at a final, constant velocity called the terminal velocity.

At the terminal velocity, \(v_t\), we have, substituting approximations to the values for \(m\), \(g\), and \(k\) assigned above: \[ \begin{aligned} mg &= kv_t^2\; \text{ transposing}\\ kv_{t}^2 &= mg\\ v_t^{2} &= \frac{mg}{k}\\ v_t &= \sqrt{\frac{mg}{k}}\\ &\approx\sqrt{\frac{70\times 10}{0.07}}\\ &\approx 100\;\mathrm{m\,s^{-1}}. \end{aligned} \qquad{(14)}\]

Slope field plot

3. The original DE is \[ m\frac{\mathrm{d}v}{\mathrm{d}t} = mg - kv^2; \; \text{$m = 70, g =10, k = 0.07$} \qquad{(15)}\] and this may be divided through by \(70\) to give the identical DE \[ \frac{\mathrm{d}v}{\mathrm{d}t} = 10 - 0.001v^2. \qquad{(16)}\] This is a first-order, nonlinear, ordinary differential equation, that belongs to a class called Riccati’s equation. We may then plot, from its very definition, the slope of the solution curve at different points. We will get a series of arrows that give a piecewise picture of the solution. The resulting plot is called a slope field, and is shown for our case in Figure 6. The arrows are changing direction—almost ninety degrees—to align with the \(t\) axis around the value \(v = 100\) which is the terminal velocity. Note that \[ \left.\frac{\mathrm{d}v}{\mathrm{d}t}\right\rvert_{v=100} = 10 - 0.001(100^2) = 0. \qquad{(17)}\]

To preserve similarity of look and feel, I have included a font initialization script font_select.py which must also be downloaded in order for the other scripts to work. The Python 3 script skydiver_slope_field.py was used to produce Figure 6.

General solution to the DE

4. The general solution to Equation 12 is the one with an arbitrary constant. For this part, we will retain the symbols in Equation 12 without substituting numerical values, as was done for part (c) above. We start by separating the variables: \[ \begin{aligned} m\frac{\mathrm{d}v}{\mathrm{d}t} &= mg -kv^2;\; \text{separate variables by dividing by $mg - kv^2$}\\ \frac{m}{mg-kv^2} \mathrm{d}v &= \mathrm{d}t; \; \text{integrate left side wrt $v$ and right side wrt $t$}\\ m\int \frac{1}{mg-kv^2}\,\mathrm{d}v &= \int \mathrm{d}t. \end{aligned} \qquad{(18)}\]

We want the denominator of the integrand in the LHS to be a difference of two squares to get a standard integral of a rational function. For this to happen, we divide through by \(k\) to require \(\left[\dfrac{mg}{k} - v^2\right]\) to be perfect square. This is easily achieved by setting the constant \(\dfrac{mg}{k}\) to be \(a^2\), i.e., by substituting \(a = \sqrt{\dfrac{mg}{k}}\) into the LHS of Equation 18.

We then have \[ \begin{aligned} m\int\frac{1}{mg-kv^2}\,\mathrm{d}v &= \int\mathrm{d}t\\ m\int\frac{1}{k\left[\frac{mg}{k} - v^2\right]}\,\mathrm{d}v &= \int\mathrm{d}t;\; \text{substitute $a = \sqrt{\dfrac{mg}{k}}$ into LHS}\\ \frac{m}{k}\int\frac{1}{a^2 - v^2}\,\mathrm{d}v &= \int\mathrm{d}t. \end{aligned} \qquad{(19)}\]

The indefinite integral on the LHS of Equation 19 is a standard integral. There is no special merit in working out standard integrals from first principles. It used to be that such standard integrals were looked up in books that classified and listed integrals. But, in this day and age, books like Dwight’s Table of Integrals and Other Mathematical Data [16] have faded into the sunset, like the slide rule. They have been replaced by online resources like Wolfram Alpha [17] and SymPy [18].

Resuming our solution, \[ \begin{aligned} \frac{m}{k}\int\frac{1}{a^2 - v^2}\,\mathrm{d}v &= \int\mathrm{d}t; \; \text{ the LHS is a standard integral}\\ \left[\frac{m}{k}\right]\left[\frac{1}{a}\right] \tanh^{-1}\frac{v}{a} &= t + C; \; \text{ substituting for $a = \sqrt{\dfrac{mg}{k}}$}\\[0.75em] \frac{m}{\sqrt{mgk}}\tanh^{-1}\frac{v}{\sqrt\frac{mg}{k}} &= t + C\\ \sqrt{\frac{m}{gk}}\tanh^{-1}\frac{v}{\sqrt\frac{mg}{k}} &= t + C\\ \end{aligned} \qquad{(20)}\] We require \(v\) to be the subject of the formula. So, we need to isolate \(v\) on the LHS and move everything else to the RHS: \[ \begin{aligned} \sqrt{\frac{m}{gk}}\tanh^{-1}\frac{v}{\sqrt\frac{mg}{k}} &= t + C;\; \text{ isolate $v$ alone on the LHS}\\ \tanh^{-1}\frac{v}{\sqrt\frac{mg}{k}} &= \left[\sqrt{\frac{gk}{m}}\right][t +C]; \; \text{ apply the $\tanh$ function to both sides}\\ v(t) &= \sqrt\frac{mg}{k}\tanh{\left[\sqrt{\frac{gk}{m}}(t + C)\right]}. \end{aligned} \qquad{(21)}\] And there we have it: the general solution to the DE in the last line of Equation 21.

The hyperbolic tangent function, \(\tanh\), is defined as: \[ \tanh x = \frac{\sinh x}{\cosh x} = \frac {e^x - e^{-x}} {e^x + e^{-x}} = \frac{e^{2x} - 1} {e^{2x} + 1}, \] and is illustrated in Figure 7.

The particular or exact solution

We have been given the initial condition \(v(0) = 0\). Using this, we get the particular solution thus: \[ \begin{aligned} v(0) &= \sqrt\frac{mg}{k}\tanh{\left[\sqrt{\frac{gk}{m}}(0 + C)\right]};\; \text{ i.e., }\\ 0 &= \sqrt\frac{mg}{k}\tanh{\left[C\sqrt{\frac{gk}{m}}\right]}\\ \end{aligned} \qquad{(22)}\] In Equation 22, we know that \(\sqrt\frac{mg}{k} \neq 0\) which means that \(\tanh{\left[C\sqrt{\frac{gk}{m}}\right]} = 0\). To solve this systematically, we need to know how the \(\tanh\) function behaves; this is shown in Figure 7 above.

It is clear that the only value at which \(\tanh(x) = 0\) is at \(x = 0\). This implies \(C=0\) because \(\sqrt{\frac{gk}{m}} \neq 0\). Therefore, our particular solution is \[ \begin{aligned} v(t) &= \sqrt\frac{mg}{k}\tanh\left[\sqrt{\frac{gk}{m}}t\right].\\ \end{aligned} \qquad{(23)}\] Substituting the values \(g \approx 10\), \(m = 70\) and \(k \approx 0.07\) into Equation 23, we get Equation 16 as expected: \[ \begin{aligned} m\frac{\mathrm{d}v}{\mathrm{d}t} &= mg - kv^2;\; \text{divide through by $m$}\\ \frac{\mathrm{d}v}{\mathrm{d}t} &= g - \frac{k}{m}v^2;\; \text{ substitute values for $m$, $g$, and $k$}\\ &= 10 - 0.001v^2, \end{aligned} \] whose exact solution is: \[ v(t) = \sqrt\frac{mg}{k}\tanh\left[\sqrt{\frac{gk}{m}}t\right] \qquad{(24)}\] which is numerically \[ v(t) = 100 \tanh(0.1t) \qquad{(25)}\]

The partial fraction route

What if we had taken the partial fraction route to the solution? Would we have gotten the same answer? Let us proceed step-by-step.

1. The starting point is \(\dfrac{dv}{dt} = 10 - 0.001v^2\) with \(v(0) = 0\).

2. Separating variables, \(\dfrac{dv}{10 - 0.001v^2} = dt\).

3. Integrating both sides, \[ \int \frac{dv}{10 - 0.001v^2} = \int dt. \]

4. Factoring out 10 from the denominator, \[ \begin{aligned} \int \frac{dv}{10(1 - 0.0001v^2)} &= \int dt\\ \frac{1}{10} \int \frac{dv}{1 - (0.01v)^2} &= \int dt. \end{aligned} \]

5. Let \(u = 0.01v\), so \(du = 0.01dv\), and \(dv = 100du\): \[\begin{aligned} \frac{1}{10} \int \frac{100du}{1 - u^2} &= \int dt\\ 10 \int \frac{du}{1 - u^2} &= \int dt. \end{aligned} \]

6. Using partial fractions: \[ \begin{aligned} \frac{1}{1 - u^2} &= \frac{1}{(1 - u)(1 + u)}\\ &= \frac{A}{1 - u} + \frac{B}{1 + u};\; \text{multiply both sides by $(1-u^2)$}\\ 1 &= A(1 + u) + B(1 - u); \; \text{ equating constants, we get}\\ 1 &= A + B; \; \text{and equating coefficients of $u$, we get}\\ 0 &= A - B\\ \end{aligned} \] This leads to \(A = B\) and \(A + B = 1\) giving \(A = B = \dfrac{1}{2}\).

7. Substituting, \[ \begin{aligned} 10 \int \left( \frac{1/2}{1-u} + \frac{1/2}{1+u} \right) du &= \int dt\\ 5 \int \left( \frac{1}{1-u} + \frac{1}{1+u} \right) du &= \int dt\\ 5 (-\ln|1-u| + \ln|1+u|) &= t + C\\ 5 \ln \left| \frac{1+u}{1-u} \right| &= t + C\\ \end{aligned} \]

8. Substituting back \(u = 0.01v\): \[ \begin{aligned} 5 \ln \left| \frac{1+0.01v}{1-0.01v} \right| &= t + C; \; \text{applying the initial condition $v(0) = 0$}\\ 5 \ln \left| \frac{1+0}{1-0} \right| &= 0 + C\\ 5 \ln(1) &= C\\ 0 &= C. \end{aligned} \]

9. So, \[ \begin{aligned} 5 \ln \left| \frac{1+0.01v}{1-0.01v} \right| &= t\\ \ln \left| \frac{1+0.01v}{1-0.01v} \right| &= \frac{t}{5}\\ \left| \frac{1+0.01v}{1-0.01v} \right| &= e^{t/5} \end{aligned} \]

10. Since \(v(0) = 0\) and \(v\) is increasing, because \(v' > 0\) initially, \(\frac{1+0.01v}{1-0.01v}\) is positive, so we can drop the absolute value: \[ \begin{aligned} \frac{1+0.01v}{1-0.01v} &= e^{t/5}\\ 1 + 0.01v = e^{t/5}(1 - 0.01v)\\ 1 + 0.01v &= e^{t/5} - 0.01v e^{t/5}\\ 0.01v + 0.01v e^{t/5} &= e^{t/5} - 1\\ 0.01v(1 + e^{t/5}) &= e^{t/5} - 1\\ v &= 100 \frac{e^{t/5} - 1}{e^{t/5} + 1}\\ v &= 100 \frac{e^{t/10} - e^{-t/10}}{e^{t/10} + e^{-t/10}}\\ v &= 100 \tanh(t/10)\\ \end{aligned} \]

11. Therefore, the exact solution is: \[ v(t) = 100\tanh(0.1)t, \; \text {as before}. \]

Verifying the solution by differentiating

Earlier in this blog I had suggested that the last step in the solution of any differential equation should be to verify the answer by differentiating it to determine if we indeed get back the original DE. If we do, we have confirmed the correctness of our result.

From Equation 25, we know that \[ \begin{aligned} v(t)&= 100\tanh(0.1t), \text{ or, transposing terms, that}\\ \tanh(0.1t) &= \frac{v(t)}{100}. \end{aligned} \] We will need this last equation to back substitute below for \(\tanh\) to get back the original DE in terms of \(v(t)\) alone, without involving the independent variable \(t\). It might seem like sleight of hand, and I am alerting you to it now so that it seems less so later.

Recall that the derivative of \(\tanh\) is \(\sech^2\), and moreover that \(\tanh^2 t + \sech^2 t = 1\). \[ \begin{aligned} \text{Since}\; v(t)&= 100\tanh(0.1t)\\ \frac{\mathrm{d}v}{\mathrm{d}t} &= 100(0.1)\sech^2 (0.1t)\\ &= 10(1 - \tanh^2(0.1t))\\ &= 10 - 10\left[\frac{v(t)}{100}\right]^2\\ &= 10 - 0.001v^2(t),\; \text{which is the original DE, and we are done.}\\ \end{aligned} \]

Verification of terminal velocity

5. The graph of the hyperbolic tangent \(\tanh(t)\) approaches the asymptote \(1\) as \(t\) approaches \(\infty\), i.e., \(\lim_{t\to\infty}\tanh(t) = 1\), as shown in Figure 7. Therefore, we may set the term inside the square brackets in Equation 21 to \(1\) and then get the terminal velocity to be \(\sqrt\frac{mg}{k}\) which accords with the answer to (b). Pedantically, \[ \lim_{t\to\infty}\sqrt\frac{mg}{k}\tanh{\left[\sqrt{\frac{gk}{m}}t\right]} = \sqrt\frac{mg}{k}. \]

For a visual appreciation of the interplay between the slope field, the analytical solution, and the asymptote, see Figure 8.

The Python 3 script slope_field_analytical_soln.py was used to generate Figure 8.

Euler’s method

Euler’s method is a numerical method for estimating the solution to a first order DE using the value of a starting point and its gradient at that point. It is a segmental piecewise linear approximation to the actual solution and its accuracy varies with the step size used. The Runge-Kutta methods include and extend the Euler method in the quest for greater accuracy.

The interested reader is referred to a detailed online explanation of the Euler method [19,20]. Typically this part of the problem will be tackled in a computing laboratory session which will allow for exploration, experimentation, and experience.

In Figure 9 above, we have four plots, three of which are identifiable from the legend. The fourth, in grey, is the slope field shown as arrows. Note how closely the Euler approximation tracks the arrows and the exact solution.

In an earlier iteration, I had failed to use floating point numbers for the variables in my code, and ended up with errors that accumulated and showed significant differences between the exact solution and the Euler approximation. Once the step size was changed from \(1\) to \(1.0\), the two curves effectively coincided and were indistinguishable. Awareness of nuances like this plays an important role in correctly using numerical methods. One should develop the knack of unearthing plausible reasons for inaccurate results, and know what remedial steps to undertake. The Python 3 script for the above plot is downloadable as eulers_method.py.

Beware the uncritical use of AI

I have relied heavily on software for generating the plots for this blog. My mainstay has usually been the TikZ-PGF suite and PGFplots in particular.

However when it came to plotting slope fields, exact solutions, and Euler approximations, I had to look further afield for software assistance. I have been using Julia for most of the programs on these blogs, but I think Julia is not ready for prime time when it comes to plotting.

Accordingly, I have used the more mature matplotlib, which is written in Python, for the later graphs.

But because I am not familiar with Python, I specified clear requirements to AI resources like ChatGPT, Copilot, Gemini and Grok for them to code. The programs they provided in response were all respectable, well commented, and generally complied with standards. They also usually compiled without error and ran as expected. One advantage of this approach was that I could check both the programs and the outputs from one AI agent against those from another.

But imagine my surprise when one of the AI engines gave a curve that looked a “bit off”. This was ultimately tracked down to a decimal point error that was made by the AI agent—this was not the \(1\) versus \(1.0\) rounding error noted above. So, be aware that AI too is as fallible as human beings, because it is trained on data generated by human beings. Do check whatever result you get from one AI source, with results from at least two other AI sources.⁶ As Ronald Regan famously said, “Trust but Verify”.

Closing Thoughts

Differential equations link the abstract world of mathematics to the three-dimensional world of space and the one-dimensional world of time. Linear algebra and differential equations are probably the two most applied areas of mathematics. Accordingly they are taught to students of engineering, the physical sciences, the biological sciences, economics, etc. A good grounding in these two areas bodes well for those planning on making a career in these fields.

The variety of ways in which a DE may be viewed—qualitatively, graphically, analytically, and numerically—confer a certain roundedness and maturity on the student of the subject. Moreover since particular solutions depend on initial or boundary conditions, the elementary school mathematical fixation on a single correct answer is gradually dispelled. Mathematics is seen more as a philosophy than as a mere calculating machine.

The ready availability of computing resources takes away much of the grunt work involving differential equations. But the student needs judgement, fostered by familiarity with the subject, to use these aids wisely. And lots of exposure and practice. No computer can stand in for that sort of mastery.

The advent of artificial intelligence (AI) might cause stirrings in your mind that we could ask a chatbot to solve the DE. These bots are not infallible, at least at present. I have encountered enough instances of errors—both trifling and catastrophic—that convince me of the need for human shepherding. Do not let some so-called smart bot take you for a ride. Do remain in control.

Recommended Books

To many students, mathematics is daunting, calculus even more so, and differential equations are downright terrifying. Nevertheless, there are authors who tame the subject rather than taunt the student. They approach the subject leisurely, allowing the reader to enjoy the view, and explain interesting sights along the way. There are others who use pictures in preference to words. And so on. Here, I have selected several books to act as guideposts on your journey to understand differential equations.

One book stands out as more of a popularization than a textbook. It is the book by the applied mathematician Steven Strogatz entitled Infinite Powers: How Calculus Reveals the Secrets of the Universe [1], exploring its history, techniques, significance, and applications. Make the effort to read this book—even if it is a little at a time—and complete it. You will become a better scholar as a result, and calculus will become your friend.

The other books are traditional textbooks on differential equations with different strengths like breadth, depth, clarity, historical context, etc. They include the texts by:

• Boyce and DiPrima [4];
• Blanchard, Devaney, and Hall [5];
• Edwards, Penney, and Calvis [6];
• Simmons [7];
• Krantz [8];
• Chasnov [9];
• Strang [10];
• Dettman [11]; and
• Howell [12].

These tomes embody a wealth of knowledge and approaches. My suggestion is that you should consult one or more of them to clear up specific doubts. If a concept or technique is not clear in any one book, look at the others. Perhaps the most reader-friendly book is the one by Howell [12]. If you are in a hurry to get an overview, the shortest presentation is Chasnov’s [9]. Some of these texts have web-based resources like videos that might better resonate with you. This way, you can gradually become self-reliant in any chosen field of study.

I particularly enjoyed an online article entitled Differential Equation Systems and Nature [13]. For those seeking a deeper understanding of the physics of skydving, I recommend this online presentation [14].

Acknowledgements

The resources provided gratis by Wolfram Alpha, ChatGPT, Copilot, Gemini, and Grok have facilitated this blog incalculably. They have also contributed to my learning experience on how (much) to depend on AI. My deepest gratitude to them.

I have freely adapted the figure TikZ Free Body Diagram of Skydiver with Parachute by Mohammed Benmiloud from his website to depict Figure 4. I express my gratitude for his work and recommend his site for those who wish to learn more about TikZ.

I also express my thanks to all the creators of the TikZlings packages for their delightful depictions of animals and beings.

Feedback

Since I work independently and alone, there is every chance that unintentional mistakes have crept into this blog, due to ignorance or carelessness. Therefore, I especially appreciate your corrective and constructive feedback.

Please email me your comments and corrections.

A PDF version of this article is available for download here:

References

•  AI (1)
•  analytical solution (1)
•  asymptote (1)
•  calculus (2)
•  differential equations (1)
•  Euler method (1)
•  numerical methods (1)
•  Simkin Tales (9)
•  slope fields (1)
•  terminal velocity (1)