text:
# Preface

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

# Autonomous Equations

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

$$\label{EqAuto.1} y' = f(y) \qquad\mbox{or}\qquad \frac{{\text d}y}{{\text d}t} = f(y) ,$$
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]
 Family of general solutions. Phase portrait. Phase line for slope function 1-y².

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(t0) = y0 "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(t0) = y0 "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 portrait for (y-1)(y-2)². 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 x0 the solution y(x) enters this interval.
• If f(y) > 0 for y(x0) ≤ y < y* or f(y) < 0 for y* < yy(x0), then $$\lim_{x\to \infty} y(x) = y^{\ast} .$$
• If f(y) < 0 for y(x0) ≤ y < y* or f(y) > 0 for y* < yy(x0), 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] -> (-a + A + a A + t - A t)/(1 - a + a A + t - A 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(t0) = y0 ≠ 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 y0 = 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
0
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 portrait for (y-1)(y-2)³. 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 V2/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.

Let y(x) be a positive function defined on the interval (x1, x2). 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 c1 is negative. Without any loss of generality, we can set α = -c1, 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 (x1, y1) 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 (x2, y2), 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 (x1, y1) and (x2, y2). The following cases are known to occur.
1. There exists no Catenary connecting the points (x1,y1) and (x2,y2).
2. There exists exactly one Catenary connecting the points (x1,y1) and (x2,y2).
3. Milk There exists exactly two Catenaries connecting the points (x1,y1) and (x2,y2).
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 y1 is large enough in a certain sense. Therefore, a surface of revolution of minimum area exists only in case iii) and only if y1 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] A chain hanging from points A and B. \[ \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(xx), respectively, at points x and xx.

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 C1 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.

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


text:
## Math 4 Unit 4: Rational Functions

#### Graphing Rational Functions

Identify the holes, intercepts, HA, VA, and SA, then graph.

#### Multiplying Expressions

Multiply and simplify.

#### Dividing Expressions

Divide and simplify.

#### Subtracting Expressions

Subtract and simplify.

#### Solving Inequalities

Solve the inequality.

## Math 4 Unit 3: Polynomials

Completely factor x^2+6x+8

= (x+4)(x+2)

Completely factor 3x^2+4x-7

= (3x+7)(x-1)

Completely factor

= x^2(x-7) – 3(x-7)

= (x^2-3)(x-7)

#### Factoring Quartics

Completely factor 2x^4-28x^2+90

= 2(x^4-14x^2+45)

= 2(x^2-9)(x^2-5)

= 2(x+3)(x-3)(x^2-5)

#### Solving Quartics

Completely factor.

= 2x^4-28x^2+90

= 2(x^4-14x^2+45)

= 2(x^2-9)(x^2-5)

= 2(x+3)(x-3)(x^2-5)

x=-3, 3, +-sqrt 5

#### Vertex Form

Rewrite in vertex form.

#### Graphing Polynomials

• Determine the End Behavior based on the Leading Coefficient Test and write the end behavior of this polynomial in Limit Notation.
• Find the zeros of the polynomial function
• Find the y-intercept of the polynomial function
• Graph the polynomial function.

#### Write a Polynomial with Roots

Find a polynomial function with the given zeros: 0, 2 mult 2

#### Long Division

Divide with Long Division

#### Synthetic Division

Divide with Synthetic Division

#### Rational Roots Theorem

Use the Rational Roots Theorem to find all possible roots and then solve.

## Math 4 Unit 2: Functions and Their Graphs

#### Vocab

Domain – the span of all the x-values of a function

Range – the span of all the y-values of a function

x-intercept(s) – the point(s) where a function touches the x-axis

y-intercept – the point where a function touches the y-axis

Horizontal asymptote(s) – the imaginary line(s) where the domain of a rational function approaches on both sides but never touches.

• There are three rules that define where the asymptote will be and how it will look like:
• If the degree of the numerator is less than the degree of the denominator, the asymptote will, by default, be y = 0.
• The asymptote will be higher or lower if there is a k value in the function.
• If the degree of the numerator is equal to the degree of the denominator, the asymptote will, by default, be at y = (degree of the numerator/degree of the denominator)
• The asymptote will be higher or lower if there is a k value in the function.
• If the degree of the numerator is greater than the degree of the denominator, there will not be a horizontal asymptote.

Vertical asymptote – the imaginary line(s) where the range of a rational function approaches above and below but never touches.

Local minimum – the point where a positive even polynomial’s vertex is, where the function does not go any lower.

Local maximum – the point where a negative even polynomial’s vertex is, where the function does not go any higher.

End behavior – the direction, attributes, and/or limit of a function as it goes left or right.

#### Application

The height of a passenger on an escalator is modeled as

where t is the time in seconds after the passenger has boarded the escalator on the second floor and h is the height in feet above the first floor.

This means it takes the passenger 16 seconds to reach the first floor (0 feet).

This means the passenger is 12 feet above the first floor when he/she first gets on the escalator. Likewise, this also means that the second floor is 12 feet above the first floor.

#### Inverse Functions

Find the inverse of the following function:

#### Composition of Functions

Find f(g(x)) and its domain if:

and:

#### Decomposition of Functions

Find what h(x) decomposes into if:

and:

## Math 4 Unit 1: Sequences and Series

#### Recursive Sequence

A sequence expressed in the form of an operation that is performed on the previous term of the sequence to find the next. To use a recursive sequence, the previous term must be given to calculate any following terms.

Example: Find the first 4 terms of the recursive sequence below.

#### Explicit Sequence

A sequence expressed in the form of an equation in terms of n, where n is the term’s position in the sequence. This method of expressing a sequence does not require any previous terms to find following terms; it only requires you have the position of the term (n) you are looking for and the first term of the sequence.

Example: Find the first 4 terms of the explicit sequence below.

#### Arithmetic Sequence

Example: Find the 4th term of an Arithmetic sequence with a first term of 7 and a common difference of -4.

#### Arithmetic Series

Example: Find the sum of the first 10 terms of an Arithmetic sequence with a first term of 7 and a common difference of -4.

#### Geometric Sequence

Example: Find the 6th term of a Geometric sequence with a first term of 2 (a1) and a common ratio of 3 (r).

#### Geometric Series

Example: Find the 6th term of a Geometric series with a first term of 2 (a1) and a common ratio of 3 (r).

#### Application Problem

A cybersecurity technician has a starting salary \$106,000. Per his contract with the company, he is given a 5% raise each year. If he works for 20 years, how much will he have earned total by the end?

a1 = 106,000 because the salary starts at \$106,000

r = 105% or 1.05 because the salary increases by 5% each year

n = 20 because the technician is working 20 years

#### Infinite Sequences

Find the sum of the following sequence, if there is one:

## 1/20, 1/40, 1/80…

#### Convergent & Divergent Series

A convergent series has an r value with an absolute value less than 1, meaning the sequence will get closer and closer to 0 as it progresses. This means you can calculate the sum of the series to the infinite term, as there will be a finite sum.

For example, the geometric series in the previous question has an r value of 1/2 (|1/2|<1), because its terms are getting closer and closer to 0 (in this case, smaller and smaller)

On the other hand, a divergent series has an r value with an absolute value greater than 1, meaning the sequence will get farther and farther away from 0 as it progresses. Because of this, you cannot find the sum to an infinite term as the sequence has no limit to how large its terms it can get.

For example, this geometric series has an r value of 2 (|2|>1), which means the terms will get larger and larger, farther and farther away from 0, meaning there is no real way to find an infinite total of the series since there is no limit to how high it will go.

#### Summation Notation

Find the sum using the correct sequence formula.

#### Binomial Theorem (all terms)

Expand completely.

#### Binomial Theorem (single term)

Find the 3rd term in expansion of

format: plain


text:
# Why are these two functions equal?

These functions are equal. But I don't understand why.

$$a \leftrightarrow f(x) =|\cos(2\pi x)|^2$$

$$b \leftrightarrow f(x) = \dfrac{\cos(4\pi x)}{2} + 0.5$$

Which results both in this plot:

-

They are the same because $$|\cos(x)|^2=\cos^2(x)=\frac{1+\cos(2x)}{2}=\frac{\cos(2x)}{2}+\frac{1}{2}.$$ (Note that the absolute value doesn't do anything, because of the square; i.e., $x^2=|x|^2$ for any $x$.) See here. This can be derived as follows:

For any $x$ and $y$, we have $$\cos(x+y)=\cos(x)\cos(y)-\sin(x)\sin(y).$$ Thus, $$\cos(2x)=\cos^2(x)-\sin^2(x)$$ $$\cos(2x)=\cos^2(x)-(1-\cos^2(x))$$ $$\cos(2x)=2\cos^2(x)-1$$ and therefore $$\cos^2(x)=\frac{1+\cos(2x)}{2}$$

Here is another, more direct derivation, using complex exponential definition of cosine: $$\cos^2(x)=\left(\frac{e^{ix}+e^{-ix}}{2}\right)^2=\frac{(e^{ix})^2+2(e^{ix})(e^{-ix})+(e^{-ix})^2}{4}=$$ $$\frac{e^{i(2x)}+e^{-i(2x)}+2}{4}=\frac{\cos(2x)+1}{2}$$

-
Thanks. I can accept, but not +1... "An error has occurred - please retry your request" I already tried 5 times. –  Martijn Courteaux Jun 27 '11 at 8:51
@Martijn: Hmm, that's a bit odd. I guess try refreshing the page, or restarting your browser? If it doesn't work, don't worry about it, one upvote is not much in the grand scheme of things :) –  Zev Chonoles Jun 27 '11 at 8:53
Really strange, I restarted the browser already. I can't down vote as well (just to try ;)). –  Martijn Courteaux Jun 27 '11 at 8:57
Aha, succeeded. I logged out from the site, and logged in again. Now, I can see my upvote was registered, but I was unable to see it. –  Martijn Courteaux Jun 27 '11 at 9:00
Another way to see they are equal is to consider the derivative of their difference. It's zero, so the difference is a constant, and the constant can be found by considering the difference at any convenient value such as x=0. –  tzs Jun 27 '11 at 11:40


text:
# A sodium chloride crystal in the shape of a cube is expanding at the rate of 60 cubic microns per second. How fast is the side of the cube growing when the volume is 1000 cubic microns?

A sodium chloride crystal in the shape of a cube is expanding at the rate of 60 cubic microns per second. How fast is the side of the cube growing when the volume is 1000 cubic microns?
You can still ask an expert for help

• Questions are typically answered in as fast as 30 minutes

Solve your problem for the price of one coffee

• Math expert for every subject
• Pay only if we can solve it

Dominique Ferrell
The sodium chloride is a cube.
$V={s}^{3}$
Take the derivative with respect to time, t.
$\frac{dV}{dt}=3{s}^{2}\frac{ds}{dt}$
We know that:
$\frac{dV}{dt}=60$
We want to find $\frac{ds}{dt}$ when V=1000, which means that s=10.
$60=3{\left(10\right)}^{2}\frac{ds}{dt}$
$60=300\frac{ds}{dt}$
$\frac{ds}{dt}=\frac{1}{5}$ microns/second
When the cube has a volume of 1000 ${\text{microns}}^{3}$, the side of the cube is growing at a rate of $\frac{1}{5}$ $\text{microns/second}$.


text:
# Horizontal alignment within itemize

My minimum working example:

\documentclass{article}

\begin{document}

\begin{itemize}

\item 1. First
\item 2. Second
\item 3a. Third
\item 3b. Third
\item 4. Fourth

\end{itemize}

\end{document}

The result:

How can I horizontally align the text (First, Second, Third etc)?

-
this is a hack, but ..., after 1. and the other single digit entries, insert \phantom{b}\  -- you do need that final space. –  barbara beeton Apr 11 at 18:53
Is it OK to assume that the first-level enumerate environment will never have more than nine items? –  Mico Apr 11 at 19:13

## 3 Answers

You can define a new command \myitem

\newcommand{\myitem}[1]{%
\item\makebox[1.5em][l]{#1}%
}

and use it instead of \item in these particular itemizes

MWE:

\documentclass{article}

\newcommand{\myitem}[1]{%
\item\makebox[1.5em][l]{#1}%
}

\begin{document}

\begin{itemize}

\myitem{1.} First
\myitem{2.} Second
\myitem{3a.} Third
\myitem{3b.} Third
\myitem{4.} Fourth

\end{itemize}

\end{document}

-

It seems more appropriate to use an enumerate environment for numeral enumeration. However, there are a number of ways to achieve this:

\documentclass{article}
\usepackage{calc}% http://ctan.org/pkg/calc
\begin{document}

\begin{itemize}
\item 1.\phantom{a}\ First
\item \leavevmode\rlap{2.}\phantom{3a.} Second
\item 3a.\ Third
\item 3b.\ Third
\item \makebox[\widthof{3a.}][l]{4.} Fourth
\end{itemize}

\end{document}

\rlap provides a right overlap, while a LaTeX2e approach would be to use \makebox (together with the help of calc's \widthof). There is a mild difference in the alignment around 3a. and 3b.. So, depending on your level of detail, you may want to use some of the other overlapping techniques in the sub-enumeration as well.

The use of \  is so that the period . doesn't produce a sentence-ending space.

The adjustbox package can provide similar "lapping" support.

-

I would recommend using an automated approach using an enumerate environment to add the numbering:

## Code:

\documentclass{article}
\usepackage{enumitem}

\newcommand*{\MakeBox}[1]{\makebox[1.5em][r]{#1}}%

\begin{document}

\begin{enumerate}[label={\textbullet\protect\MakeBox{\theenumi.}}]

\item First
\item Second
\begin{enumerate}[label={\textbullet\protect\MakeBox{\arabic{enumi}\alph*.}}, align=right, before={\stepcounter{enumi}}, leftmargin=0pt]
\item Third
\item Third
\end{enumerate}
\item Fourth

\end{enumerate}

\end{document}

-


text:
# R/spicefp.R In SpiceFP: Sparse Method to Identify Joint Effects of Functional Predictors

#### Documented in spicefp

#' spicefp
#'
#' This function is used to implement the spiceFP approach. This approach
#' transforms 2 (by default) or 3 functional predictors into candidate explonatory matrices in order to
#' identify joint classes of influence. It can take functional predictors and
#' partitioning functions as inputs in order to create candidate matrices to be
#' evaluated. The user can choose among the existing partitioning functions (as logbreaks) or provide
#' his own partitioning functions specific to the functional predictors under consideration. The user
#' can also directly provide candidate matrices already constructed as desired.
#'
#' @param y a numerical vector. Contains the dependent variable. This vector will
#' be used as response variable in the construction of models involving each
#' candidate matrix.
#' @param fp1 a numerical matrix with in columns observations of one statistical
#' individual to partition. Each column corresponds to the functional
#' predictor observation for one statistical individual. The order of the
#' statistical individuals is the same as in fp2. It is assumed that no
#' data are missing and that all functional predictors are observed on an equidistant (time) scale.
#' @param fp2 a numerical matrix with the same number of columns and rows as fp1.
#' Columns are also observations.
#' The order of the statistical individuals is the same as in fp1.
#' @param fp3 NULL by default. A numerical matrix with the same number of columns and
#' rows as fp1 and fp2.
#' The order of the statistical individuals is the same as in fp1 and fp2.
#' @param fun1 a function object with 2 arguments. First argument is fp1 and
#' the second is a list of parameters that will help to partition fp1, such as
#' the number of class intervals, etc. For example using the logbreaks function,
#' the list of parameters is equivalent to list(alpha, J). All the arguments
#' to be varied for the creation of different candidate matrices must
#' be stored in the parameter list. The other arguments must be set by default.
#' @param fun2 a function object with 2 arguments. First argument is fp2 and the
#'  second is a list of parameters.
#' @param fun3 NULL by default. Same as fun1 and fun2, a function with 2
#' arguments fp3 and a list of parameters.
#' @param parlists a list of 2 elements when fp3 and fun3 are equal to NULL or of 3
#' elements when fp3 and fun3 are provided.
#' All the elements of parlists are lists that have the same length. Each list
#' contains all the lists of parameters that have to be used to create different
#'  candidates. The first element of parlists concerns the first functional
#'  predictor fp1, the second element is relative to fp2 and the third to fp3.
#' @param xcentering TRUE by default. Defined whether or not the variables in
#' the new candidate matrices should be centered.
#' @param xscaling FALSE by default. Defined whether or not the variables in
#' the candidate matrices should be scaled.
#' @param candmatrices NULL by default. List. Output of the "candidates"
#' function. The spiceFP dimension is its first element. The second contains
#' many lists of one candidate matrix and related vector with index and numbers
#' of class intervals used per predictor. The other elements of the lists are
#' the inputs of "candidates" function. If the user does not need the
#' "candidates" function for the creation of candmatrices, it is possible to
#' build a list while making sure that it respects the same structure as well as
#' the names of the outputs of the "candidates" function. In this case, only the
#' first two elements of the list are essential: spicefp.dimension and
#' candidates. The remaining elements can be NULL.
#' @param K number of iterations of the spiceFP approach. Equal to 2 by default.
#' @param criterion character. One of "AIC_", "BIC_", "Cp_". The
#' criterion to be used in each iteration in order to identify the best
#' candidate matrix and to estimate the regulation parameters. This criterion
#' is used to perform model selection as well as variable selection.
#' @param penratios a numeric vector with values greater than or equal to 0. It
#' represents the ratio between the regularization parameters of parsimony and
#' fusion. When penratios=0, it corresponds to the pure fusion. The higher its
#' value, the more parsimonious the model is.
#' @param nknots integer. For one value in penratios vector, it represents the
#' number of models that will be constructed for each candidate matrix. It is
#' the argument "nlam" of \code{\link[genlasso:coef.genlasso]{coef.genlasso}}
#' function.
#' This argument can be also NULL. In this case, the argument appropriate.df
#' must be provided.
#' @param appropriate.df (appropriate degree of freedom) NULL by default. When used,
#' nknots must be NULL. It is
#' the argument "df" of \code{\link[genlasso:coef.genlasso]{coef.genlasso}}
#' function. When the user has a prior idea of the number of zones of influence that the
#' solution could contain, it is advisable to provide appropriate.df,
#' a vector of appropriate degrees of freedom.
#' appropriate.df is a numerical vector with values greater than or equal to 1.
#' The degree of freedom of generalized fused Lasso models is equal to the number of
#' connected components. A connected component gives information on a group of
#' non-zero coefficients sharing the same value and connected by a contiguity
#' matrix. More simply, it can be interpreted as a group of coefficients that
#' have a unique influence.
#'
#' @param penfun function with 2 arguments (dim1, dim2) when dealing with 2
#' dimensional spiceFP, or with 3 arguments (dim1, dim2, dim3) when dealing with
#' 3 dimensional spiceFP. The argument order in the penalty function is
#' associated with the order of numbers of class intervals used per predictor
#' in the second element of candmatrices argument. NULL by default.
#' When penfun=NULL, getD2dSparse of genlasso or getD3dSparse is used according
#' to the dimension of spiceFP.
#' @param dim.finemesh numeric vector of length 2 or 3. This vector informs
#' about the dimension of the fine-mesh arrays (or matrices) that will be used
#' for the visualization of the sum of the coefficients selected at different
#' iterations.
#' @param file_name character vector. Of length K, it contains the list of
#' names that will be used to name the files containing informations on the
#' candidate matrix models
#' @param ncores numbers of cores that will be used for parallel computation. By default, it is equal
#' to detectCores()-1.
#' @param write.external.file logical. indicates whether the result table related
#' to each iteration should be written as a file (txt) in your working directory.
#' It is recommended to use write.external.file=TRUE when evaluating a large
#' number of candidate matrices (more than 100) in order to keep memory available.
#'
#'
#' @details  Three main steps are involved to implement spiceFP: transformation
#' of functional predictors, creation of a graph of contiguity constraints and
#' identification of the best class intervals and related regression
#' coefficients.
#'
#' @return Returns a list with:
#'
#' \describe{
#' \item{Candidate.Matrices}{a list with candidate matrices and their
#' characteristics. same as candmatrices if it has been provided.}
#' \item{Evaluations}{List of length less than or equal to K. Each element of
#' the list contains information about an iteration. Contains the results
#' related to the evaluation of the candidate matrices. These include the name
#' of the file where the model information is stored, the best candidate matrix
#' and related coefficients, the partition vector that indexes it, the \eqn{X \beta}
#' estimation, the residuals, etc.}
#' \item{coef.NA}{List of length less than or equal to K. For each iteration,
#' it contains the coefficient vector where the coefficient value of
#' never-observed joint modalities is NA}
#' \item{coef.NA.finemeshed}{List of length less than or equal to K. For each
#' iteration, the coefficient vector is transformed into fine-mesh array or
#' matrix allowing arithmetic operations to be performed between coefficients
#' coming from different partitions}
#' \item{spicefp.coef}{fine-mesh array or matrix. Sum of the coefficients
#' selected at all iterations}
#' }
#'
#'
#' @importFrom doParallel registerDoParallel
#' @importFrom foreach foreach %dopar%
#' @importFrom stringr str_c
#' @importFrom parallel detectCores
#' @importFrom purrr pmap map
#' @importFrom genlasso fusedlasso fusedlasso2d coef.genlasso
#' @importFrom stats sd lm var
#'
#' @examples
#'\donttest{
#' ##linbreaks: a function allowing to obtain breaks linearly
#' linbreaks<-function(x,n){
#'     sort(round(seq(trunc(min(x)),
#'                 ceiling(max(x)+0.001),
#'                 length.out =unlist(n)+1),
#'             1)
#'         )
#' }
#'
#' # In this example, we will evaluate 2 candidates with 14 temperature
#' # classes and 15 irradiance classes. The irradiance breaks are obtained
#' # according to a log scale (logbreaks function) with different alpha
#' # parameters for each candidate (0.005, 0.01).
#' ## Data and inputs
#' tpr.nclass=14
#' irdc.nclass=15
#' irdc.alpha=c(0.005, 0.01)
#' p2<-expand.grid(tpr.nclass, irdc.alpha, irdc.nclass)
#' parlist.tpr<-split(p2[,1], seq(nrow(p2)))
#' parlist.irdc<-split(p2[,2:3], seq(nrow(p2)))
#' parlist.irdc<-lapply(
#'    parlist.irdc,function(x){
#'    list(x[[1]],x[[2]])}
#' )
#' m.tpr <- as.matrix(Temperature[,-c(1)])
#'
#' # For the constructed models, only two regularization parameter ratios
#' # penratios=c(1/25,5) are used. In a real case, we will have to evaluate
#' # more candidates and regularization parameters ratio.
#' start_time_sp <- Sys.time()
#' ex_sp<-spicefp(y=FerariIndex_Difference$fi_dif, #' fp1=m.irdc, #' fp2=m.tpr, #' fun1=logbreaks, #' fun2=linbreaks, #' parlists=list(parlist.irdc, #' parlist.tpr), #' penratios=c(1/25,5), #' appropriate.df=NULL, #' nknots = 100, #' ncores =2, #' write.external.file=FALSE) #' #' duration_sp <- Sys.time() - start_time_sp #' # View(ex_sp$Evaluations[[1]]$Evaluation.results$evaluation.result)
#' # View(ex_sp$Evaluations[[2]]$Evaluation.results$evaluation.result) #' # Visualization of the coefficients #' g<-ex_sp$spicefp.coef
#' g.x<-as.numeric(rownames(g))
#' g.y<-as.numeric(colnames(g))
#'
#' library(fields)
#' plot(c(10,2000),c(15,45),type= "n", axes = FALSE,
#'      xlab = "Irradiance (mmol/m²/s - Logarithmic scale)",
#'      ylab = "Temperature (°C)",log = "x")
#' rect(min(g.x),min(g.y),max(g.x),max(g.y), col="black", border=NA)
#' image.plot(g.x,g.y,g, horizontal = FALSE,
#'            col=designer.colors(256, c("blue","white","red")),
#' axis(1) ; axis(2)
#'
#' closeAllConnections()
#'
#'}
#'
#' @export

spicefp<-function(y,
fp1, fp2, fp3=NULL,
fun1, fun2, fun3=NULL,
parlists,
xcentering=TRUE,
xscaling=FALSE,
candmatrices=NULL, K=2,
criterion="AIC_",
penratios=c(1/10,1/5,1/2,1,2,5,10),
nknots=50,
appropriate.df=NULL,
penfun=NULL,
dim.finemesh=c(1000,1000),
file_name=paste0("parametertable",1:2),
ncores=parallel::detectCores()-1,
write.external.file=TRUE){

try(if( !is.logical(write.external.file)) stop("write.external.file must be logical "))

try(if(sum(criterion %in% c("AIC_","BIC_","AICc_","Cp_"))!=1)
stop("Only one criterion has to be chosen among AIC_, BIC_, AICc_, Cp_"))
try(if(length(unique(file_name))!=K)
stop("Different filenames must be provided for iterations of spicefp.
length(unique(file_name)) should be equal to K"))
if (is.null(candmatrices)) {
# stop candidates
try(if(is.null(fp3)+is.null(fun3)==1)
stop("fp3 and fun3 must be provided both or equal to NULL if only 2 functional predictors are used"))
try(if(is.null(fp3) & length(parlists)!=2)
stop("The length of parlists should be equivalent to the number of functional predictors provided"))
try(if(ncol(fp1)!=ncol(fp2))
stop("fp1 and fp2 must have the same number of columns "))
}
if(!is.null(candmatrices)) {
try(if(!is.list(candmatrices)|sum(names(candmatrices)==c("spicefp.dimension","candidates","fp1","fp2","fp3","fun1","fun2",
"fun3","parlists","xcentering","xscaling"))!=11)
stop("candmatrices x must be a list containing 11 items named spicefp.dimension,candidates,fp1,fp2,fp3,fun1,fun2,
fun3,parlists,xcentering,xscaling. It corresponds to an object of class xxxxxx"))
}
try(if(is.null(appropriate.df)+is.null(nknots)>1)
stop("Only one of appropriate.df or nknots must be provided. The other must be NULL"))

if(is.null(candmatrices)){try(if((is.null(fp3) & length(dim.finemesh)!=2) | (!is.null(fp3) & length(dim.finemesh)!=3))
stop("The length of the dim.finemesh vector must correspond to the spicefp dimension used"))} else  {
try(if(length(dim.finemesh)!=candmatrices$spicefp.dimension) stop("The length of the dim.finemesh vector must correspond to the spicefp dimension used"))} try(if(isTRUE(write.external.file) & !is.character(file_name)) stop("If write.external.file=TRUE then filenames (one per iteration) must be provided")) # Candidate construction candmatrices <- if(is.null(candmatrices)){ candidates(fp1,fp2,fp3, fun1,fun2,fun3, parlists, ncores, xcentering, xscaling) } else { candmatrices } wef<-write.external.file ; rm(write.external.file) response_var<-y tokeep <-list() iter_=1 check_coef_zeros <- FALSE while(iter_ <= K & check_coef_zeros==FALSE){ rescd <- evaluate.candidates(candmatrices, y=response_var, penratios, nknots, appropriate.df, ncores, penfun, file_name[iter_], write.external.file=wef) rescd_file<-if(wef){ read.table(rescd$evaluation.result, header = TRUE)} else {
rescd$evaluation.result } try(if(sum(rescd_file[,"Slope_"]>0.05 & rescd_file[,"Slope_"]<0.95)==0) stop("At iteration ", iter_ ,", none of the candidate matrices produces a model with a slope between 0.05 and 0.95")) thecand.par<-as.list(rescd_file[rescd_file[,"Slope_"]>0.05 & rescd_file[,"Slope_"]<0.95, ][which.min(rescd_file[rescd_file[,"Slope_"]>0.05 & rescd_file[,"Slope_"]<0.95, ][,criterion]),]) rm(rescd_file) ## Coef for thecand.par thecand = candmatrices$candidates[[thecand.par$Candidate_id]][[1]] dim.thecand = candmatrices$candidates[[thecand.par$Candidate_id]][[2]] suppressWarnings({ thecand.mdl<- if (candmatrices$spicefp.dimension==2){
fusedlasso(y=response_var,
X=thecand,
D=rescd$penalty.function(dim.thecand[2], dim.thecand[3]), gamma = thecand.par$Pen_ratio)
} else {
fusedlasso(y=response_var,
X=thecand,
D=rescd$penalty.function(dim.thecand[2], dim.thecand[3], dim.thecand[4]), gamma = thecand.par$Pen_ratio)
}
})

thecand.coef <-coef.genlasso(thecand.mdl,lambda=thecand.par$PenPar_fusion) rownames(thecand.coef$beta) <- colnames(thecand)
xbeta<- thecand%*%thecand.coef$beta thecand.residuals<-response_var-xbeta tokeep[[iter_]]<-list(Evaluation.results=rescd, thecandidate.parameters=thecand.par, thecandidate=thecand, dim.thecandidate=dim.thecand, thecandidate.model=thecand.mdl, thecandidate.coefficients=thecand.coef, XBeta=xbeta, thecandidate.residuals=thecand.residuals) ## Check and next if(isTRUE(all.equal(c(thecand.coef$beta), rep(0,nrow( thecand.coef$beta ))))) print(paste("Best coefficients chosen by", criterion,"at iteration",iter_,"are all equal to 0.", "Stop after iteration",iter_)) if(isTRUE(all.equal(c(thecand.coef$beta), rep(0,nrow( thecand.coef$beta ))))){ check_coef_zeros=T } iter_=iter_+1 response_var<-thecand.residuals } ## Treat the coef for sum coef.NA=list() coef.NA.finemeshed=list() for (g in 1:length(tokeep)) { amop<-tokeep[[g]] nrowcand<-nrow(amop$thecandidate)
tobeNA<- unlist(map(1:ncol(amop$thecandidate), function(x){ if(sum(amop$thecandidate[,x]==0)==nrowcand){x}}))
coefcand<-amop$thecandidate.coefficients$beta
coefcand[tobeNA,1]<-NA
coefcand.finemeshed <- if(candmatrices$spicefp.dimension==2){ finemeshed2d(coefcand, n.breaks1=dim.finemesh[1], n.breaks2=dim.finemesh[2])$finemeshed.matrix} else {
finemeshed3d(coefcand,
n.breaks1=dim.finemesh[1],
n.breaks2=dim.finemesh[2],
n.breaks3=dim.finemesh[3])$finemeshed.array} coef.NA[[g]]<-coefcand coef.NA.finemeshed[[g]]<- coefcand.finemeshed } spicefp_coefsum <- function(x){ ifelse(setequal(x, rep(NA,length(x))), NA, sum(x,na.rm = T)) } if(candmatrices$spicefp.dimension==2){
acoef<-array(NA,c(length(coef.NA.finemeshed),dim.finemesh[1],dim.finemesh[2]))
for (i in 1:length(coef.NA.finemeshed)) {
acoef[i,,]<-coef.NA.finemeshed[[i]]
}
spicefp.coef= apply(acoef, c(2,3),FUN = function(x){spicefp_coefsum(x)})
dimnames(spicefp.coef)=dimnames(coef.NA.finemeshed[[1]])} else {
acoef<-array(NA,c(length(coef.NA.finemeshed),
dim.finemesh[1],
dim.finemesh[2],
dim.finemesh[3]))
for (i in 1:length(coef.NA.finemeshed)) {
acoef[i,,,]<-coef.NA.finemeshed[[i]]
}
spicefp.coef= apply(acoef, c(2,3,4),FUN = function(x){spicefp_coefsum(x)})
dimnames(spicefp.coef)=dimnames(coef.NA.finemeshed[[1]])}

## Global result
return(list(Candidate.Matrices=candmatrices,
Evaluations=tokeep,
coef.NA=coef.NA,
coef.NA.finemeshed=coef.NA.finemeshed,
spicefp.coef=spicefp.coef))
}

## Try the SpiceFP package in your browser

Any scripts or data that you put into this service are public.

SpiceFP documentation built on Sept. 15, 2021, 9:07 a.m.


text:
# Calling NAG Library Routines from Scilab - Example 4

## 3.4. Example 4

### Evaluating a multi-dimensional integral (d01wcc)

In this example we will demonstrate two ways of providing a user-defined function to a NAG routine in Scilab. The function we shall use is d01wcc, which evaluates a multidimensional integral of the function f(x) where x is a multidimensional vector. The user-supplied function will evaluate this function at any x.

The two ways of providing this function to the NAG routine are either to write the function in C, or to write it as a Scilab function. Both approaches have their advantages and disadvantages. The C code is a lot easier to write and requires fewer lines, but in order to change it, each time you would need to rebuild the interface. Writing the function in Scilab code is slightly harder, as you need to write a second interface for the Scilab function in C, but it does give you the flexibility to change your function without having to rebuild the Scilab-NAG interface.

1. ### Function prototype from the NAG C Library Manual

According to the C Library Manual, the prototype for function d01wcc looks like this:

    #include <nag.h>
#include <nagd02.h>

const double a[], const double b[], Integer *minpts, Integer maxpts, double eps,
double *finval, double *acc, Nag_User *comm, NagError *fail)

Note that the user-supplied function is named f, and is of type double. For more information about the argument list of nag_multid_quad_adapt_1, refer to the NAG C Library manual.

1. ### C Interface

Here is the source code of our interface written in C nag_intext4a.c with the user-supplied function also written in C:

    /* Example 4a
=========

Shows how to implement a Scilab wrapper to a NAG function using a User
defined function written in C */

#include "stack-c.h"
#undef Complex
#define Complex NagComplex
#include <nag.h>
#include <nag_stdlib.h>

#include <nagd01.h>
#define SciComplex doublecomplex

#include "stdio.h"

static double f(Integer n, double z[], Nag_User *comm);

int nag_intext4a(char *fname)
{
// to call this function in scilab use:
int m1,n1,l1;
int m2,n2,l2;
int m3,n3,l3;
int m4,n4,l4;
int m5,n5,l5;
int m6,n6,l6;
int l7;
int l8;
int l9;
int l10;

int n, i, min, max;

Integer ndim, minpts, maxpts;
double eps, finval, acc;
double *a=0, *b=0;
Nag_User comm;
NagError ifail;

// define minimum and maximum left and right hand arguments
int minlhs=1, minrhs=6, maxlhs=4, maxrhs=6;
Nbvars = 0;
CheckRhs(minrhs, maxrhs);
CheckLhs(minlhs,maxlhs);

// Initialise ifail
INIT_FAIL(ifail);

GetRhsVar(1, "i", &m1, &n1, &l1); // input ndim
GetRhsVar(2, "d", &m2, &n2, &l2); // input a
GetRhsVar(3, "d", &m3, &n3, &l3); // input b
GetRhsVar(4, "i", &m4, &n4, &l4); // input minpts
GetRhsVar(5, "i", &m5, &n5, &l5); // input maxpts
GetRhsVar(6, "d", &m6, &n6, &l6); // input eps

//******************************************************************
// Read in input variables from Scilab and convert to NAG variables
//******************************************************************

//*============*
if (m1!=1 || n1!=1)
{
sciprint("%s: Dimension should be 1x1 character for arg 1\r\n",fname);
Error(999); goto END;
}
else
{
n = *istk(l1);
ndim = (Integer) n;
}

// Allocate memory for a(ndim) and b(ndim)
//========================================
if ( !(a = NAG_ALLOC(n,double) ) ||
!(b = NAG_ALLOC(n,double)))
{
sciprint("Allocation failure\n");
Error(999); goto END;
}

//===========
if (m2!= n || n2!=1)
{
sciprint("%s: Dimension should be ndim by 1 for arg 2\r\n",fname);
sciprint("n = %i", n);
Error(999); goto END;
}
else
{
for (i=0; i<n; i++)
{
a[i] = *stk(l2+i);
}
}

//==============
if (m3!=n || n3!=1)
{
sciprint("%s: Dimension should be ndim by 1 for arg 3\r\n",fname);
Error(999); goto END;
}
else
{
for(i=0; i<n; i++)
{
b[i] = *stk(l3+i);
}
}

//============
if (m4!=1 || n4!=1)
{
sciprint("%s: Dimension should be 1x1 for arg 4\r\n",fname);
Error(999); goto END;
}
else
{
min = *istk(l4);
minpts = (Integer)min;
}

//============
if (m5!=1 || n5!=1)
{
sciprint("%s: Dimension should be 1x1 for arg 5\r\n",fname);
Error(999); goto END;
}
else
{
max = *istk(l5);
maxpts = (Integer) max;
}

//==============
if (m6!=1 || n6!=1)
{
sciprint("%s: Dimension should be 1x1 for arg 6\r\n",fname);
Error(999); goto END;
}
else
{
eps = *stk(l6);
}
//***********************************
// End of reading in input variables
//***********************************

//*====================================================
&comm, &ifail);

// Print NAG error message if routine fails
//=========================================
if (ifail.code != NE_NOERROR)
sciprint("Ifail message is %s", ifail.message);

// Create output variables for Scilab
//====================================
CreateVar(7, "d", &n1, &m1, &l7); // finval
CreateVar(8, "d", &n1, &m1, &l8); // acc
CreateVar(9, "i", &n1, &m1, &l9); // minpts
CreateVar(10, "i", &n1, &m1, &l10); // ifail

// Assign scilab output variables
//================================
*stk(l7) = finval;
*stk(l8) = acc;
*istk(l9) = (int) minpts;
*istk(l10) = (int) ifail.code;

// Assign order for output variables
//===================================
LhsVar(1) = 7;
LhsVar(2) = 8;
LhsVar(3) = 9;
LhsVar(4) = 10;

END:

// Deallocate memory
//===================
if (a) NAG_FREE(a);
if (b) NAG_FREE(b);

return(0);
}

static double f(Integer n, double z[], Nag_User *comm)
{
// Function describing the integrand for the NAG routine used above
//==================================================================
double tmp_pwr;
tmp_pwr = z[1]+1.0+z[3];
return z[0]*4.0*z[2]*z[2]*exp(z[0]*2.0*z[2])/(tmp_pwr*tmp_pwr);
}

• The input function for the NAG routine has a specific interface that has to be adhered to, which is documented in the NAG C Library manual [5]
• The interface doesn't take in the input arguments of the function as its arguments, but only takes fname the name of the function in Scilab. In this example we have chosen fname to be the same as the NAG C library name.
• The Scilab function GetRhsVar is used to retrieve the input arguments of fname
• CreateVar creates a variable into which the NAG output is placed
• The interface doesn't actually return the outputs of the function fname - instead the Scilab array LhsVar stores the variable numbers that are stored on Scilab's stack stk()
2. ### Compiling with the Builder Script

To compile and link this function, you can run the build script nag_builder4a.sce which links this interface to the function name fname that will be seen in Scilab. This must be run in the same directory as the interface file, and when run successfully, generates a loader script loader.sce that needs to be executed in order to dynamically link the newly created library into Scilab.

    exec nag_builder4a.sce

3. ### Calling the function

If the function built and the loader.sce file loaded the function without problem, then we can use the function within Scilab, either on the command-line or in a script. Here we run the example script, nag_example4a.sce.

    --> exec nag_example4a.sce
finval = 0.5753623

Tip: If you get any error messages in Scilab check the troubleshooting section.

3. ### Writing the user-supplied functions in Scilab

1. #### C Interface for the Scilab code

The file nag_intext4b.c contains similar code to the example above.

The only thing that is different between the two source codes nag_intext4a.c and nag_intext4b.c is that in the second one, instead of the user-supplied function being written in C, we have an interface to the user-supplied function written in Scilab. This is so that any variable types can be converted and so that the interface to the NAG routine is correct. Below we just show the part of the file containing the new interface, where we have named the function fSci.

    static double fSci(Integer n, double z[], Nag_User *comm);

static double fSci(Integer n, double z[], Nag_User *comm)
{
int mlhs, mrhs, ibegin;
int m_u, n_u, l_u;
int m_z, n_z, l_z;
int i;
double retval;
char name[] = "fUser";

// Put input variables for scilab function onto stack in the right order
// On input u (stored in space 7) is a dummy input to be overwritten with
// the output from fUser
// The inputs for fUser must be in the correct order in the last places of
// the stack.  They are then scrapped and the outputs are created starting
// at the same place as the first input variable.
m_u = 1;
n_u = 1;
m_z = (int) n;
n_z = 1;

CreateVar(7, "d", &m_u, &n_u, &l_u); // first input
CreateVar(8, "d", &m_z, &n_z, &l_z); // second input

// Write the input vector z onto the stack in the second input
for (i=0; i<m_z; i++)
{
*stk(l_z+i) = z[i];
}

ibegin = 7; // the first input variable on the stack
mlhs = 1;  // number of left hand side outputs to write onto the stack
mrhs = 2; // number of right hand side inputs to read off the stack
SciString(&ibegin,name,&mlhs,&mrhs); // Call scilab function 'name'

// output from scilab function is now on the stack so return the result
retval = *stk(l_u);

return retval;
}

• The User-defined function must still have the same interface as specified in the C Library manual and therefore we must write an interface to the Scilab function.
• The Scilab function fUser is not called directly, but through a Scilab function called SciString. The arguments for SciString are
• ibegin: the variable number of the first input argument. Note that the input arguments must be next to each other and in the order they are to passed into fUser with the first input argument being at ibegin.
• name: the name of the Scilab function.
• mlhs and mrhs: the number of input and output variables of the Scilab function. The input arguments are destroyed after use, and the output arguments are then put on the stack in order starting at ibegin.
2. ### User-supplied functions in Scilab

Source code of our Scilab function is in the example script nag_example4b.sce but could equally have been in a separate script. The source also demonstrates two different ways of writing the function.

    deff('u=fUser(uu,z)',
['tmp_pwr = z(2)+1.0+z(4)';
'u = z(1)*4.0*z(3)*z(3)*exp(z(1)*2.0*z(3))/(tmp_pwr*tmp_pwr)']);

3. ### Compiling with the Builder Script

To compile and link this function, you need to run the build script nag_builder4b.sce which links this interface to the function name fname that will be seen in Scilab. This must be run in the same directory as the interface file, and when run successfully, generates a loader script loader.sce that needs to be executed in order to dynamically link the newly created library into Scilab.

exec nag_builder4b.sce

      --> exec nag_example4b.sce


text:
# Relationship between these two probability mass functions

If I have two different discrete distributions of random variables X and Y, such that their probability mass functions are related as follows:

$P(X=x_i) = \lambda\frac{P (Y=x_i)}{x_i}$

what can I infer from this equation? Any observations or interesting properties that you see based on this relation?

What if,
P($X=x_i$) = $\lambda\sqrt{\frac{P (Y=x_i)}{x_i} }$

In both cases, $\lambda$ is a constant.

-
It would be clearer to write this with two random variables and one measure $P$; i.e. $P(X = x_i) = \lambda P(Y=x_i)/x_i$. Also, "what can I infer" is very vague; what sorts of things would you like to infer? –  Nate Eldredge Nov 4 '10 at 22:58
@Nate:Any observation will do. I am applying a polynomial decay on a pmf to get another pmf. Something like 'A uniform Y will result in a geometric distribution(?) for X.' –  SkypeMeSM Nov 4 '10 at 23:35
Also on MathOverflow: mathoverflow.net/questions/44902/… (repaired clipboard failure) –  Nate Eldredge Nov 5 '10 at 3:06

if you rewrite the equation as $P(Y = x) = xP(X = x)/\lambda$, then the distribution of $Y$ is called the length-biased distribution for $X$. it arises, for example, if one has a bunch of sticks in a bag and reaches in and selects one at random - where the probability a particular stick is selected is proportional to its length. if the lengths of the sticks are realizations of the random variable $X$, the distribution of the length of the selected stick is that of $Y$.

-

I don't have enough reputation yet, or this would be a comment.

A random variable can have only one probability mass function, so it is not clear what you are asking. Where do these equations come from?

-
I am considering two different systems which operate on the same random variable X and have two different pmfs' u and v. Edited the question to reflect this. –  SkypeMeSM Nov 4 '10 at 22:05

First, you don't have the freedom to choose $\lambda$. (Maybe you already realize that?) In order for $P(X = x_i)$ to be a true probability distribution, you must have $$\lambda = \left(\sum_{x_i} \frac{P(Y = x_i)}{x_i}\right)^{-1}.$$

Given that, if $Y$ is geometric with success probability $p$, and $P(X = x_i) = \lambda P(Y = x_i)/x_i$ then $X$ does have the logarithmic distribution with parameter $q = 1-p$.

Here's why: $P(Y = k) = (1-p)^{k-1}p$, for $k = 1, 2, \ldots$. Thus $$P(X = k) = \lambda \frac{P(Y = k)}{k} = \lambda \frac{q^{k-1} p}{k} = \frac{\lambda p}{q} \frac{q^k}{k}.$$ The expression on the right is the logarithmic probability mass function with $\lambda = \frac{-q}{p \ln p}$.

-


text:
## Hardness of Computing Individual Bits for One-way Functions on Elliptic Curves

We prove that if one can predict any of the bits of the input to an elliptic curve based one-way function over a finite field, then we can invert the function. In particular, our result implies that if one can predict any of the bits of the input to a classical pairing-based one-way function with non-negligible advantage over a random guess then one can efficiently invert this function and thus, solve the Fixed Argument Pairing Inversion problem (FAPI-1/FAPI-2). The latter has implications on the security of various pairing-based schemes such as the identity-based encryption scheme of Boneh–Franklin, Hess’ identity-based signature scheme, as well as Joux’s three-party one-round key agreement protocol. Moreover, if one can solve FAPI-1 and FAPI-2 in polynomial time then one can solve the Computational Diffie--Hellman problem (CDH) in polynomial time. Our result implies that all the bits of the functions defined above are hard-to-compute assuming these functions are one-way. The argument is based on a list-decoding technique via discrete Fourier transforms due to Akavia--Goldwasser–Safra as well as an idea due to Boneh–Shparlinski.

Editor(s):
Safavi-Naini, Reihaneh
Canetti, Ran
Published in:
Proceedings of the 32nd Annual Cryptology Conference, 7417, 832-849
Presented at:
CRYPTO 2012 - 32nd Annual Cryptology Conference, Santa Barbara, California, USA, August 19-23, 2012
Year:
2012
Publisher:
Springer
Keywords:
Laboratories:


text:
Waiting for answer This question has not been answered yet. You can hire a professional tutor to get the answer.

QUESTION

# Calculate the reaction quotient, Q, for mixing 100ml of 0.00028M "Pb"("NO"_3)_2 with 200ml of 0.0012M "NaCl". The K_(sp) for "PbCl"_2 is 1.9 xx10^(-4). Does a precipitate form?

No, a precipitate does not form.

Your goal here is to determine if the concentrations of lead(II) cations, "Pb"^(2+), and chloride anions, "Cl"^(-), are high enough to allow for the formation of lead(II) chloride, "PbCl"_2, an insoluble solid that precipitates out of solution.

Lead(II) nitrate, "Pb"("NO"_3)_2, dissociates in a 1:1 mole ratio to form lead(II) cations

"Pb"("NO"_ 3)_ (color(red)(2)(aq)) -> "Pb"_ ((aq))^(2+) + color(red)(2)"NO"_ (3(aq))^(-)

For every mole of lead(II) nitrate present in the solution, you will have one mole of lead(II) cations. Use the and volume of the solution to determine how many moles of lead(II) cations you have

color(purple)(|bar(ul(color(white)(a/a)color(black)(c = n_"solute"/V_"solution" implies n_"solute" = c * V_"solution")color(white)(a/a)|)))

n_(Pb^(2+)) = "0.00028 mol" color(red)(cancel(color(black)("L"^(-1)))) * overbrace(100 * 10^(-3)color(red)(cancel(color(black)("L"))))^(color(blue)("volume in liters"))

n_(Pb^(2+)) = 2.8 * 10^(-5)"moles Pb"^(2+)

Do the same for the sodium chloride solution, "NaCl"

"NaCl"_ ((aq)) -> "Na"_ ((aq))^(+) + "Cl"_((aq))^(-)

You will have

n_(Cl^(-)) = "0.0012 mol" color(red)(cancel(color(black)("L"^(-1)))) * 200 * 10^(-3)color(red)(cancel(color(black)("L")))

n_(Cl^(-)) = 2.4 * 10^(-4)"moles Cl"^(-)

Now, the volume of the resulting solution (assuming volumes are additive) will be

V_ "total" = V_(Pb(NO_ 3)_ 2) + V_(NaCl)

V_"total" = "100 mL" + "200 mL" = "300 mL"

The concentrations of the lead(II) cations and of the chloride anions (not yet incorporating "PbCl"_2) in the resulting solution will be

["Pb"^(2+)] = (2.8 * 10^(-5)"moles")/(300 * 10^(-3)"L") = 9.33 * 10^(-5)"M"

["Cl"^(-)] = (2.4 * 10^(-4)"moles")/(300 * 10^(-3)"L") = 8.0 * 10^(-4)"M"

Now, lead(II) chloride dissociates according to the following equilibrium reaction

"PbCl"_ (color(red)(2)(s)) rightleftharpoons "Pb"_ ((aq))^(2+) + color(red)(2)"Cl"_((aq))^(-)

Its solubility product constant, K_(sp), is equal to

K_(sp) = ["Pb"^(2+)] * ["Cl"^(-)]^color(red)(2)

Keep in mind that the solubility product constant is calculated using equilibrium concentrations. In order to determine whether or not a precipitate will form, you need to use the initial concentrations of the two ions.

This is what the reaction quotient, Q_(sp), is used for

Q_(sp) = ["Pb"^(2+)]_"initial" * ["Cl"^(-)]_"initial"^color(red)(2)

A precipitate will form if

color(purple)(|bar(ul(color(white)(a/a)color(black)(Q_(sp) > K_(sp))color(white)(a/a)|)))

Plug in your values to find

Q_(sp) = 9.33 * 10^(-5) * (8.0 * 10^(-4))^color(red)(2)

Q_(sp) = 6.0 * 10^(-11) " << " stackrel(K_(sp))overbrace(1.9xx10^(-4))

Since Q_(sp) < K_(sp), and actually quite significantly smaller, a precipitate will not be formed when these two are mixed.

color(green)(|bar(ul(color(white)(a/a)color(black)(Q_(sp) < K_(sp) implies "NO precipitate")color(white)(a/a)|)))


text:
user54412
12:11 AM
@BrandonEnright looks like you got it just in time :)

user54412
and without my upvotes even ;)

user54412
... i keep forgetting what UT it is around here

12:43 AM
@ChrisWhite yeah I made it by ~4 min. Lucky.

3 hours later…
3:40 AM
0

The facts: No antigravity has been observed ... and no DE has been observed also. It was discovered (measured) in 1998 an acceleration in the expansion of space. The mass (M) is the source of gravity (G). Every force has a source of it, charges, masses,... . The source of antigravity (AG) is ...

Whoa. That is one doozey of an answer.

2 hours later…
6:06 AM
@BrandonEnright yeah. Awesome answer btw :)
@BrandonEnright I have the ability to shift sleep schedules with relative ease. So when I'm at home, I sleep on time because everyone else does. When in the institute during vacations, I sleep late and wake up late because my friends do too. So...yeah, screwed up sleep schedule :P
But quite common amongst students
@BrandonEnright First time?
It's been a while since I rep capped on Phys.
It used to be quite hard a year a go, probably could do it easily now if I have the time

1 hour later…
7:21 AM
3

According to this Wikipedia page, Mendeleev originally believed that the inert gases belonged in Group 0 (to the left of the alkali metals). Thus, helium would be placed in the second period to the left of lithium and argon would be placed in the third period to the left of sodium. Apparently thi...

7 hours later…
1:54 PM
Hey all, I posted a question here, and I think I know only one of the answers (what we are solving for), should I edit the message? Or post it as a reply to myself?
http://physics.stackexchange.com/questions/65023/explain-the-fopplvon-karman-equations

2 hours later…
3:48 PM
@GuySoft Self answers are not only allowed but encouraged.

4:02 PM
@GuySoft Also if you get a good answer you must add it and also accept it!
like :
1

Now I have got a method to get it directly.And again it came out to be an easy problem. The answer comes out to be $$I= \dfrac{m}{12}(a^2+c^2)$$ See if we add another such plate along it's side $AC$ , then it comes out to be a parallelogram plate , whose MOI is known, same as rectangular plate ...

4:27 PM
@007 I'd love to accept it once I get one :)

5:14 PM
@ManishEarth Thanks! I'm always amazed at how much detail must be left out to provide a casual answer.

@BrandonEnright ?
ah

I never knew there were so many quack theories out there until I started searching for topics based on questions on this site. The amount of energy going into crazy stuff is really staggering.

5:49 PM
Anyone here?

6:17 PM
@BrandonEnright same here :s
@Alyosha me. though i need to sleep soon, gotta wake up early tomorrow

1 hour later…
7:20 PM
@Ϛѓăʑɏ why do you want to close this perfectly good and even nice theoretical physics question? Why do you think anything is wrong with it? Please dont closevote (perfectly good) interesting theoretical physics questions, if you dont understand the physics well enough they are asking about to judge their validity :-(. Just leave theme alone or skip them in the review queues ... BTW cant you change your name back? You are no longer pingable...
This is a valid question that shoul not be closed either, so maybe you are a bit too generous with your closevotes sometimes ... ;-)? Please dont try to kill valid questions like this, if you are not sure about the physics they are asking about just leave them alone :-/
Darn @CrazyBudy this is a perfectly valid theoretical physics question too, WTF are you doing? Why are you closevoting perfectly good and even nice questions like this?


text:
# The Universe of Discourse

Sat, 28 Aug 2021

The most important combinator in combinatory logic is the !!S!! combinator, defined simply:

$$S x y z ⇒ (x z)(y z)$$

or in !!\lambda!!-calculus terms:

$$S = \lambda x y z. (x z)(y z).$$

A wonderful theorem states that any !!\lambda!!-expression with no free variables can be converted into a combinator expression that contains only the combinators !!S, K,!! and !!I!!, where !!S!! is really the only interesting one of the three, !!I!! being merely the identity function, and !!K!! a constructor of constant functions:

\begin{align} I x & = x \\ K x y & = x \\ \end{align}

In fact one can get along without !!I!! since !!S K K = I!!.

A not-too-infrequently-asked question is why the three combinators are named as they are. The !!I!! is an identity function and pretty obvious stands for “identity”.

Similarly the !!K!! constructs constant functions: !!K x!! is the combinator which ignores its argument and yields !!x!!. So it's not hard to imagine that !!K!! is short for Konstant, which is German for “constant”; no mystery there.

But why !!S!!? People typically guess that it stands for “substitution”, the idea being that if you have some application $$A\,B$$ then !!S!! allows one to substitute some term !!T!! for a free variable !!v!! in both !!A!! and !!B!! prior to the application:

$$S\, A\, B\, T = A[v/T]\, B[v/T].$$

Although this seems plausible, it's not correct.

Combinatory logic was introduced in a 1924 paper of Moses Schönfinkel. In it, he defines a family of combinators including the standard !!S!!, !!K!!, and !!I!!; he shows that only !!S!! and !!K!! are required. His initial set of combinators comprises the following:

$$\begin{array}{cllrl} I & \textit{Identitätsfunktion} & \text{“identity function”}& I\,x =& x \\ C & \textit{Konstanzfunktion} & \text{“constancy function”} & C\,x\,y =& x \\ T & \textit{Vertauschungsfunktion} & \text{“swap function”} & T\,x\,y\,z=& x\,z\,y \\ Z & \textit{Zusammensetzungsfunktion} & \text{“composition function”} & Z\,x\,y\,z=& x\,(y\,z) \\ S & \textit{Verschmelzungsfunktion} & \text{“fusion function”} & S\,x\,y\,z=& x\,z\,(y\,z) \end{array}$$

(Schönfinkel also had combinators representing logical operations (one corresponding to the Sheffer stroke, which had been discovered in 1913), and to quantification, but those don't concern us right now.)

!!T!! and !!Z!! are now usually called !!C!! and !!B!!. These names probably originated in Curry's Grundlagen der kombinatorischen Logik (1930). Curry 1930 is probably also the origin of the change from !!C!! to !!K!!. I have no idea why Schönfinkel chose to abbreviate Konstanzfunktion as !!C!! instead of !!K!!. Curry notes that for !!I, K, B, C, S!! Schönfinkel has !!I, C, Z, T, S!!, but does not explain his changes. In Curry and Feys’ influential 1958 book on combinatory logic, the !!B!! and !!C!! combinators given names that are are literal translations of Schönfinkel's: “elementary permutator” and “elementary compositor”.

Returning to the !!S!! combinator, one sees that its German name in Schönfinkel's paper, Verschmelzungsfunktion, begins with the letter V, but so does Vertauschungsfunktion, so abbreviating either with V would have been ambiguous. Schönfinkel instead chose to abbreviate Verschmelzungsfunktion with S for its root schmelzen, “fusion”, and Vertauschungsfunktion with T for its root tauschen, “swap”. The word schmelzen is akin to English words “melt” and “smelt”.

The “swap” is straightforward: the !!T!! combinator swaps the order of the arguments to !!x!! in !!x\,y\,z!!: $$T\,x\,y\,z = x\,z\,y$$ but does not otherwise alter the structure of the expression.

But why is !!S!! the “melting” or “fusion” combinator? It's because Schönfinkel was interested in reducing abitrary mathematical expressions to combinators. He will sometimes have an application !!(f\, x)(g\, x)!! and he wants to ‘fuse’ the two occurrences of !!x!!. He can do this by rewriting the expression as !!S\, f\, g\, x!!. Schönfinkel says:

Der praktische Nutzen der Function !!S!! besteht ersichtlich darin, daß sie es ermöglicht, mehrmals auftresnde Veränderliche — und bis zu einem gewissen Grade auch individuelle Functionen — nur einmal auftreten zu lassen.

Clearly, the practical use of the function !!S!! will be to enable us to reduce the number of occurrences of a variable — and to some extent also of a particular function — from several to a single one.

(Translation from van Heijenoort, p. 362.)

So there you have it: the !!S!! combinator is so-named not for substitution, but because S is the first letter of schmelzen, ‘to fuse’.

## References

• Schönfinkel, M. “Über die Bausteine der mathematischen Logik” (“On the building-blocks of mathematical logic”), Mathematische Annalen (1969), p. 305–316; Berlin, Göttingen, Heidelberg.

English translation in Van Heijenoort, Jean (ed.) From Frege to Gödel: a Source Book in Mathematical Logic, 1879–1931 (1967) pp. 355–366 Harvard University Press; Cambridge and London.

• Curry, H.B. “Grundlagen der kombinatorischen Logik” (“Fundamentals of combinatory logic”), American Journal of Mathematics Vol. 52, No. 3 (Jul., 1930), pp. 509-536.

• Curry, H.B. and Robert Feys Combinatory Logic (1958) p. 152 North-Holland Publishing Company, Amsterdam.