# First Order Differential Equations
Dowling16.wxmx

load(draw)$  set_draw_defaults(line_width=2, grid = [2,2], 
             background_color = light_gray,  draw_realpart=false)$ 
            fpprintprec:5$   ratprint:false$     kill(all)$

load ("Econ1.mac");

## Preface

This worksheet is one of a number of wxMaxima files available in the section
    Economic Analysis with Maxima
on my CSULB webpage.

We use Maxima to solve a few of the problems (and draw some of the plots) in Ch. 16 of the 
supplemental text: Introduction to  Mathematical Economics, 3rd ed, (Schaum's Outline 
Series), by Edward T. Dowling (1992, 2001), McGraw-Hill. This modestly priced text is a 
bargain with many worked out examples. You should compare the examples worked out 
"by hand" in this text with what we do using Maxima. Section numbers [x.y] refer to sections 
in Dowling's text.

We have slightly changed some of the symbols used in particular problems.

Ted Woollett
https://web.csulb.edu/~woollett/
Sept. 8, 2021

## References

Fundamental Methods of Mathematical Economics, 4th ed, A.C. Chiang, K. Wainwright,
McGraw-Hill, 2005, Ch. 15.

Introductory Mathematical Economics, 2nd ed., D.W. Hands, Oxford Univ. Press, 2004,
Ch. 5.

Mathematical Economics, 2nd ed., Baldani, Bradfield, & Turner, Thomson-SouthWestern,
2005, Ch. 15.

## Definitions and Concepts

Quoting Dowling:

"A "differential equation" is an equation which expresses an explicit or implicit relationship
between a function y = f(t) and one or more of its derivatives. Examples of differential
equations are
        dy/dt = 5 t + 9,       y' = 12 y,          y'' - 2 y' + 19 = 0."

"Equations involving a single independent variable (such as above) are called "ordinary
differential equations". The "solution" or "integral" of a differential equation is any equation,
without derivative or differential, that is defined over an interval and satisfies the differential
equation for all values of the independent variable in the interval."

"The "order" of a differential equation is the order of the highest derivative in the [differential]
equation. The "degree" of the differential equation is the highest power to which the derivative
of highest order is raised."

### Example 1: y'' = 7

To solve the second order differential equation y''(t) = 7 for all the functions y(t) which
satisfy this equation, simply integrate both sides of this equation to find the integrals:

   y'(t) = ∫7 dt = 7 t + c1,
   y(t)  = ∫(7 t + c1) dt = (7/2) t^2 + c1 t + c

This is called a "general solution" which indicates that when c and c1 are unspecified, 
a differential equation has an infinite number of possible solutions. If c and c1 can be 
specified,  the differential equation has a "particular solution" or "definite solution" which
alone of all possible solutions is relevant.

### Example 2: Degree vs Order

The order and degree of some differential equation examples are shown here:

dy/dt = 2 x + 6:    first order ode, first degree ode

(dy/dt)^4  - 5 t^5 = 0:   first order ode, 4th degree ode

d²y/dt²  + (dy/dt)^3  + y^2  = 0:    second order, first degree ode

(d²y/dt²)^7  + (d³y/dt³)^5 = 75 y:   third order, fifth degree ode

### Common First Order ODE's: ode2 (eqn, dvar, ivar)

### dA/dt = - k A

'diff (A,t) = - k*A;
ode2 (%, A, t);

Check solution.

A : c*exp ( -k*t);
diff (A, t) = - k*A;

### dP/dt = A - B*P

killAB()$
'diff (P, t) = A - B*P;
ode2 (%, P, t);

## First Order Linear ODE: dy/dt + v*y = z

For a first order linear ode, both dy/dt and y must be of first degree and
no product y*(dy/dt) can occur.  For such an equation
         dy/dt + v y(t) = z
where v and z may be constants or functions of t, the formula for a general
solution is

        y(t) = exp (- ∫v*dt) [ A  + ∫ z*exp (∫ v*dt)*dt ]

where A is an arbitrary constant. This solution is composed of two parts:
          yc =   A * exp (- ∫v*dt) 
is called the "complementary function"  and is the solution to the given ode
when z = 0:  dy/dt + v*y = 0, with one arbitrary constant of integration A.
The complementary function yc(t) represents the "deviation from the equilibrium".

The second part of the general solution,
           yp =   exp (- ∫v*dt) [ ∫ z*exp (∫ v*dt)*dt ]
is the  "particular integral"  and is the "intertemporal equilibrium level of y(t)".

For y(t) to be "dynalically stable", yc(t) must approach zero as t  → ∞.

