Preface

The set of autonomous equations $ y' = f(y) $ is a subclass of separable differential equations. We discuss some specific properties of these equations.

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

Glossary

Autonomous Equations

Autonomous (meaning independent of time variable) equations are equations of the form

\begin{equation} \label{EqAuto.1} y' = f(y) \qquad\mbox{or}\qquad \frac{{\text d}y}{{\text d}t} = f(y) , \end{equation}

where slope function f(y) is a function of only dependent variable and does not involve independent variable explicitly. This equation encompasses, for example, all one-dimensional motions subject to conservative forces, but also many one-dimensional physical systems. This equation \eqref{EqAuto.1} occurs in cosmology, fluid mechanics, glaciology, hydrology, oceanography, and seismology.

We can make some observations regarding solutions to autonomous equations because their solutions have very limited types of behavior. Their certain mathematical properties are translated into physical properties relevant for natural phenomena. The integral curves (solutions) are invariant under translation along independent variable line, i.e., a horizontal shift of a solution is another solution. For example, consider an autonomous equation

\[ y' = 4 - y^2 . \]

We can plot a family of solutions using Mathematica:

sola = DSolve[{y'[x] == (4 - (y[x])^2)}, y[x], x];
g[x_, c_] = sola[[1, 1, 2]] /. C[1] -> c;
Plot[Table[g[x, c], {c, -5, 5}], {x, -2, 2}, PlotStyle -> Thick, PlotRange -> {-3, 3}, ImageSize -> 250]
StreamDensityPlot[{1, 4 - y^2}, {x, 0, 5}, {y, -6, 6}]
line = Graphics[Line[{{0, -2}, {0, 2}}], Epilog -> Inset[Text[Style["-2", Black, Medium]], {-0.2, -1}]];
p1 = Graphics[{PointSize[Large], Red, Point[{0, -1}]}];
p2 = Graphics[{PointSize[Large], Red, Point[{0, 1}]}];
a1 = Graphics[{Arrowheads[0.1], Arrow[{{0, -1.2}, {0, -2.2}}]}];
a2 = Graphics[{Arrowheads[0.1], Arrow[{{0, -1}, {0, 0.5}}]}];
a3 = Graphics[{Arrowheads[0.1], Arrow[{{0, 2.3}, {0, 1.3}}]}];
Show[line, a1, a2, a3, p1, p2]

All autonomous equations can be solved at least implicitly by separating variables:

\[ \frac{{\text d}y}{f(y)} = {\text d}x \qquad \Longrightarrow \qquad \int \frac{{\text d}y}{f(y)} = x+C, \quad f(y) \ne 0, \]

where C is an arbitrary constant. Since the above formula involves the reciprocal of the slope function, we have to exclude from our consideration all points where f vanishes---we have analyses this situation separately.

A point y = y^* where $ f\left( y^{\ast} \right) = 0 $ is called a critical point or equilibrium solution for the given differential equation $ y' = f(y) . $

Equilibrium solutions can be classified as follows.

Stable: The equilibrium solution y(t) = y^* is stable if all solutions with initial condition y(t₀) = y₀ "near" the critical point y^* approach it as t → ∞. This point is also called a sink.
Unstable: The equilibrium solution y(t) = y^* is unstable if all solutions with initial conditions "near" the critical point y^* do not approach it as t → ∞. This point is also called a source.
Semi-stable: The equilibrium solution y(t) = y^* is semistable if initial conditions y(t₀) = y₀ "near" the critical point on one side of it lead to solutions y(t) that approach y^* as t → ∞, while initial conditions on the other side of the critical point do not approach it. This point is also called a node.

A critical point is called globally stable if every solution approaches the equilibrium solution independently of the initial conditions.

If the slope function f(x) is continuous, the behavior of solutions of the autonomous equation can be determined from the slope lines along the vertical axis. This leads to construction of what is called a phase line for the differential equation.

Every critical point is a solution of the autonomous differential equation $ y' = f(y) , $ but solutions cannot cross it. If a solution touches the equilibrium solution, then this stationary solution is a singular solution. According to the existence theorem, if f(y) and ∂f(y)/ ∂y are continuous in some domain, then in this domain we are guaranteed that no other solution can intersect an equilibrium solution.

The nature of autonomous equations makes spotting constant solutions and interpreting the general behavior of solutions fairly straightforward by analyzing the slope function.

If the slope function is increasing below and decreasing above the critical point ⇒ Stable
If the slope function is decreasing below and increasing above the critical point ⇒ Unstable
If the slope function is decreasing below and above OR increasing below and above ⇒ Semi-stable

Theorem 1: Suppose that the slope function f(y) of the autonomous differential equation $ {\text d}y/{\text d}x = f(y) $ has a continuous derivative in some interval containing an equilibrium point y^*, that is, f(y^*) = 0.

If the derivative of f(y) of the slope function evaluated at the critical point is negative, i.e., $ f'\left( y^{\ast} \right) < 0, $ then the equilibrium solution y = y^* is stable and the equilibrium solution y^* is an attractor. This is a situation when neighbouring solutions that start above with y > y^* decrease toward it, but those neighbouring solutions that start below it with y < y^* increase towards the constant solution y^*.
If this derivative is positive, that is, $ f'\left( y^{\ast} \right) > 0, $ then the equilibrium solution y = y^* is unstable and the critical point is a repellor. This is a situation when neighbouring solutions that start above with y > y^* depart it, but those neighbouring solutions that start below it with y < y^* decrease away from the constant solution y^*.

When the derivative of the slope function evaluated at the critical point is zero, $ f'\left( y^{\ast} \right) = 0, $ then the above test is inconclusive. For example, the autonomous differential equation

\[ \dot{y} = f(y) = (y-1)\left( y-2 \right)^2 \]

has two critical points y = 1 (unstable and y = 2 semistable. The derivatives at these points are $ f' (1) = 1 > 0 $ and $ f' (2) = 0 . $

p1 = Graphics[{PointSize[Large], Red, Point[{0, 1}]}];
p2 = Graphics[{PointSize[Large], Red, Point[{0, 2}]}];
txt = Graphics[{Text[Style["1", Black, Medium], {-0.2, 1.01}],
Text[Style["unstable", Black, Medium], {2, 1.01}],
Text[Style["2", Black, Medium], {-0.2, 3.1}],
Text[Style["semi-stable", Black, Medium], {2.1, 3.1}]}];
line = Graphics[Line[{{0, 0}, {0, 3}}], Epilog -> Inset[txt]];
a1 = Graphics[{Arrowheads[0.1], Arrow[{{0, 0.8}, {0, -0.2}}]}];
a2 = Graphics[{Arrowheads[0.1], Arrow[{{0, 1.2}, {0, 1.9}}]}];
a3 = Graphics[{Arrowheads[0.1], Arrow[{{0, 2.2}, {0, 3.2}}]}];
Show[line, p1, p2, a1, a2, a3]
StreamPlot[{1, (y - 1)*(y - 2)^2}, {x, -3, 3}, {y, -1, 4}, StreamStyle -> Yellow, FrameStyle -> LightGray, Background -> Black]

Phase line for slope function (y-1)(y-2)².

On the other hand, the differential equation $ \dot{y} = f(y) = (y-1)^3 $ has one unstable critical point y = 1.

Theorem 2: Suppose that the slope function f(y) of the autonomous differential equation $ {\text d}y/{\text d}x = f(y) $ has two continuous derivatives in some interval containing an equilibrium point y^*, at which the slope function has a doublke root: $ f\left( y^{\ast} \right) = f'\left( y^{\ast} \right) = 0 , \quad\mbox{but} \quad f''\left( y^{\ast} \right) \ne 0. $

If $ f''\left( y^{\ast} \right) < 0, $ then the function f(y) has a maximum at y^* and $ y' = f(y) < 0 $ in neighbourhood of y^*. If neighbouring solutions that start below it y < y^*, they are running away from it, but neighbouring solutions that start above it, approach it in phase space. The equilibrium solution is semi-stable.
If $ f''\left( y^{\ast} \right) > 0, $ then the function f(y) has a minimum at y^* and $ y' = f(y) > 0 $ in neighbourhood of y^*. The reagion of stability and instability are reversed with respect to the previous case. The equilibrium solution is still semi-stable.

A critical point y = y^* of the autonomous differential equation $ y' = f(y) $ is referred to as being hyperbolic if the derivative of the slope function is not zero at the stationary point: $ f' \left( y^{\ast} \right) \ne 0 . $ Correspondingly, when $ f' \left( y^{\ast} \right) = 0 $ the case is called nonhyperbolic.

When you observe a higher order root of the slope function f(y), you must look at the lowest order derivative $ f^{(n)} (y) $ that does not vanish at the critical point y^*.

If n is odd, the sign of $ f^{(n)} \left( y^{\ast} \right) $ determines stability. So if $ f^{(n)} \left( y^{\ast} \right) < 0,$ it is stable, and if $ f^{(n)} \left( y^{\ast} \right) > 0,$ it is unstable.
If n is even, there is stability on one side and instability on the other side of the constant solution y = y^*.

Theorem 3: Let y(x) be a solution of the autonomous differential equation $ {\text d}y/{\text d}x = f(y) . $ Suppose that the slope function has a continuous derivative in some interval containing an equilibrium point y^* and that at some point x₀ the solution y(x) enters this interval.

If f(y) > 0 for y(x₀) ≤ y < y^* or f(y) < 0 for y^* < y ≤ y(x₀), then $ \lim_{x\to \infty} y(x) = y^{\ast} . $
If f(y) < 0 for y(x₀) ≤ y < y^* or f(y) > 0 for y^* < y ≤ y(x₀), then $ \lim_{x\to -\infty} y(x) = y^{\ast} . $

Here are three basic examples for illustration of each type of stability:

\[ \frac{{\text d}y}{{\text d}t} = \dot{y} = 1- y , \qquad \mbox{critical point $y=1$ is stable}. \]

\[ \dot{y} = y-1 , \qquad \mbox{critical point $y=1$ is unstable}. \]

\[ \dot{y} = \left( 1- y \right)^2 , \qquad \mbox{critical point $y=1$ is semi-stable}. \]

Theorem 4: Let f( y) be a continuous function on the closed interval [𝑎,b] that has only one null y^* ∈ (𝑎,b) and f( y) ≠ 0 for all other points y ∈ ([𝑎,b). If the integral

\[ \int_{y}^{y^\ast} \frac{{\text d}y}{f(y)} \]

diverges, then the initial value problem for the autonomous differential equation $ y' = f(y), \quad y(t_0 ) = y^{\ast} $ has the unique solution y( t) = y^*. If the above integral converges, then the initial value problem has multiple solutions.

Example: We want to find the velocity of the falling parachutist as a function of time t and are particularly interested in the constant limiting velocity, v_∞, due to air resistance. We assume that air drag is proportional to the square of the velocity, -k v², and opposing the force of the gravitational attraction, mg, of the Earth. We choose a coordinate system in which the positive direction is downward so that the gravitational force is positive. For simplicity we assume that the parachute opens immediately, that is, at time t = 0, where v = 0. Newton’s law applied to the falling parachutist of total mass m (including chute) gives

\[ m\,\dot{v} = mg - k\,v^2 , \qquad v(0 ) = 0 , \]

where g is the acceleration due to gravity. As usual, the dot indicates the derivative with respect to time variable: $ \dot{v} = {\text d}v/{\text d}t . $

The terminal velocity, v_∞, can be found from the equation of motion as t → ∞. When there is no acceleration, its velocity is zero, so

\[ m\,v_{\infty}^2 = mg , \qquad \mbox{or} \qquad v_{\infty} = \sqrt{\frac{m\,g}{k}} . \]

Separating variables t and v, we get

\[ \frac{{\text d}v}{g - k\,v^2 /m} = {\text d} t \qquad \Longrightarrow \qquad \int \frac{{\text d}v}{g - k\,v^2 /m} = t + c . \]

Using partial fraction decomposition, we get

\[ \frac{1}{g - k\,v^2 /m} = \frac{m}{2\,v_{\infty} k} \left( \frac{1}{v + v_{\infty}} - \frac{1}{v - v_{\infty}} \right) . \]

Integrating both terms yields

\[ \frac{{\text d}v}{g - k\,v^2 /m} = \frac{1}{2}\,\sqrt{\frac{m}{k\,g}} \, \ln \frac{v_{\infty} +v}{v_{\infty} -v} = t , \]

where the constant of integration was set to zero in order to satisfy the initial condition. Solving for the velocity gives

\[ v = v_{\infty} \frac{e^{2t/s} -1}{e^{2t/s} +1} = v_{\infty} \frac{\sinh \frac{t}{s}}{\cosh \frac{t}{s}} = v_{\infty} \tanh \frac{t}{s} , \]

where $ \displaystyle s = \sqrt{\frac{m}{k\,g}} $ is the time constant governing the asymptotic approach of the velocity to the limiting velocity v_∞. ■

Example: Solve the initial value problem

\[ \dot{y} = \left( y - 1 \right)^2 , \qquad y(t_0 ) = y_0 . \]

We find its solution with Mathematica:

DSolve[{y'[t] == (y[t] - 1)^2 , y[a] == A}, y[t], t]

\[ y(t) = \frac{y(t_0 ) - t_0 + t_0 y(t_0 ) + t - t\,y(t_0 )}{1- t_0 + t_0 y(t_0 ) + t - t\,y(t_0 )} . \]

When y(t₀) = y₀ ≠ 0, then the integral curve is a hyperbola. This solution blows up at

\[ t = \frac{t_0 -1 -t_0 y(t_0 ) }{1- y(t_0 )} \qquad \mbox{if}\quad y(t_0 ) \ne 1. \]

When y₀ = 1, we get the integral curve to be nonsingular equilibrium solution. ■

Example: Consider the differential equation

\[ \dot{y} ≝ {\text d}y/{\text d}t = \left( y-1 \right) \left( y-2 \right)^3 . \]

The slope function f(y) = (y-1)(y-2)³ has two critical points y = 1 (stable) and y = 2 (unstable). Evaluating its derivatives at these points, we get

\[ f'(1) = \]

f[y_] = (y - 1)*(y - 2)^3;
D[f[y], y] /. y -> 1

-1

D[f[y], y] /. y -> 2

Since the derivative test at y = 2 is inconclusive, we draw a phase line:

StreamPlot[{1, (y - 1)*(y - 2)^3}, {x, -2, 2}, {y, -1, 4}, StreamPoints -> 24];
p1 = Graphics[{PointSize[Large], Red, Point[{0, 1}]}];
p2 = Graphics[{PointSize[Large], Red, Point[{0, 2}]}];
txt = Graphics[{Text[Style["1", Black, Medium], {-0.2, 1.01}], Text[Style["stable", Black, Medium], {2, 1.01}],
Text[Style["2", Black, Medium], {-0.2, 3.1}], Text[Style["unstable", Black, Medium], {2.1, 3.1}]}];
line = Graphics[Line[{{0, 0}, {0, 3}}], Epilog -> Inset[txt]];
a1 = Graphics[{Arrowheads[0.1], Arrow[{{0, -0.2}, {0, 0.8}}]}];
a2 = Graphics[{Arrowheads[0.1], Arrow[{{0, 1.9}, {0, 1.1}}]}];
a3 = Graphics[{Arrowheads[0.1], Arrow[{{0, 2.2}, {0, 3.2}}]}];
Show[line, p1, p2, a1, a2, a3]

Phase line for slope function (y-1)(y-2)³.

■

Example:

Karl Ludwig von Bertalanffy (1901--1972) was an Austrian biologist known as one of the founders of general systems theory. This is an interdisciplinary practice that describes systems with interacting components, applicable to biology, cybernetics and other fields. He made valuable contributions to the study of organism growth, including the growth of tumors, based on the relationship between body size of the organism and the metabolic rate. Ludwig theorized that weight is directly proportional to the volume and that the metabolic rate is proportional to surface area. This indicates that surface area is proportional to V^2/3 because V is measured in cubic units and its area in square units. Hence, Bertalanffy studied the initial value problem

\[ \frac{{\text d}V}{{\text d}t} = a\,V^{2/3} - b\,V , \qquad V(0) = V_0 , \]

where 𝑎 and b are positive constants. Separation variables yields

\[ \frac{{\text d}V}{a\,V^{2/3} - b\,V} = {\text d}t \qquad \Longrightarrow \int \frac{{\text d}V}{a\,V^{2/3} - b\,V} = t+C . \]

We ask Mathematica for help with integration:

Integrate[1/(a*v^(2/3) - b*v), v]

-((3 Log[a - b v^(1/3)])/b)

This gives us the general solution in implicit form:

\[ \ln \left\vert a - b\,V^{1/3} \right\vert = C - \frac{b}{3}\,t . \]

Raising to exponent and using the initial condition, we obtain the explicit solution

\[ V(t) = \left[ \frac{a}{b} - \left( \frac{a}{b} - V_0^{1/3} \right) e^{-bt/3} \right]^3 . \]

Since the solution approaches (𝑎/b)³ as t → ∞, the horizontal line $ V(t) = (a/b)^3 $ is the upper fence for solutions of the Bertalanffy model. We plot the corresponding phase portrait:

ludwig = StreamPlot[{1, 2*y^(2/3) - y}, {x, 0, 10}, {y, -0.3, 8.1}, Frame -> False, Axes -> True, AspectRatio -> 1];
line = Plot[8, {x, 0, 10}, PlotStyle -> {Thick, Green}, PlotLabels -> {"upper fence"}]
Show[line, ludwig, PlotRange -> {0, 8.2}]

■

Example: The bone marrow of a human constantly produces red blood cells (erythrocytes) because they have a finite life span (approximately 120 days). An ordinary differential equations model of bone marrow erythrocytes was developed in the literature presented below. The model consists of concatenated cell compartments describing the dynamics of the bone marrow cell stages of erythrocytes. Differentiation of cells is modeled by fluxes from one compartment to the next stage. Proliferation is modeled by amplification in the corresponding compartment. Compartment sizes are modeled by balance equations of the form

\[ \frac{\text d}{{\text d}t}\,C = C^{in} - \frac{A}{T}\, C , \]

where A is the amplification, T is the transit time, C,sup>in is the influx from the previous compartment, and C is the content of the compartment. The oxygen level in kidney tissue determines the production of erythropoietin (Epo), and the level of Epo is proportional to the rate of red cell production. This means that the oxygen level is inverse proportional to the rate of erythrocyte production. This leads to the following differential equation model:

\[ \frac{\text d}{{\text d}t}\,u = -a\, u(t) + \frac{b}{u(t)} = \frac{b-a\,u^2 }{u} \qquad (u \ne 0), \]

with positive constants 𝑎 and b. This autonomous differential equation has two stable equilibrium solutions (of which only positive one has a physical meaning) y = ±b/𝑎 because the derivative of the slope function evaluated at these points is negative.

StreamPlot[{1, -(y^2 - 1)/y}, {x, -3, 3}, {y, -2, 2}, StreamStyle -> Yellow, FrameStyle -> LightGray, Background -> Black]
Simplify[D[(b - a*u^2)/u, u] /. u -> (b/a)]

-((a (a + b))/b)

Note that u = 0 corresponds to the vertical line and this point is not a singular point for the red cell model because at this point the differential equation has just horizontal slope:

\[ \frac{{\text d}t}{{\text d}u} = \frac{u}{b - a\, u^2} . \]

We separate variables to obtain

\[ \frac{{\text d}u}{-a\,u + b/u} = {\text d}t \qquad (-a\,u + b/u \ne 0). \]

Integration yields

\[ t + C = \int \frac{u\,{\text d}u}{b - a\,u^2} = -\frac{1}{2a} \int \frac{{\text d} \left( b - a\, u^2 \right)}{b - a\,u^2} = -\frac{1}{2a} \, \ln \left\vert b - a\, u^2 \right\vert , \]

with a constant of integration C. Raising both sides of the above equation to an exponent, we get

\[ b - a\, u^2 = a\,c\, e^{-2a\,t} , \]

with another arbitrary constant c. Hence,

\[ u (t) = \left[ \frac{b}{a} - c\,e^{-2a\,t} \right]^{1/2} \]

is the general solution of the differential equation that serves as a mathematical model of the red cell system.

Loeffler M, Pantel K, Wulff H, Wichmann H (1989) A mathematical model of erythropoiesis in mice and rats. part 1: Structure of the model. Cell Tissue Kinet 22: 13–30.
Pantel K, Loeffler M, Bungart B, Wichmann H (1990) A mathematical model of erythropoiesis in mice and rats. part 4: Differences between bone marrow and spleen. Cell Tissue Kinet 23: 283–297.
Schirm, S., Engel, C., Loeffler, M., Scholz, M., A Biomathematical Model of Human Erythropoiesis under Erythropoietin and Chemotherapy Administration, 2013, PLoS One, Vol. 8, No. 6, e65630. doi: 10.1371/journal
Wichmann H (1983) Computer modeling of erythropoiesis. In: Current Concepts in Erythropoiesis, Chapter V., John Wiley and Sons.
Wichmann H, Loeffler M (1985) Mathematical modeling of cell proliferation: Stem cell regulation in hemopoiesis, Vol. 1, 2. CRC Press.

■

Example: A catenary is the curve that an idealized hanging chain or cable assumes under its own weight when supported only at its ends. The word "catenary" is derived from the Latin word catēna, which means "chain". The English word "catenary" is usually attributed to Thomas Jefferson. The catenary is also the surface of revolution, the catenoid, having a minimal surface of revolution. The mathematical properties of the catenary curve were first studied by Robert Hooke in the 1670s, and its equation was derived by Leibniz, Huygens, and Johann Bernoulli in 1691.

A chain hanging from points forms a catenary.

Let y(x) be a positive function defined on the interval (x₁, x₂). Let A(y) be the area of the surface of revolution that we get by rotating the curve y(x) around the x-axis. From calculus, we know that

\[ A(y) = 2\pi \int_{x_1}^{x_2} {\text d} x \,y(x)\,\sqrt{1 + y'(x)^2} . \]

Therefore, we need to find a curve y(x) that minimizes the area A(y and satisfies the boundary conditions $ y(x_1 ) = y_1 , \quad y(x_2 ) = y_2 . $

Let us introduce the Lagrangian:

\[ {\cal L} (y, y' ) = 2\pi\, y\,\sqrt{1 + y'^2} . \]

Then

\[ \frac{\partial {\cal L}}{\partial y} = 2\pi\, \sqrt{1 + y'^2} , \qquad \frac{\partial {\cal L}}{\partial y'} = 2\pi\,y\, \frac{y'}{\sqrt{1 + y'^2}}, \]

and from the Euler--Lagrange equation

\[ \frac{\partial {\cal L}}{\partial y} = \frac{\text d}{{\text d}t} \,\frac{\partial {\cal L}}{\partial y'} , \]

it follows

\[ 2\pi\,\sqrt{1 + y'^2} = \frac{\text d}{{\text d}t}\left( 2\pi\,y\, \frac{y'}{\sqrt{1 + y'^2}} \right) . \]

Any stationary curve y(x) is a solution of the above nonlinear second order differential equation. Since the above Lagrangian does not explicitly depend on x, we have

\[ y'\, \frac{\partial {\cal L}}{\partial y'} - {\cal L} = c , \]

for some constant c. Now we immediately obtain the catenary equation

\[ - \frac{y}{\sqrt{1 + y'^2}} = c_1 ≝ \frac{c}{2\pi} . \]

Since y(x) is positive, the latter can have a solution only when c₁ is negative. Without any loss of generality, we can set α = -c₁, where α > 0. Using this, we have

\[ \frac{y}{\sqrt{1 + y'^2}} = \alpha \qquad \Longrightarrow \qquad y' = \pm \sqrt{\left( \frac{y}{\alpha} \right)^2 -1} . \]

This equation is separable and we can find its its solution to be

\[ y(x) = \alpha\, \cosh \left( \frac{x}{\alpha} + \beta \right) , \qquad \alpha > 0, \quad \beta \in \mathbb{R} , \]

where the constants α and β are determined from the constraints:

\[ y \left( x_1 \right) = y_1 , \qquad y \left( x_2 \right) = y_2 . \]

The graph of the above cosine hyperbolic curve is called a catenary. Then the catenary curve that goes through the specified point (x₁, y₁) becomes

\[ y(x) = \alpha\, \cosh \left( \frac{x-x_1}{\alpha} + \mbox{arccosh} \frac{y_1}{\alpha} \right) . \]

In order to make sure that y(x) also pass through the point (x₂, y₂), we must determine an α such that

\[ y_2 = \alpha\, \cosh \left( \frac{x_2 -x_1}{\alpha} + \mbox{arccosh} \frac{y_1}{\alpha} \right) . \]

This is a transcendent equation for α and its solution depends on coordinates of given two points (x₁, y₁) and (x₂, y₂). The following cases are known to occur.

There exists no Catenary connecting the points (x₁,y₁) and (x₂,y₂).
There exists exactly one Catenary connecting the points (x₁,y₁) and (x₂,y₂).
Milk There exists exactly two Catenaries connecting the points (x₁,y₁) and (x₂,y₂).

In case ii) the unique catenary is neither a global nor local minimum. In case iii), one of the solutions is a local minimum, and this is also a global minimum, if y₁ is large enough in a certain sense. Therefore, a surface of revolution of minimum area exists only in case iii) and only if y₁ is large enough.

Chain's equation

Suppose that a heavy uniform chain is suspended at points Mi>A,B, which may be at different heights (see figure). Consider equilibrium of a small element of the chain of length Δs. The forces acting on the section of the chain are the distributed force of gravity

a = Plot[Cosh[x*2] + 1, {x, -1, 1.2}, Axes -> False, PlotRange -> {{-1.5, 2.5}, {-0.5, 7.5}}, PlotStyle -> {Thickness[0.015], Blue}];
b = Plot[Cosh[x*2] + 1, {x, 0.5, 0.75}, Axes -> False, PlotRange -> Full, PlotStyle -> {Thickness[0.015], Green}]
line1 = Graphics[Line[{{1.1, 2.5430806}, {-0.1, 2.5430806}}]]
line2 = Graphics[Line[{{0.75, 3.3524}, {1.5, 3.3524}}]];
arrow1 = Graphics[Arrow[{{0.5, 2.5430806}, {-0.1, 0.99}}]];
arrow2 = Graphics[Arrow[{{0.75, 3.3524}, {1.5, 5.9}}]];
disk1 = Graphics[{Orange, Disk[{-1, 4.7633}, {0.05, 0.2}]}];
disk2 = Graphics[{Orange, Disk[{1.2, 6.55695}, {0.05, 0.2}]}];
textA = Graphics[Text[Style["A", FontSize -> 18, Red], {-1.0, 5.5}]];
textB = Graphics[Text[Style["B", FontSize -> 18, Red], {1.35, 6.95}]];
txt1 = Graphics[ Text[Style["\[Alpha](x)", FontSize -> 16, Black], {-0.35, 1.7}]];
txt2 = Graphics[ Text[Style["\[Alpha](x+\[CapitalDelta]x)", FontSize -> 16, Black], {1.35, 3.9}]];
textB = Graphics[Text[Style["B", FontSize -> 18, Red], {1.35, 6.75}]];
textB = Graphics[Text[Style["B", FontSize -> 18, Red], {1.35, 6.75}]];
txt3 = Graphics[ Text[Style["\[CapitalDelta]s", FontSize -> 16, Black], {0.45, 3.1}]];
txt4 = Graphics[ Text[Style["\[CapitalDelta]y", FontSize -> 16, Black], {1.0, 3.0}]];
arrow3 = Graphics[{Arrowheads[Large], Arrow[{{0.0, -0.2}, {0.0, 7.5}}]}];
arrow4 = Graphics[{Arrowheads[Large], Arrow[{{-0.5, 0.0}, {1.8, 0.0}}]}];
textX = Graphics[ Text[Style["x", FontSize -> 18, Black], {1.75, 0.5}]];
textY = Graphics[Text[Style["y", FontSize -> 18, Black], {0.2, 7.0}]];
line3 = Graphics[{Dashed, Line[{{0.5, 2.5430806}, {0.5, 0}}]}];
line4 = Graphics[{Dashed, Line[{{0.75, 3.5}, {0.75, 0}}]}];
txt5 = Graphics[Text[Style["x", FontSize -> 14, Black], {0.4, 0.3}]];
txt6 = Graphics[ Text[Style["x+\[CapitalDelta]x", FontSize -> 14, Black], {1.0, 0.3}]];
arrow5 = Graphics[{, Green, Arrowheads[Large], Arrow[{{0.65, 3.0 \[CapitalDelta]x}, {0.65, 1.4}}]}];
txt7 = Graphics[ Text[Style["\[CapitalDelta]p", FontSize -> 16, Black], {0.65, 1.0}]];
txt8 = Graphics[ Text[Style["T(x)", FontSize -> 16, Black], {-0.3, 0.6}]];
txt8 = Graphics[ Text[Style["T(x)", FontSize -> 16, Black], {-0.3, 0.6}]];
Show[a, b, line1, line2, arrow1, arrow2, disk1, disk2, textA, textB, \ txt1, txt2, txt3, txt4, txt5,
txt6, txt7, txt8 , txt9, arrow3, \ arrow4, arrow5, textX, textY, line3, line4]

\[ \Delta P = \rho\,g\,\Delta s , \]

where ρ is the density of the chain material, g is the acceleration of gravity, A is the cross sectional area of the thread, and the tension forces T(x) and T(x+Δx), respectively, at points x and x+Δx.

The equilibrium conditions of the element of the length &Deltas for projections on the axes Ox and Oy are turned out to be

\[ \begin{split} -T\left( x\right) \cos\alpha (x) + T\left( x + \Delta x \right) \cos \alpha (x+\Delta x) &= 0, \\ -T\left( x\right) \sin\alpha (x) + T\left( x + \Delta x \right) \sin \alpha (x+\Delta x) - \Delta P &= 0 \end{split} \]

It follows from the first equation that the horizontal component of the tension force T(x) is always a constant:

\[ T\left( x\right) \cos\alpha (x) = T_0 = \mbox{constant} . \]

Using differentials in the second equation we can rewrite it in differentials:

\[ {\text d}\left( T(x)\,\sin \alpha(x)\right) = {\text d} P(x) . \]

Since $ \displaystyle T(x) = \frac{T_0}{\cos \alpha (x)} , $ we have

\[ {\text d}\left( T_0 \tan \alpha (x) \right) = {\text d} P(x) \qquad \Longrightarrow \qquad T_0 {\text d}\left( \tan \alpha (x) \right) = {\text d} P(x) . \]

Since the tangent is just the derivative of its curve function, we get

\[ T_0 {\text d}\left( y' \right) = {\text d} P(x) \qquad \Longrightarrow \qquad T_0 {\text d}\left( y' \right) = \rho g\,A\,{\text d}s . \]

The element of of curve length can be expressed as

\[ {\text d} s = \sqrt{1 + \left( y' \right)^2}\,{\text d}x . \]

As a result we obtain the differential equation of the catenary:

\[ T_0 \frac{{\text d} y'}{{\text d} x} = \rho g A \sqrt{1 + \left( y' \right)^2} \qquad \Longrightarrow \qquad T_0 y'' = \rho g A \sqrt{1 + \left( y' \right)^2} . \]

The order of this equation can be reduced by introducing anew variable $ \displaystyle u= y' . $ Then we have

\[ T_0 u' = \rho g A \sqrt{1 + u^2} . \]

Separation of variables yields

\[ T_0 {\text d}u = \rho g A \sqrt{1 + u^2} \,{\text d} x \qquad \Longrightarrow \qquad \frac{{\text d} u}{\sqrt{1 + u^2}} = \frac{\rho g A}{T_0}\,{\text d}x . \]

Then we integrate and obtain

\[ \ln \left( u + \sqrt{1+u^2} \right) = \frac{x}{a} + C_1 , \qquad \mbox{where } a= \frac{T_0}{\rho g A} . \]

The arbitrary constant C₁ can be determined from the condition that the catenary at the lowest point is parallel to the x-axis. Hence,

\[ u(x=0) = y' (x=0) = 0 \qquad \Longrightarrow \qquad C_1 = 0 . \]

Thus, we get the following equation:

\[ u + \sqrt{1+ u^2} = e^{x/a} . \]

Multiplying both sides of the equation by the conjugate expression $ \displaystyle u - \sqrt{1+u^2} $ gives

\[ \left( u + \sqrt{1+u^2} \right) \left( u - \sqrt{1+u^2} \right) = \left( u - \sqrt{1+u^2} \right)e^{/a} . \]

Upon simplification, we have

\[ u - \sqrt{1+u^2} = e e^{-x/a} . \]

Adding to the previous equation, we find the expression for $ \displaystyle u = y' : $

\[ u + \sqrt{1+u^2} + u - \sqrt{1+u^2} = e^{x/a} - e^{-x/a} \qquad \Longrightarrow \qquad u = y' = \sinh \frac{x}{a} . \]

Integrating once more gives the final nice expression for the shape of the catenary:

\[ y = a\,\cosh \frac{x}{a} . \]

The catenary has another interesting feature. When revolved about the x-axis, the catenary gives the surface called catenoid. This surface has minimum surface area, i.e. any part of the catenoid will be less than any other surface bounded by the same contour. Catenoid was first described in 1744 by the mathematician Leonhard Euler. In particular, the soap film between two circles trying to minimize the free energy takes the form of a catenoid. ■

Example: Consider the initial value problem

\[ \frac{{\text d}y}{{\text d}x} = e^y , \qquad y(0) = c \quad (\mbox{a constant}). \]

Separating variales and integrating, we obtain its explicit solution

\[ y(x) = - \ln \left( e^{-c} - x \right) , \qquad -\infty < x < e^{-c} . \]

So choosing a different value c in the initial condition, we obtain a different domain of existence. ■

Example: Consider a simple initial value problem for an autonomous equation

\[ \frac{{\text d}y}{{\text d}x} = (y-1)^2 , \qquad y(0) =1. \]

Using separation of variables

\[ \frac{{\text d}y}{(y-1)^2} = {\text d}x \qquad \Longrightarrow \qquad -\frac{1}{y-1} = x + C , \]

we find the general solution in explicit form (which is a very rare case)

\[ y = 1 + \frac{1}{x+C} . \]

To satisfy the initial condition, we are forced to set the arbitrary constant to be infinity and the required solution y = 1. If we make a 1% change in the initial condition to y(0) = 1.01, we obtain a hyperbola

\[ y = 1 + \frac{1}{x-100} , \]

with one vertical asymptote. If we perturb the driving input again by 1%,

\[ y' = (y-1)^2 + 0.01, \]

we get a periodic solution

\[ y = 1 + 0.1\,\tan \left( \frac{x}{10} \right) , \]

with vertical asymptotes at (2n+1)5π, n = 0, ±1, ±2, .... Hence, we see how sensitive a nonlnear differential equation can be to small perturbations.

■

Arino, O. and Kimmel, M., Stability analysis of models of cell production systems, Mathematical Modeling, 1986, Vol. 7, pp. 1269--1300.

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

MATHEMATICA
TUTORIAL

under the terms of the GNU General Public License ( GPL)

for the First Course. Part II: Autonomous Equations

Email: Prof. Vladimir Dobrushkin ()

Preface

Contents [hide]

Glossary

Autonomous Equations

MATHEMATICA TUTORIAL under the terms of the GNU General Public License ( GPL) for the First Course. Part II: Autonomous Equations

Email: Prof. Vladimir Dobrushkin ()

Preface

Contents [hide]

Glossary

Autonomous Equations

MATHEMATICA
TUTORIAL

under the terms of the GNU General Public License ( GPL)

for the First Course. Part II: Autonomous Equations