The solution of a differential equation can always be checked by differentiation.

### Example 3: dy/dt + 4*y = 12

The general solution for the differential equation
   dy/dt + 4 y = 12
is calculated as follows:
Since v = 4, z = 12,
   y(t) = exp (- ∫4*dt) [ A  + ∫ 12*exp (∫ 4*dt)*dt ].

Now  ∫4*dt = 4 t + c, and when used in our general solution formula, the constant
c is always ignored and subsumed under the constant A. We then have

  y(t) = exp (- 4 t) [ A  + 12* ∫ exp ( 4 t)*dt ]

integrate(exp (4*t), t);

Hence,  y(t) = exp (- 4 t) [ A  + 3*  exp ( 4 t) ] = 3 + A exp(-4 t),  and the
complementary function --> 0 as t --> inf (infinity), y(t) --> 3, the intertemporal
equilibrium level. y(t) is dynamically stable.

Check the solution by substituting into the left hand side of the differential
equation:

y : 3 + A*exp (-4*t);
diff (y, t) + 4*y, expand;

### Example 4: dy/dt + 3*t^2*y(t) = t^2

Given the ode 
     dy/dt + 3 t^2 y(t) = t^2, 
v = 3 t^2 and z = t^2, 

integrate (3*t^2, t);

Thus the general solution is
   y(t) = exp (- t^3) [ A  + ∫ t^2*exp (t^3)*dt ]

integrate (t^2*exp (t^3), t);

So we get
   y(t) = A*exp (- t^3) + 1/3, and y --> 1/3 as t --> inf.

To check the general solution y(t) over the time interval [0, inf], we again
substitute the general solution into the left hand side (lhs) of the given
ode.

y : 1/3 + A*exp (-t^3);
diff (y, t) + 3*t^2*y, expand;

### Example 5

Suppose that y(0) = 1 in Example 4. Find the definite solution.

Since exp(0) = 1, from our general solution, we get
     y(0) = 1 = A*exp(0) + 1/3 = A + 1/3.
Solve for A = 1 - (1/3) = 2/3, then the definite solution is
     y(t) = 1/3 + (2/3)*exp(- t^3).

## First Order Nonlinear ODE, ode2 (eqn, dvar, ivar)

### Example 6:  (6 y t  + 9 y^2) dy + (3 y^2 + 8 t) dt = 0

Solve the "exact" nonlinear ode
  (6 y t  + 9 y^2) dy + (3 y^2 + 8 t) dt = 0.

This ode is an example from Dowling, Ch. 16, Sec. 3, called:
  "Exact Differential Equations and Partial Integration".
We divide both terms by dt and are presented with the 1st order ode
 of the form g(t,y)*dy/dt + f(t,y) = 0.

Instead of following along with Dowling's "by hand" methods, we use Maxima's
ode2 function, with the syntax
   ode2 (eqn, dvar, ivar)
in which eqn is a symbolic ode, dvar is the dependent variable, and
ivar is the independent variable. 

Note carefully that we use 'diff in our statement of the differential equation
(rather than diff). The single quote (') preceding diff prevents execution and
turns diff into the "noun form"   'diff.

killAB()$
eqn : (6*y*t + 9*y^2)* 'diff (y, t) + 3*y^2 + 8*t = 0;
gsoln : ode2 (eqn, y, t);

ode2 uses symbols such as %c to represent an arbitrary constant.

Take the differential of both sides using diff (gsoln), and use %c = constant
so del(%c) = 0, dt #  0, and we get back our original ode.

diff (gsoln);

ic1 (eqn, dval = d0, ival = i0) can be used to get an implicit solution in which the
arbitrary constant %c takes on a specific value.
Or you can use  ic1 (eqn, ival = i0, dval = d0)
Here we evaluate %c so that y = 1 when t = 0.

ic1(gsoln, t = 0, y = 1);

ic1(gsoln, y = 1, t = 0);

### Example 7: 5 y t y' + 5 y^2 + 8 t = 0

Solve for y(t) if the ode is
 5 y t y' + 5 y^2 + 8 t = 0.

killAB()$
eqn : 5*t*y* 'diff (y, t) + 5*y^2 + 8*t = 0;
gsoln : ode2 (eqn, y, t);

Suppose when t = 1, y = 1.

psoln : ic1 (gsoln, t = 1, y = 1);

psoln : 6*psoln;

### Example 9: y'(t) = t y^2

Solve the nonlinear ode y'(t) = t y^2.

killAB()$
eqn : 'diff (y, t) = t*y^2;
gsoln : ode2 (eqn, y, t);

Suppose when t = 0, y = 5.

psoln : ic1 (gsoln, t = 0, y = 5);

solve (psoln, y);

### Example 10:  t^2 y'(t) + y^3 = 0

Solve the nonlinear ode t^2 y'(t) + y^3 = 0.

killAB()$
eqn : t^2 * 'diff (y, t) + y^3 = 0;
gsoln : ode2 (eqn, y, t);

psoln : ic1 (gsoln, t = 1, y = 1/2);

solve (psoln, y^2);

## Economic Applications of First Order ODE's

Differential equations serve many functions in economics. They are used to determine the
conditions for dynamic stability in microeconomic models of market equilibria and to trace 
the time  path of growth under various conditions in macroeconomics. 

Given the growth rate of a function, differential equations enable the economist to find the 
function whose growth is described.

From the point elasticity, differential equations enable the economist to estimate the demand 
function  (see Example 11 and Problems 16.38 to 16.47). 

In Section 14.6 they were used to estimate capital functions from investment functions 
and total cost and total revenue functions from marginal cost and marginal revenue 
functions. 

Quoting Baldani, Bradfield, & Turner, pp 380-381:
"...our interest is in determining whether there is a steady state and, if there is, its value;
describing whether the steady state is stable, and describing the dynamics path that the
variable takes in getting to its steady-state value."

### Example 11: 

Given the demand function Qd = c + b P and the supply function Qs = g + h P, the
equilibrium price is the solution of Qd = Qs.

killAB()$
Qd : c + b*P;
Qs : g + h*P;
solns : solve (Qd = Qs, P);

Pe : at (P, solns);

We can get the numerator in the form we want by the following trick:

Pe : num(Pe)/denom(Pe);

The equilibrium price Pe = (c - g)/(h - b).

Assume the rate of change of price is proportional to the instantaneous "excess demand":
      dP/dt = m ( Qd(t) - Qs(t))
where m is a positive constant so that if demand is greater than supply, Qd > Qs, dP/dt > 0 
and price increases with time.

eqn : 'diff (P, t) = m* (Qd - Qs);

Write the ode as dP/dt = A - B*P. We will leave A and B undefined for now, but will find
and record those values in the symbols A1 and B1.

RHS : rhs(eqn), expand;

B1 : - coeff (RHS, P), factor;

A1 : RHS + B1*P, expand;

We want to write this in the form m*(c - g), and will do this the long way just to show you
how you could do it when you are dealing with a large and long expression.

Using factor doesn't work:

factor(A1);

So try using coeff.

coeff(A1,m);

A1 : m*%;

Now redefine eqn in terms of the (undefined) constants A and B.

eqn : 'diff (P, t) = A - B*P;
gsoln : ode2 (eqn, P, t), expand;

Since m > 0, B1 > 0 if h > b, and with B taking on the value of B1 > 0, P(t) approaches
A/B = A1/B1 as t gets large (and positive). Note that A1/B1 = Pe.

is (equal (Pe, A1/B1));

For normal cases where demand Qd(P) is negatively sloped ( so b < 0 ) and supply
Qs(P) is positively sloped (so h > 0), B = B1 > 0, "dynamic stability" is assured, that is
P(t) will approach the value Pe which reflects market equilibrium. The normal case has
dQd(P)/dP < 0 and dQs(P)/dP > 0.

### Prob. 16.38

Find the demand function Q(P) if point elasticity ε is -1 for all P > 0.

The dimensionless point elasticity of demand is defined by the equation
   ε = (P/Q) dQ/dP. 
Hence we are told that 
   dQ/dP = - Q/P.

killAB()$
eqn : 'diff (Q, P) = - Q/P;
gsoln : ode2 (eqn, Q, P);

The Maxima function ode2 has returned a solution which has one adjustable
constant called %c.

### Prob. 16.39;  radcan

Find the demand function Q(P) if the point elasticity of demand is
    ε = - k
where k is a constant.

killAB()$
eqn : 'diff (Q, P) = - k*Q/P;
gsoln : ode2 (eqn, Q, P);

radcan (gsoln);

The Maxima function radcan recognises that  x*log(a) = log (a^x) and
also that exp (log (A)) = A. Remember that log(A) is Maxima's expression
for the natural log of A.

exp (log(A));

Our solution is then Q(P) = constant/P^k.

### Prob. 16.40

Find the demand function Q(P) if the point elasticity of demand ε is given by
     ε = - (5 P + 2 P^2)/Q
and if Q = 500 when P = 10.

ε = (P/Q) dQ/dP = - (5 P + 2 P^2)/Q, so cancelling on both sides, we get
  dQ/dP = - 5 - 2 P.  Multiply both sides by dP and integrate both sides to get
  Q = - 5 P - P^2 + c. Evaluating both sides for P = 10, Q = 500,
 500 = - 50 - 100 + c, c = 650,
                   Q(P) = 650 - 5 P - P^2

killAB()$
gsoln : ode2 ('diff (Q,P) = - 5 - 2*P, Q, P);

solns : solve (subst ([P = 10, Q = 500], gsoln), %c);

psoln : at (gsoln, solns);

### Prob. 16.41

Derive the formula P = P(0) e^(r*t) for the total value of an initial sum of money P(0) set out
for t years at (decimal) interest rate r, when r is compounded continuously.

If r is compounded continuously, then
   dP/dt = r P

killAB()$
gsoln : ode2 ('diff (P, t) = r*P, P, t);

solns : solve (subst ([t = 0, P = P0], gsoln), %c);

psoln : at (gsoln, solns);

grind(%)$

### Prob. 16.42,  scanmap ('factor, expr)

Determine the stability condition for a two sector income determination model Y = C + I
in which we assume the rate of increase with time of the national income Y is proportional
to the "excess demand" C + I - Y, so
   dY/dt = a (C + I - Y).
We denote equilibrium income, consumption and investment are Ye, Ce, and Ie, respectively,
so Ye = Ce + Ie.  

Assume consumption and investment are both proportional to Y:
  C(t) = g Y(t),  I(t) = b Y(t), and assume 0 < a, g, b <  1.
Find Y(t) = f(Y(0), a, g, b).

killAB()$
gsoln : ode2 ('diff (Y, t) = a*(g + b - 1)*Y, Y, t);

solns : solve (subst ( [t = 0, Y = Y0], gsoln), %c);

psoln : at (gsoln, solns);

scanmap ('factor, psoln);

Dowling defines deviations dY, dC, dI in income, consumption and investment
using
  Y(t)  = Ye + dY(t), C(t) = Ce + dC(t), I(t) = Ie + dI(t).
If we substitute  into our assumed equation for dY/dt, we get
   d(dY)/dt = a (dC + dI - dY).
Dowling goes for a solution in terms of dY(t) first, before returning to get Y(t).
Since we have assumed C(t) = g Y(t), we have Ce = g Ye, 
Ce + dC = g (Ye + dY), or dC = g dY, and likewise I(t) = b Y(t), Ie = b Ye, 
Ie + dI = b ( Ye + dY), dI = b dY. Hence we get
  d(dY)/dt = a (g + b - 1) dY.

### Prob. 16.44: The Domar Model

A change in the rate of investment will affect both aggregate demand and the productive 
capability of an economy. The Domar model seeks to find the time path along which an 
economy can grow while maintaining full utilization of its productive capacity. If the marginal
propensity to save s and the marginal capital-output ratio k are constant, find the investment 
function needed for the desired growth. 

Aggregate demand Y is calculated as the sum of consumer spending, investment spending, 
government spending, and the difference between exports and imports. Whenever one of 
these factors changes and when aggregate supply remains constant, then there is a shift in 
aggregate demand.

In Dowling Ch. 14, Sec 6 is a discussion of investment and capital formation, in which
investment I(t) and capital K(t) are related by
      I(t) = dK/dt,
and we assume capital K and output Q are related by K = k Q, so  k = ∂K/∂Q = marginal
capital-output ratio, so I(t) = dK/dt = k dQ/dt, if k is constant. Thus dQ/dt = I(t)/k. We assume 
investment I = I0 + s Y, with s = ∂I/∂Y = marginal propensity to save = constant and
I0 is another constant. So
      dY/dt = (1/s) dI(t)/dt. 

We want to find I(t) needed for full-utilization growth, in which increase in demand, which
we take as dY, equals increase in productive capacity dQ, dY = dQ,  or
  dY/dt = dQ/dt,  (1/s) dI(t)/dt = I(t)/k, or dI(t)/dt = (s/k) I(t).

killAB()$
gsoln : ode2 ( 'diff (I, t) = (s/k)*I, I, t);

solns : solve ( subst ([t = 0, I = I0], gsoln), %c);

psoln : at (gsoln, solns);

Investment must grow at a constant rate r = s/k, the savings rate divided by the capital-output
ratio, for full-utilization growth of the economy. We assumed both s and k were constant
numbers in this derivation.

### Prob. 16.44: The Solow Model

The Solow model examines equilibrium growth paths with full employment of both capital 
and  labor. Based on the assumptions that  

1. Output is a linearly homogeneous function of capital and labor exhibiting constant  
returns to scale: Y = f(K, L), with (a > 0):   f(a K, a L) = a f(K, L),  which implies
          f(K, L) = a f(K/a, L/a),

2. A constant proportion s of output is saved and invested: 
             I(t) = dK/dt = s Y,  s = constant,

3. The supply of labor is growing at a constant rate r: 
                 L = L0 exp (r t),

derive the differential equation in terms of the single variable K/L, which serves as the basis 
of  the model.  

I = dK/dt = s Y = s f(K, L).
Assume z(t) = K/L, then K = z L. Then 
  dK/dt = d (z L)/dt = z dL/dt + L dz/dt  = s f(K, L) = s L f(K/L, 1),
  or z r L + L dz/dt = s L f(z,1).
Divide both sides by L, and find an ordinary differential equation for z(t) = K/L in
terms of parameters r and s:
            dz/dt = s f(z,1) - r z.

### Prob. 16.46:  ratsubst, limit

Assume that the demand for money is for transaction purposes only. Thus, 
   Md = k P(t) Q
where P(t) is the price level, k is a constant, Q is output, assumed constant.
Assume the money supply Ms is exogenously determined by the Fed and money supply
equals money demand when price level equilibrium is reached: Ms = k Pe Q, so
the equilibrium price level is Pe = Ms/(k*Q).

If inflation or the rate of change of prices is proportional to excess demand for goods, and
from Walras' law, an excess demand for goods is the same thing as an excess supply of
money, then we can assume
     dP/dt = b (Ms - Md), or dP/dt = b Ms - b k Q P(t)
and we can solve for the path of prices in time, beginning with P(0) = P0.

killAB()$
gsoln : ode2 ('diff (P,t) = b*Ms - b*k*Q*P, P, t), expand;

gsoln1 : ratsubst (Pe, Ms/(Q*k), gsoln), expand;

solns : solve (subst ([t=0, P = P0], gsoln1), %c);

psoln : at (gsoln1, solns);

Since each of Q, b, and k is positive, the price level approaches Pe as t --> infinity.

assume (A > 0)$
limit (exp (- A*t), t, inf);

### Prob. 16.47

If the expectation of inflation is a positive function of the present rate of inflation
   [dP/dt]_E = h dP/dt, h > 0,
and the expectation of inflation reduces people's desire to hold money, so that
   Md = k P(t) Q - g [dP/dt]_E = k P Q - g h dP/dt, g > 0,
check the price stability conditions assuming that the rate of inflation is proportional to the
excess supply of money (as in Prob. 16.46):
  dP/dt = b (Ms - Md), b > 0, Ms > 0 and fixed, then
dP/dt = b Ms - b ( k P Q - g h dP/dt ).

killAB()$
gsoln : ode2 ('diff (P, t) = b*Ms - b*k*Q*P + b*g*h*'diff (P,t), P, t);

gsoln : ratsubst (Pe, Ms/(Q*k), gsoln);

solns : solve ( subst ( [t = 0, P = P0], gsoln), %c);

psoln : at (gsoln, solns);

For a stable approach to the equilibium price Pe = Ms/(Q*k), we need the product b*g*h < 1.
Even if h > 1, so that people expect inflation to accelerate, the economy need not be
unstable, as long as b and g are sufficiently small.

## Direction Fields, plotdf (dydx), plotdf (dvdu, [u,v])

### Direction Fields of a Single First-Order ODE, plotdf (dydx)

A "direction field" associated with the ode
      'diff(y,x) = f(x,y)
is a plot of the short line segments which show the local slope line dy/dx at a large set
of points (x,y).

To plot the direction field of a single ODE written in a form in which y is the dependent
variable and x is the independent variable, such that f(x,y) gives the local slope at (x,y):
      dy/dx = f(x, y)
the simplest invocation (without other possible options) is:
   plotdf ( f(x,y) );

If you instead let v(u) be the dependent variable and u the independent variable, you must
include a list [u, v] as the second argument, without other options, just:
   plotdf (f(u,v), [u,v] );
which shows the short line segments which give the local slope dv/du at (u,v), and are
governed by the first order ode 'diff (v,u) = f(u,v).

When called, plotdf opens a new window that contains the direction field plot. The length 
of the draw vectors is proportional to the local absolute value of the first derivative.

A direction field for the differential equation
    dy/dx  = exp (- x) + y
is created by the following Maxima code.

killAB()$
plotdf(exp(-x)+y);

This Maxima command will work as long as the Xmaxima interface files have been selected
when installing Maxima (as well as wxMaxima files). A separate plotdf window opens which
shows the coordinates (x,y) (in the lower right hand corner) of your cursor as you move the 
cursor across the direction field plot.

There is no present way to create a plotdf window image directly as part of a normal cell in 
wxMaxima, but one can import an image into a cell (as long as the original graphics file is 
in your work folder). The author's work folder is c:\work5. Inside wxMaxima, use the Cell
menu from the menu bar at the top, then choose "import image", and then browse to the
image file, select the file,  and then choose open.

maxima_userdir;

In that folder is the file pdf01.jpg.

file_search ("pdf01.jpg");

If you choose Cell, Insert image... from the wxMaxima menu along the top, you can browse
to the graphics file c:\work5\pdf01.jpg, and insert that image here in the worksheet, and it
will continue to show up as long as that file remains present. This jpg image was created
in Windows by pressing the PrintScreen button once, and opening the freely available
graphics editor Irfanview, and choosing Edit, Paste, and then choosing File, Save As, to
save to a temporary pics folder.

C:\work5\pdf01.jpg

I double clicked on that file in the file listing of the temporary pics folder to have the default 
Windows Photos app open and display the graphics, then selected Edit & Create, Edit, and
grabbed the corners to crop down to just the plotdf window, and when the image looked
suitably cropped, chose Save a Copy in the lower right corner. This allowed the cropped
image to be saved finally to c:\work5\pdf02.jpg.

file_search ("pdf02.jpg");

C:\work5\pdf02.jpg

Since dy/dx = exp(- x) + y, when (x,y) = (0, 0), we get dy/dx = 1, and you see a small arrow
pointing at a 45 degree angle with the base of the arrow at (x=0, y=0).

At x = 5, y = -5, dy/dx = exp(- 5) - 5 ~ -5.

exp(-5) -5, numer;

and at (5, -5) we see an arrow five times in length as the (0,0) arrow, and pointing almost
straight down. At that location, y(x) is a rapidly decreasing function of x.

In contrast, at (x = 5, y = 5), dy/dx = exp( - 5) + 5 ~ 5,

exp(-5) + 5, numer;

and we see an arrow pointing almost straight up; y(x) is rapidly increasing with x at (5, 5).

### plotdf Window Menu Bar

The menu bar of the plotdf window has the following seven icons:

(1)   An X. 
          Used to close the plot window.

(2)  A wrench and a screwdriver. 
       This opens the configuration menu with several fields that show the ODE(s) 
           in use and various other settings. If a pair of coordinates are entered in 
           the field "Trajectory at" and the keyboard Enter key is pressed, a new integral 
           curve will be shown, in addition to the ones already shown. You might want to
           increase the linewidth of the line which shows a trajectory, and perhaps the color,
           and then use icon (3) to replot using your choices.

(3)  Two arrows forming a circle. 
         Replots the direction field with the new settings defined in the configuration 
         menu and replots only the last integral curve that was previously plotted.

(4)  Hard disk drive with an arrow. 
         Used to save a copy of the plot, in Postscript format, in the file specified in a 
           field of the box that appears when that icon is clicked.

(5)  Magnifying glass with a plus sign. 
          Zooms into the plot.

(6)  Magnifying glass with a minus sign. 
       Zooms out of the plot. The plot can be displaced by holding down the right 
           mouse button while the mouse is moved.

(7)  Icon of a plot. 
        Opens another window with a plot of the two variables in terms of time, 
            for the last integral curve that was plotted.

By zooming in and displacing the plot using the right mouse button held down, a larger
plot called pdf03.jpg was created after some experimentation.

file_search ("pdf03.jpg");

C:\work5\pdf03.jpg

Working with the same ode, let the independent variable be t and the dependent variable
be u, so our ode looks like du/dt = exp (-t) + u(t). We have to tell plotdf what the variables
are called in this case, with the list [t, u] as the second argument. We are also going to
include an optional trajectory, starting at t = 2 and u = - 0.1. We are going to tell plotdf
the extent of the t axis to show and the extent of the u axis, using [t, t1, t2], [u, u1, u2].
We are also telling plotdf to proceed in the direction of increasing t for the trajectory,
by including [direction, forward].

plotdf(exp(-t) + u, [t, u], [trajectory_at, 2, -0.1],
[direction, forward], [t, 0, 7], [u, -6, 1] );

C:\work5\pdf04.jpg

## Phase Diagram Analysis of Stability

Many nonlinear differential equations cannot be solved explicitly as functions of time. 
Phase  diagrams, however, offer qualitative information about the stability of equations that is 
helpful in determining whether the solution of the differential equation will converge to an 
intertemporal (steady-state) equilibrium or not (as time passes).

A phase diagram of the differential equation dx(t)/dt = f(x) depicts the derivative dx/dt which 
we now express as vx, for simplicity of notation, as a function of x. 

vx (= f(x) = dx/dt) is plotted on the vertical axis versus the value of x (which naturally 
corresponds with a specific time t) on the horizontal axis. The number vx can be thought of
as the x-component of the velocity vector of a particle executing one dimensional motion 
x(t) along an x-axis.

When vx > 0, the velocity vector points to the right, and the "particle" is moving to the right. 
When vx < 0, the velocity vector points to the left, and the particle is moving to the left.

A possible steady-state solution is easily identified on this phase diagram as any point at 
which the graph of vx versus x crosses the horizontal x-axis, because at that point vx = 0, 
the particle is momentarily at rest. Under the right conditions, once the particle reaches
that value of x for which vx = 0, it may stay there.  For some equations there may be more 
than one intersection and hence  more than one value of x for which vx = 0.  

Diagrammatically, the stability of each of the values of x for which vx = 0 is revealed by the
directions of the velocity vector just to the left and just to the right of that point x.

Dowling calls these velocity vectors "arrows of motion", and we can continue to use his
language when we replace x(t) by P(t), the market price as a function of time, for example.

The arrow of motion will point to the RIGHT (indicating x is increasing with time) any time 
the graph of vx = dx/dt = f(x) is ABOVE the horizontal x-axis and hence vx > 0.  The x-axis, 
by definition, is the locus of points for which vx = 0.

The arrow of motion will point to the LEFT (indicating x is decreasing with time) any time 
the graph of vx = dx/dt = f(x) is BELOW the horizontal x-axis and hence vx < 0.
 
If the arrows of motion (on both the left and right sides) point TOWARDS a location (x value)
where vx = 0, that solution is stable (the implied "forces" are pushing the particle back to
that value of x).

If the arrows of motion (on both the left and right sides) point AWAY FROM a location 
(x value) where vx = 0, that solution is unstable (the implied "forces" are pushing the particle
away from that value of x).
 
Mathematically, the SLOPE of the phase diagram path as it passes through a possible 
steady-state equilibrium point tells us if the candidate point is stable or not. 

When evaluated at a possible steady-state equilibrium point x*,

  (1)  IF  d(vx)/dx < 0, the critical point x* is a point of stable equilibrium.

  (2)  IF  d(vx)/dx > 0, the critical point x* is a point of unstable equilibrium.

Consider case (1). At x* we have vx = dxdt = 0. Since vx is decreasing as x increases
(in the neighborhood of x*), for x > x* (to the right of x*) we must have vx < 0 and the 
velocity vector (arrow of motion) points to the left, toward x*. And for x < x* (to the left of x*) 
we must have vx > 0 and the velocity vector points to the right, toward x*. On both sides 
the "forces" are pushing the "particle" back to x*.

Likewise in case (2), vx = 0 at x = x*, and since vx is increasing as x increases in this
neighborhood, for x > x* (to the right of x*) we must have vx > 0 and the velocity vector
points away from x*. To the left of x*, x < x*, we must have vx < 0, and the velocity 
points to the left, away from x*. On both sides the forces are pushing the particle away
from x  = x*.

Here is an example in which vx = dx/dt is given as a function of x.

Let vx = x^2 - 9. 

The points where vx = 0 (the particle is momentarily at rest) are potential
points of stable equilibrium for one dimensional motion. Points where vx > 0 are points
where x is growing with time and getting larger. Points where vx < 0 are points where
dx/dt < 0 and x is shrinking in value and getting smaller with time.

For vx = x^2 - 9, vx = 0 at x = +/- 3, so we have two candidate points x* where vx = 0.

Since the crucial stability test depends on the sign of d(vx)/dx at x* = (-3,3), we can take
the first derivative to get slope = d(vx)/dx = 2 x, which will be negative (stable) for x < 0
and positive (unstable) for x > 0.

We can also calculate the x-component ax of the "particle" acceleration:
    ax = d(vx)/dt = (dx/dt) * [d(vx)/dx] = vx*slope, using the "chain rule" of calculus.

killAB()$
[vx : x^2 - 9, slope : 2*x, ax : vx*slope];
wxdraw2d (xlabel = "x", ylabel = "vx, slope, ax", key_pos = top_center, color = red,
   yrange = [-50, 50], key = "vx = dx/dt", explicit (vx, x, -5, 5), 
    color = blue, key = "ax = dvx/dt", explicit (ax, x, -5,5),
    color = purple, key = "slope = d(vx)/dx", explicit (slope, x, -5, 5),
    key = "", color = black, line_width = 1, explicit (0,x,-5,5),
    parametric (-3,yy,yy,-50,0), parametric (0,yy,yy,-50,0), parametric (3,yy,yy,-50,0));

at (slope, x = -3);

So x* = -3 is a point of stable equilibrium. Look at the values of vx and ax to the LEFT
of x* = -3, say at x = -4.

at ([vx, ax], x = -4);

So at x = -4 the velocity vector (x-component) is vx =  +7 > 0, the particle is moving 
toward x*, and the acceleration vector (x-component) is ax = -56, so the acceleration
vector points to the left (slowing down the speed of the particle) and the dot product
of the velocity vector with the acceleration vector is vx*ax < 0, and vectors point in opposite
directions.

Look at the values of vx and ax to the right of the stable point x* = -3, say x = -2.

at ([vx, ax], x = -2);

So at x = -2, the velocity vector points to the left, towards x = 3, and the acceleration vector
points to the right, again slowing the speed of the particle as it approaches x = -3. Again
the dot product of v_vec and a_vec is a negative number, they point in opposite directions.

Let's add some "arrows of motion" (particle velocity vectors) to our plot.

wxdraw2d (xlabel = "x", ylabel = "vx, slope", key_pos = top_center, color = red,
   yrange = [-50, 50], key = "vx = dx/dt", explicit (vx, x, -5, 5), 
    color = purple, key = "slope = d(vx)/dx", explicit (slope, x, -5, 5),
    key = "", color = black, line_width = 1, explicit (0,x,-5,5),
    parametric (-3,yy,yy,-50,0), parametric (0,yy,yy,-50,0), parametric (3,yy,yy,-50,0),
    color = blue, line_width = 5, head_length = 0.2, vector ([- 4.5, 0], [1, 0]),
     vector ( [-1.5, 0], [-1, 0]), vector ( [2.5, 0], [-1, 0]), vector ( [3.5, 0], [1, 0]))$

### Dowling Example 12

Consider the ode dy/dt = 8 y - 2 y^2. 
Let vy = dy/dt = 8 y - 2 y^2.

killAB()$
[vy : 8*y - 2*y^2, slope : diff (vy, y), ay : expand (vy*slope)];

solns : solve (vy, y);

[y1, y2] : map ('rhs, solns);

at (slope, y = y1);

at (slope,y = y2);

So y* = 4 is a stable equilibrium point and y* = 0 is not.

wxdraw2d (xlabel = "y", ylabel = "vy, slope/10, ay/10", key_pos = top_center, color = red,
   yrange = [-4, 12], key = "vy = dy/dt", explicit (vy, y, -1, 5), 
    color = blue, key = "ay/10", explicit (ay/10, y, -1, 5),
    color = purple, key = "slope/10 ", explicit (slope/10, y, -1, 5),
    key = "", color = black, line_width = 1, explicit (0,x,-1,5),
    parametric (0,yy,yy,-4,0), parametric (4,yy,yy,-4,0),
    color = blue, line_width = 5, head_length = 0.1, vector ([- 0.2, 0], [-0.5, 0]),
     vector ( [ 0.2, 0], [ 0.5, 0]), vector ( [3.2, 0], [ 0.5, 0]), vector ( [ 4.8, 0], [ -0.5, 0]))$

The arrows of motion complete the graph. As explained above, where the graph lies above 
the horizontal y-axis, vy > 0 and the arrows must point to the right; where the graph lies 
below the horizontal y-axis, vy < 0  and the arrows must point to the left.
  
The stability of the points y* for which vy = 0 can now be read from the graph. Since 
the arrows of motion point away from y* = 0, that y is an unstable equilibrium. 
With the arrows of motion pointing toward y* = 4, that value of y is a stable equilibrium point
for the time evolution determined by the given first order autonomous ode.

The slope of the phase diagram at the steady-state solutions can also be used to test 
stability independently of the arrows of motion. Since the slope of the phase diagram is 
positive at y1 = 0, we can conclude y1 is an unstable equilibrium. Since the slope of the 
phase diagram is negative at y2 = 4, we  know y2 must be stable. 