# Second Order Difference Equations
Dowling18B.wxmx

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

load ("Econ2.mac");

## Preface

Dowling18B.wxmx is one of a number of wxMaxima files available in the section
    Economic Analysis with Maxima
on my CSULB webpage. Some of the calculations done in this worksheet need to have
the software file Econ2.mac loaded first. Use load ("Econ2.mac");

We use Maxima to solve a few of the problems (and draw some of the plots) in Sec. 18.2
[Second Order Difference Equations] and in Sec. 18.8 [Stability Conditions] of Ch. 18 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 supplemented Dowling's problem 18.37 (Samuelson's Multiplier-Accelerator Model)
with the more extensive treatment of this model by Chiang and Wainwright (their Ch. 18,
Sec. 2).

We have slightly changed some of the symbols used in particular problems. An approximate
pdf translation (using Microsoft print to pdf) is available as Dowling18Bfit.pdf. This pdf file
can be searched using Ctrl-F.

Ted Woollett
https://home.csulb.edu/~woollett/
Dec. 9, 2021

## References

Fundamental Methods of Mathematical Economics, Alpha C. Chiang and Kevin Wainwright,
4th ed., 2005, McGraw-Hill (Ch. 18)

https://en.wikipedia.org/wiki/Linear_difference_equation

https://en.wikipedia.org/wiki/Multiplier-accelerator_model

"How to Solve Difference Equations? - A Complete Video Tutorial", 35 min,
 Dr. Mostapha Kalami Heris,
https://www.youtube.com/watch?v=gI2JIwIVk8E

Math. Methods for Economics, 22 min,
https://www.youtube.com/watch?v=u1ec0t1doAM

## Concepts and Definitions

We first review notational conventions which were presented in Dowling17.wxmx which dealt 
solely with first order difference equations.

A difference equation expresses a relationship between a dependent variable and a lagged 
independent variable (or variables) which changes at discrete intervals of time, for example, 
    I[t] = F( Y[t - 1] )  
with F(x) being some arbitrary function, where I and Y are each measured at the end of 
each year. Investment during year 4 is some function of national income during year 3, for
example: I[4] = F (Y[3]), etc.

The order of a difference equation is determined by the greatest number of periods 
lagged. A first-order difference equation contains only a time lag of one period, as above.
 
A second-order difference equation involves a time lag of two periods, in addition to
 (perhaps) also a time lag of one period, etc. 

The change in y as t changes from t to t +1  is called the "first difference of y". 
The FIRST DIFFERENCE of y is written :
  Δy/Δt  =  Δy[t] = y[t + 1] - y[t]  = change in y per period over the time (t, t+1).

In Maxima, y[t] is a (hash) array and wxMaxima will display the symbol using subscript
notation ( if display2d is set to true, which is the default).

y[t];

display2d;

### Examples of Difference Equations of Various Orders

Each of the following is a difference equation of the order indicated: 

Order 2:

I[t] = a* ( Y[t-1] - Y[t-2]);

Order 1:

Q[t] = a + b*P[t-1];

Order 3:

y[t+3] - 9*y[t+2] + 2*y[t+1] + 6*y[t] = 8;

Order 1:

Δy[t] = 5*y[t];

The last equation is the same as:

y[t+1] - y[t] = 5*y[t];

which is the same as:

y[t+1] = 6*y[t];

## Linear Second Order Difference Equations

### Total Solution: y[t] + b1*y[t-1] + b2*y[t-2] = a

The total solution y[t] of a second order linear difference equation is the sum
of the particular solution yp[t] and the complimentary solution yc[t]:
  y[t] = yp[t] + yc[t],
where yc is the general solution of the homogeneous equation

### ypart(b1,b2,a,t), ycompl(b1,b2,t), ysoln (b1,b2,a,t)

Given the 2nd order linear difference equation
    y[t] + b1*y[t-1] + b2*y[t-2] = a, 
with constants b1,b2, a, the particular solution yp[t] is returned by the Maxima function
ypart (b1,b2,a,t) defined in Econ2.mac.

(ypart, ycompl, and ysoln are all defined in Econ2.mac.)

fundef (ypart);

ypart (b1,b2,a,t);

ypart (-2, 1, 8, t);

yc[t] is given by the Maxima function ycompl (b1, b2, t).

ysoln (b1,b2,a,t) returns yp + yc in terms of %A1 and %A2, an "indefinite" solution
for y[t].

ycompl(-2,1,t);

ysoln (-2,1,8,t);

## Examples 5, 6 & 7

### Eqn. 1: Real Roots,  b1^2  >  4*b2

Find the particular solution and the complementary solution of the second order
linear differential equation
  y[t] - 10*y[t-1] + 16*y[t-2] = 14, in which b1 = -10, b2 = 16, a = 14.

ypart (-10, 16, 14, t);

ycompl (-10, 16, t);

ysoln (-10, 16, 14, t);

### Eqn. 2: Real Roots,  b1^2  >  4*b2

Find the particular solution and the complementary solution of the second order
linear differential equation
  y[t] - 6*y[t-1] + 5*y[t-2] = 12, in which b1 = -6, b2 = 5, a = 12.

ypart (-6, 5, 12, t);

ycompl (-6, 5, t);

ysoln (-6, 5, 12, t);

### Eqn. 3, Double Real Root, b1^2 = 4*b2

Find the particular solution and the complementary solution of the second order
linear differential equation
  y[t] - 2*y[t-1] + y[t-2] = 8, in which b1 = -2, b2 = 1, a = 8.

ypart (-2, 1, 8, t);

ycompl (-2, 1, t);

ysoln (-2, 1, 8, t);

## Example 8, Satisfying y[0] & y[1] Conditions

We return to Eqn. 1 above and require y[0] = 10 and y[1] = 36 and look for a definite
solution respecting that data.  We first define a Maxima function ys(t). We need to use
two single quotes  (') in front of yindef here.

yindef : ysoln (-10, 16, 14, t);

ys(t) := ''yindef;

solns : solve ([ ys(0) = 10, ys(1) = 36], [%A1, %A2]);

ydef : at (yindef, solns[1]);

We can check ydef by evaluating at t=0,1,2. It is easiest to define a Maxima function yd(t).

yd(t) := ''ydef;

map ('yd, [0,1,2,3]);

Alternatively, you can use the Maxima at( ...) function to check the solution.

at (ydef, t = 0);

at (ydef, t = 1);

## Prob. 18.16: Repeated Real Roots

Find the solution y[t] for the difference equation
y[t] - 10 y[t-1] + 25 y[t-2] = 8, for which b1 = -10, b2 = 25, a = 8.

ysoln (-10, 25, 8, t);

grind(%)$

## Definite Solutions, ytlist(b1,b2,a,y0,y1,tmax)

### Stability of Solutions

If the roots are real, the characteristic root with the largest magnitude (absolute value) is 
called the dominant root because it dominates the time path. For convergence, the 
absolute value of the dominant root must be less than 1.

The "characteristic equation" corresponding to the difference equation
 y[t] + b1*y[t-1] + b2*y[t-1] = a
is
   r^2 + b1*r + b2 = 0.
The solutions of this quadratic equation are the "characteristic roots", and can either be 
    two real distinct roots r1, r2, if b1^2 > 4*b2,
    a double real root r,  if b1^2 = 4*b2,
    or a pair of complex roots,  if b1^2 < 4*b2, which then also means b2 > 0.

For the case of complex roots, stability is determined by the magnitude of k = sqrt (b2).

If |k| < 1 the solution converges, if |k| > 1 the solution diverges.

### Prob. 18.18: Distinct Real Roots: b1^2 > 4*b2

Find the definite solution of
y[t] + 7*y[t-1]  + 6*y[t-2] = 42, with y[0] = 16, y[1] = -35.
Comment on stability.

For this problem, b1 = 7, b2 = 6, a = 42, y0 = 16, y1 = -35.

We first look at this problem the "long way", as a pedagogical exercise in Maxima.

We then use Maxima our do loop code ytlist, which is much faster.

Look at the characteristic roots with b1 = 7 and b2 = 6:

solve (r^2 + 7*r + 6);

yindef : ysoln ( 7, 6, 42, t);

Define a Maxima function ys(t) by placing two single quotes (') in front of yindef.

ys(t) := ''yindef;

Find the definite solution by using the values of y[0] and y[1].

solns : solve ([ ys(0) = 16, ys(1) = -35], [%A1, %A2]);

ydef : at (yindef, solns[1]);

The characteristic roots are -1 and -6. The dominant root is -6. Since |-6| = 6 > 1, the 
time path is divergent.

Again, use two single quote marks in front of ydef to define the Maxima function yd(t).

yd(t) := ''ydef;

Make a list of times and a list of values of y[t].

tmax : 5$
tL : makelist(t,t,0.0,tmax);

ytL : map ('yd, tL), numer;

We see oscillation between positive and negative values, each point further from the
intertemporal equilibrium value ye = 3; explosive instability.

wxdraw2d (xlabel = "t", ylabel = "y[t]", xrange = [0,tmax],yrange = [-1.1e3,1.1e3],
    title = "b1 = 7, b2 = 6, a = 42, y[0] = 16, y[1] = -35",
    points (tL, ytL), line_width = 1, color = red, key = "ye = 3", explicit(3, t, 0, tmax))$

Now for the faster path. A Maxima function 
      ytlist (b1, b2, a, y0, y1, tmax)
is defined in Econ2.mac. ytlist produces a list of values of y[t] from the given
second order difference equation.

tmax : 5$
tL : makelist (t, t, 0, tmax);
ytL : ytlist (7, 6, 42, 16, -35, tmax);

wxdraw2d (xlabel = "t", ylabel = "y[t]", xrange = [0,tmax],yrange = [-1.1e3,1.1e3],
    title = "b1 = 7, b2 = 6, a = 42, y[0] = 16, y[1] = -35",
    points (tL, ytL), line_width = 1, color = red, key = "ye = 3", explicit(ye, t, 0, tmax))$

A look at longer times (without a plot):

float ( ytlist (7, 6, 42, 16, -35, 30) );

### Prob. 18.20: Repeated Real Roots: b1^2 = 4*b2

Find the definite solution of
y[t] - 10*y[t-1]  + 25*y[t-2] = 8, with y[0] = 1, y[1] = 5.
Comment on stability.

b1 = -10, b2 = 25, a = 8, y0 = 1, y1 = 5.

We take the fast path using the do loop code ytlist(b1,b2,a,y0,y1,tmax)
defined in Econ2.mac.

tmax : 5$
tL : makelist (t, t, 0, tmax);
ytL : ytlist (-10, 25, 8, 1, 5, tmax);

wxdraw2d (xlabel = "t", ylabel = "y[t]", xrange = [0,tmax],yrange = [-10,1.5e3],
    title = "b1 = -10, b2 = 25, a = 8, y[0] = 1, y[1] = 5", key_pos = top_left,
    points (tL, ytL), line_width = 1, color = red, key = "ye = 0.5", explicit(ye, t, 0, tmax))$

A look at longer times without a plot:

float (ytlist (-10, 25, 8, 1, 5, 15) );

### Prob. 18.34: Complex Roots:  b1^2 < 4*b2

Find the definite solution of
y[t] + 4*y[t-2] = 15, with y[0] = 12, y[1] = 11.
Comment on stability.

For this problem, b1 = 0, b2 = 4, a = 15. The characteristic equation is r^2 + 4 = 0.
So r = +/- sqrt(- 4) = +/- (- 4)^(1/2) = +/- %i*(4)^(1/2) = +/- 2*%i.

Stability for this complex root model depends on the magnitude of 
  k = sqrt(b2) = sqrt(4) = 2 > 1 so unstable.

tmax : 7$
tL : makelist (t, t, 0, tmax);
ytL : ytlist ( 0, 4, 15, 12, 11, tmax);

wxdraw2d (xlabel = "t", ylabel = "y[t]", xrange = [0,tmax],yrange = [-700, 600],
    title = "b1 = 0, b2 = 4, a = 15, y[0] = 12, y[1] = 11", 
    points (tL, ytL), line_width = 1, color = red, key = "ye = 3", explicit(ye, t, 0, tmax))$

Let's look at longer times.

float (ytlist ( 0, 4, 15, 12, 11, 15));

We see continued oscillation between positive and negative values, with the steps increasing
in size. We see oscillating explosive instability.

## Maxima Solutions Using solve_rec (eqn, var, [init])

The contributed Maxima package solve_rec.mac can be used to solve recurrence
relation equations.

load (solve_rec);

### General Case

kill (b1, b2, a)$
gsoln : solve_rec (y[t] + b1*y[t-1] + b2*y[t-2] = a, y[t]);

### Numerical Examples

ysoln (-10, 16, 14, t);

gs1 : rhs (at (gsoln, [b1 = -10,b2 = 16,a = 14]));

With this version of Maxima, we have to combine the numerator and denominator
A^t/B^t = (A/B)^t by hand in the first two terms to get factors 8^t and 2^t.

Here is a second numerical example:

ysoln (-2, 1, 8, t);

If you try: gs2 : rhs (at (gsoln, [b1 = -2, b2 = 1, a = 8]));
you will get a "divide by zero"  error message: undefined: 0 to a negative exponent.
You need to start again with solve_rec (eqn,var) from scratch.

solve_rec (y[t] + -2*y[t-1] + y[t-2] = 8, y[t]);

Redoing Example 8 above, with initial conditions,

solve_rec (y[t] -10*y[t-1] + 16*y[t-2] = 14, y[t], y[0] = 10, y[1] = 36);

grind(%)$

Which agrees with our expression found above for ydef  in Example 8.

## Multiplier-Accelerator model: a = acceleration,
b = propensity to consume (MPC), IaG0 = Ia+G0

The multiplier–accelerator model (also known as Keynes-Hansen–Samuelson model) is a 
macroeconomic model which analyzes the business cycle. This model was developed by
Paul Samuelson, who credited Alvin Hansen for the inspiration. This model is based on 
the Keynesian multiplier, which is a consequence of assuming that consumption intentions 
depend on the level of economic activity, and the accelerator theory of investment, which 
assumes that investment intentions depend on the pace of growth in economic activity.

One version of the multiplier–accelerator model follows:

First, the market-clearing level of economic activity is defined as that at which production 
exactly matches the total of government spending intentions, households' consumption 
intentions and firms' investing intentions (see below):
  Y[t] = G[t] + C[t] + I[t],
in which investment is the sum of autonomous investment Ia (a constant) and induced
investment Id[t];

then an equation to express the idea that households' consumption intentions depend upon 
some measure of economic activity, possibly with a lag (b is the marginal propensity to
consume (MPC),  (0 < b < 1):

  C[t] = b Y[t-1] = consumption expenditure in year t;

then an equation that makes firms' investment intentions react to the pace of change of 
economic activity ('a' stands for the accelerator, the "acceleration coefficient", with a > 0):

  Id[t] = a (C[t] - C[t-1] ) = a (b Y[t-1] - b Y[t-2] ) = a*b (Y[t-1] - Y[t-2]) = induced private firm 
investment in plant and capital in year t;

and finally a statement that government spending intentions are not influenced by any of 
the other variables in the model. We will take G[t] = G0, and combine the constant
government spending with the autonomous investment (by firms) Ia, calling the sum IaG0,
    IaG0 = Ia + G0.
We could take IaG0 = 1, and then all other numbers are multiplies of IaG0.

We can rearrange these equations and rewrite them as a second-order linear difference 
equation in the standard form we have been using.

      Y[t] - b*(1+a)*Y[t-1] + a*b*Y[t-2] = IaG0,

in which our constants b1 = -b(1+a), b2 = a*b, and "a" = IaG0.

Based on the signs of a and b: 0 < b < 1, a > 0, (1+a) > 1, b*(1+a) > b > 0, 
b1 = - b*(1+ a) < 0, and  b2 = a*b > 0.
 
We have two distinct real roots of the characteristic equation r^2 - b*(1+a)*r + a*b = 0
 if b1^2 > 4*b2, which here translates to  b > 4*a/(1+a)^2.

We have a double real root if b1^2 = 4*b2,  which translates to  b = 4*a/(1+a)^2.

We have complex roots if b1^2 < 4*b2 (possible here since b2 > 0), so if b < 4*a/(1+a)^2.
In this complex root case, k = b2^(1/2) = sqrt (a*b) > 0, and 
   we need |k| < 1 for convergence.

Let yp be the particular solution of the complete equation.

kill(a,b)$
yp : ypart (- b*(1+a), a*b, IaG0, t);

The denominator simplifies if you expand it out.

yp : num(yp)/expand(denom(yp));

The expression 1/(1-b) is the multiplier that would multiply IaG0 the absence of
induced investment. The quantity IaG0/(1- b) is the "intertemporal equilibrium income" Ye.

### yMA (a, b, IaG0, t) Multiplier-Accel. Indefinite Solution

We can piggy-back on our previous work and define a solution for Y[t] using a new Maxima
function yMA (a, b, IaG0, t), since we know how to relate b1 and b2 to the model parameters
a and b. Recall a is the acceleration and b is the propensity to consume and 
 b1 = - b*(1+ a) < 0,   b2 = a*b > 0.

yMA (aa, bb,iag0, tt) := ysoln (-bb*(1 + aa), aa*bb, iag0, tt)$

### yMA(a, b, IaG0, t)  Method Example

As an example, we assume a = 0.5, b = 0.9, IaG0 = 1, Y[0] = 1, Y[1] = 2.
The first step is to use yMA to calculate the indefinite solution; that solution includes two
arbitrary constants %A1 and %A2.

yindef : yMA (0.5, 0.9, 1, t);

To define a Maxima function ys(t), use two single quotes (') in front of yindef:

ys(t) := ''yindef;

solns : solve ([ys(0) = 1, ys(1) = 2]), numer;

Now define the expression ydef (the definite solution, taking into account the two known
values of Y[t]).

ydef : at (yindef, solns[1]);

To easily generate a list of values of Y[t], define the Maxima function yd(t), again using two
single quotes (') in front of the expression ydef:

yd(t) := ''ydef;

tmax : 10$
tL : makelist (t, t, 0, tmax);

ytL : map ('yd, tL);

An alternative route to a list of values of Y[t] is the following:

makelist (subst (tL[j], t, ydef), j, 1, length (tL));

We can then make a plot of the values of Y[t] over the chosen interval of times. We make
use of the draw2d function points (xL,yL).

wxdraw2d (xlabel = "t", ylabel = "Y[t]", xrange = [0,tmax],yrange = [0, 12],
    title = " a = 0.5, b = 0.9, IaG0 = 1, Y[0] = 1, Y[1] = 2",
    key_pos = bottom_right, points (tL, ytL), line_width = 1, color = red, 
    key = "Ye = 10", explicit(10, t, 0, tmax))$

We can map yd onto any list of times to see later behavior: for example from t = 10 to
t = 20:

map ('yd, makelist (t, t, 10, 20));

### MAlist(a,b,IaG0,Y0,Y1,tmax), makeLists(L) Method Example

To examine details (consumption C[t] and induced investment Id[t]), we use the function
MAlist (a, b, IaG0, Y[0], Y[1], tmax), defined in Econ2.mac, which uses Maxima do loop
code. This Maxima function creates a list of lists which can then be used to make a table 
(using matrix) and also a plot of C[t] and Id[t] as well as Y[t]. Indeed, this do loop code
is a faster route to see what is happening to just Y[t] itself. The function MAlist also 
produces a global value for Ye = IaG0/(1-b).

As an example, we choose a = 0.5 = 1/2, b = 0.9, IaG0 = 1, then Ye = 1/(1-b) = 10. 
We assume Y[0] = 1, Y[1] = 2, and then use:
         Y[t] = C[t] + Id[t] + IaG0 = b*Y[t-1] + a*b*(Y[t-1] - Y[t-2]) + IaG0, 
to advance the solution to the next integral time.

Since we are going to do two things with the output of MAlist, we assign the output (a
list of lists) to the name maL. The sublists start with a list of titles for each column, followed
by lists of numbers holding succeeding values [t, C[t], Id[t], IaG0[t], Y[t]], in which
IaG0 is the sum Ia + G0 (autonomous investment + government spending). The value of
the "intertemporal equilibrium" Ye is calculated and printed to the screen.

maL : MAlist (0.5, 0.9, 1, 1, 2, 11);

Here is how to use the matrix function with a list of lists to make a table, letting the matrix
function do all the work of spacing the columns so the entries can be easily understood.

apply ('matrix, maL);

Inside MAlist, the first row is a list of column titles, which is entered using double quotes (")
to create Maxima strings, in the form
   ["t", "Ct", "Idt", "Ia + G0", "Yt" ].

A second function defined in Econ2.mac is  
      makeLists(L). 
This is a Maxima function which turns the list produced by MAlist into globally known 
sublists which are defined as (globally known names): tL (a list of times), 
CtL (a list of values of C[t]), IdtL (a list of values of Id[t]), IaG0L (a list of values of IaG0), 
and finally YtL (a list of values of Y[t]).

makeLists (maL);

The list names are now globally known and we can separately examine each of these lists.

tL;

CtL;

IdtL;

IaG0L;

YtL;

We now make a plot, ignoring the list IaG0L, to show values of Y[t], C[t], and Id[t] vs.
time t.

tmax : 11$
wxdraw2d(xlabel = "t", xrange = [0,tmax], yrange = [-3,16], key_pos = top_left,
    title = "a = 0.5, b = 0.9, IaG0 = 1, Y[0] = 1, Y[1] = 2",
    point_type = filled_circle,  key = "Yt",points(tL,YtL),color = red,key = "Ct",
    points (tL, CtL), color = brown,key = "Idt",points(tL,IdtL),
     color = magenta, line_width = 1,key = sconcat( "Ye = ",Ye),  explicit(Ye, t, 0,tmax),
     color = black, key = "",  explicit (0, t, 0, tmax))$

Of course we can just plot Y[t] vs. t, if we are not interested in details.

tmax : 11$
wxdraw2d (xlabel = "t", ylabel = "Y[t]", xrange = [0, tmax],yrange = [0,16],
    title = " a = 0.5, b = 0.9, IaG0 = 1,Y[0]=1,Y[1]=2",
    points (tL, YtL), line_width = 1, color = red, 
    key = sconcat( "Ye = ", Ye),  explicit(Ye, t, 0, tmax))$
    

Here is a brief example which ignores details but looks at Y[t] behavior for longer times. We
are continuing to use the same model paramters and initial conditions, but increasing tmax
from 11 to 20.

MAlist (0.5, 0.9, 1, 1, 2, 20)$
makeLists (%)$
YtL;

Let's check tL also:

tL;

tmax : 20$
wxdraw2d (xlabel = "t", ylabel = "Y[t]", xrange = [0, tmax],yrange = [0,16],
    title = " a = 0.5, b = 0.9, IaG0 = 1,Y[0]=1,Y[1]=2",
    points (tL, YtL), line_width = 1, color = red, 
    key = sconcat( "Ye = ", Ye),  explicit(Ye, t, 0, tmax))$

### Properties of the Characteristic Roots 

If r1 and r2 are the characteristic roots, r1 and r2 are always related to each other
by the following three relations: (our definitions of r1 and r2 agree with Chiang
and Wainwright (they use (b1,b2) as symbols for the characteristic roots, we use (r1,r2)).

1.)    r1 + r2 = b*(1 + a) > 0,
2.)    r1*r2 = a*b > 0, 
3.)    0 <   (1 - r1)*(1 - r2)    <  1.

Here we derive these properties using Maxima.

kill (r,a,b)$
solns : solve (r^2 - b*(1+a)*r + a*b = 0, r);

[r2,  r1] : map ('rhs, solns);

The numerator and denominator of r2 are:

num(r2);

denom(r2);

It should be clear from this definition that r2 < r1, (provided r1 and r2 are distinct real 
numbers) since for r2 we are subtracting the square root quantity from the positive 
number b*(1+a), whereas for r1 we are adding.

rsum : r1 + r2, expand, factor;

rprod : expand (r1*r2);

grind(rprod)$

rprod2 : (1 - r1)*(1 - r2), expand;

Since 0 < b < 1, 
   (1-b) > 0 and (1-b) < 1.
Hence  0 <   (1 - r1)*(1 - r2)    <  1.

### Case 1: Two Distinct Real Roots: b > 4*a/(1+a)^2

Assuming r1 and r2 are two distinct real roots, relation (2) above implies that neither
root can be equal to zero, and both r1 and r2 must have the same sign. Then relation (1)
above implies that r1 and r2 cannot both be negative, hence r1 and r2 must both be
positive: r1 > 0, r2 > 0. Hence there can be no oscillation in the time path.

### Case 1A: a < 1, 0 < r2 < r1 < 1,  convergence

For Case 1A, r1 and r2 are both real positive numbers less than one, 0 <  r2  <  r1  < 1. 
This occurs if  0 < a < 1 and b > 4*a/(1+a)^2.     

Let's make a plot which paints the area of the (a,b) plane for this real root convergence 
case in dark gray.

[xmax,ymax] : [5,3]$
wxdraw2d(xlabel = "a",ylabel = "b",xrange = [0,xmax],  yrange = [0,ymax],
    key = "1/a",  explicit (1/a,a,0.3,xmax), key = "real root convergence",
    filled_func = true, fill_color = dark_gray, explicit(1,a,0,1), color = red,
    fill_color = light_gray, key = "",  explicit(4*a/(1+a)^2, a, 0,1),
    key = "4*a/(1+a)^2",filled_func = false, explicit(4*a/(1+a)^2, a, 0,xmax),
    color = black, line_width = 1,key = "",explicit (1,a,0,xmax), parametric(1,yy,yy,0,1))$

We can numerically survey values of the roots r1 and r2 for values of b belonging to that 
part of the range 0 < b < 1 for which we get distinct real roots, for given values of a, using 
two functions: bmin(a) and survey(a) defined here.

The function survey(a) constructs a list of b values which respect b > 4*a/(1+a)^2
(so we are dealing with distinct real roots), and then evaluates r1 and r2 at each (a,b) point.

It is easier to type 1e-10 than 10^(-10), and they produce the same result in Maxima.

is(equal (1e-10, 10^(-10) ));

bmin(x) := 4*x/(1+x)^2$

survey(%a) := block ([blist, bstart,bend,delb,xx,yy,ss,%r1,%r2],
        bstart : bmin (%a) + 1e-10,
        bend : 1 - 1e-10,
        delb : (bend - bstart)/10,
        blist : makelist (bstart + n*delb,n, 0,10),      
        print (" [b,  r1 ,  r2]"),
        for %b in blist do (
           xx : %b*(1+%a),
           yy : xx^2 - 4*%a*%b,
           ss : sqrt(yy),
           %r1 : (xx+ss)/2,
           %r2 : (xx - ss)/2,
           print ([%b, %r1, %r2])),
       done)$

For example, for a = 0.5 (we increase the number of digits printed out):

fpprintprec : 16$
survey (0.5)$

We first use the more long-winded multiplier-acceleration model function yMA(a,b,IaG0,t)
method to find the indefinite solution for Y[t], taking a = 0.5 and b = 0.9, IaG0 = 1, and then
the definite solution assuming Y[0] = 20, and Y[1] = 18. (We then use the shorter path
using MAlist and makeLists).

bmin (0.5);

fpprintprec : 5$
yindef : yMA (0.5, 0.9, 1, t);

Use two single quotes (') in front of yindef to define a Maxima function ys(t).

ys(t) := ''yindef;

We take Y[0] to be equal to 20, and Y[1] = 18, and solve for the adjustable constants.

solns : solve ([ys(0) = 20, ys(1) = 18]),numer;

ydef : at (yindef, solns [1] );

yd(t) := ''ydef;

tL : makelist (t, t, 0, 11);
YtL : map ('yd, tL);

tmax : 11$
wxdraw2d (xlabel = "t", ylabel = "Y[t]", xrange = [0, tmax],yrange = [0,30],
    title = " discrete real roots, a = 0.5, b = 0.9, IaG0 = 1, Y[0] = 20, Y[1] = 18",
    points (tL, YtL), line_width = 1, color = red, key = "Ye = 10", explicit(10, t, 0, tmax))$

In this case in which a < 1, and b > 4*a/(1+a)^2, we see steady progression to the 
equilibrium value Ye = 10, sign of a stable solution.

It is quicker to use our do-loop code MAlist (a,b,IaG0,Y0,Y1, tmax) defined in Econ2.mac.

MAlist (0.5, 0.9, 1, 20, 18, 11)$
makeLists(%)$
YtL;

tmax : 11$
wxdraw2d (xlabel = "t", ylabel = "Y[t]", xrange = [0, tmax],yrange = [0,30],
    title = " discrete real roots, a = 0.5, b = 0.9, IaG0 = 1, Y[0] = 20, Y[1] = 18",
    points (tL, YtL), line_width = 1, color = red, 
    key = sconcat( "Ye = ", Ye), explicit(Ye, t, 0, tmax))$

For the details of consumption and induced investment, we use the already defined lists
CtL, IdtL.

[CtL, IdtL];

tmax : 11$
wxdraw2d(xlabel = "t", xrange = [0, tmax], yrange = [-3,30], 
    point_type = filled_circle,  key = "Yt", points(tL, YtL), color = red, key = "Ct",
    points (tL, CtL), color = brown, key = "Idt", points(tL, IdtL),
     color = magenta, line_width = 1, key = "Ye = 10",  explicit(Ye, t, 0, tmax),
     color = black, key = "",  explicit (0, t, 0, tmax))$

Nothing dramatic here.

### Case 1B:  a > 1,  r1 > r2 > 1, divergence

Case 1B assumes r2 and r1 are both numbers greater than 1. r1 > r2 > 1.

We satisfy relation (3) provided (1 - r1)*(1 - r2) [a product of two negative numbers and 
hence positive] is less than 1, so we need (1 - r1)*(1 - r2) < 1. But we have shown above 
that with 0 < b < 1, (1 - r1)*(1 - r2) = 1 - b < 1.  Relation (2) becomes  a*b = r1*r2 > 1
since both r1 and r2 are greater than 1.

Let's make a plot which paints the real root divergence region of the (a,b) plane in dark gray.

[xmax,ymax] : [5, 3]$
wxdraw2d(xlabel = "a", ylabel = "b", xrange = [0, xmax], yrange = [0, ymax],
    filled_func = true, fill_color = dark_gray, key = "real root divergence",
      explicit (1, a, 1, xmax), key = "",   fill_color = light_gray,
    explicit (4*a/(1+a)^2, a, 1, xmax), filled_func = false, 
     key = "4*a/(1+a)^2",  explicit (4*a/(1+a)^2, a, 0, xmax), color = red,key = "1/a",
    explicit (1/a, a, 0.3, xmax), color = black, line_width = 1, key = "",
    explicit (1,a, 0, xmax), parametric (1, yy, yy, 0, 1))$

As an example, we consider a = 1.5.

The value of b must then lie in the range: bmin(1.5) < b < 1.

bmin(1.5);

fpprintprec : 16$
survey (1.5)$

As an example of this divergent region of the (a,b) plane, we choose a = 1.5 and b = 0.97,
IaG0 = 1, Y[0] = 10, Y[1] = 15.

fpprintprec : 5$
MAlist (1.5, 0.97, 1, 10, 15, 7)$
makeLists(%)$
YtL;

This already shows an increase in step size each succeeding time interval, with no sign
of slowing down near 33.3.

tmax : 7$
wxdraw2d (xlabel = "t", ylabel = "Y[t]", xrange = [0, tmax],yrange = [0,200],
    title = " a = 1.5, b = 0.97, IaG0 = 1, Y[0] = 10, Y[1] = 15", key_pos = top_left,
    points (tL, YtL), line_width = 1, 
    color = red, key = sconcat( "Ye = ", Ye), explicit(Ye, t, 0, tmax))$

In this example in which a > 1, and b > 4*a/(1+a)^2, we see larger jumps each succeeding
time interval, blowing past the intertemporal equilibrium value Ye = 33.3, an example 
of "explosive instability".

For the details:

tmax : 7$
wxdraw2d(xlabel = "t", xrange = [0, tmax], yrange = [0, 250], key_pos = top_left,
    title = "a = 1.5, b = 0.97, IaG0 = 1, Y[0] = 10, Y[1] = 15",
    point_type = filled_circle,  key = "Yt",points(tL,YtL),color = red,key = "Ct",
    points (tL, CtL), color = brown,key = "Idt",points(tL,IdtL),
     color = magenta, line_width = 1,key = sconcat ( "Ye = ", Ye), explicit(Ye, t, 0, tmax))$

This unregulated economy is going gangbusters, with both consumption and induced
investment rising steadily, after we start with a rising GDP Y[t].

But with the same acceleration and propensity to consume parameters (a,b), if we start
with a large but declining GDP, it's all downhill.

MAlist (1.5, 0.97, 1, 54, 50, 7)$
makeLists (%)$
YtL;

tmax : 7$
wxdraw2d(xlabel = "t", xrange = [0, tmax], yrange = [-100, 150], 
    title = "a = 1.5, b = 0.97, IaG0 = 1, Y[0] = 54, Y[1] = 50",
    point_type = filled_circle,  key = "Yt",points(tL,YtL),color = red,key = "Ct",
    points (tL, CtL), color = brown,key = "Idt",points(tL,IdtL),
     color = magenta, line_width = 1,key = sconcat ( "Ye = ", Ye), explicit(Ye, t, 0, tmax),
     color = black, key = "",  explicit (0, t, 0, tmax))$

### Case 2: Double Roots: b = 4*a/(1+a)^2, rval(a), ydouble(a,IaG0,t)

This is the case in which b1^2 = 4*b2, or b = 4*a/(1+a)^2. From our general relations
for (r1, r2), the relation r1 + r2 = b*(1+a), and using r = r1 = r2, we get r = b*(1+a)/2.

Because a > 0 and 0 < b < 1, we have the single root r > 0, so no oscillation.

Replacing b with 4*a/(1+a)^2, we get:
        r = 2*a/(1+a).

If r < 1 then 2*a < 1 + a, and a < 1.

Likewise if r > 1 then 2*a > 1 + a, and  a > 1.

rval(_a) := 2*_a/(1 + _a)$

[rval(0.5), rval(1.5)];

We can define ydouble(a,IaG0, t) in terms of yMA(a,b,IaG0, t):

ydouble(aa,iag0,tt) := yMA(aa, bmin(aa), iag0, tt)$

### Case 2A: a < 1,  r < 1, convergence

As an example of a converging solution, we take a = 0.5, and IaG0 = 1.
We also take Y[0] = 15, Y[1] = 13.

Built into this case is b = bmin(a) = 4*a/(1+a)^2.

bmin(0.5);

yindef : ydouble (0.5, 1, t);

ys(t) := ''yindef;

solns : solve ([ys(0) = 15, ys(1) = 13]),numer;

ydef : at (yindef, solns[1]);

grind(%)$

yd(t) := ''ydef;

tmax : 10$
tL : makelist (t, t, 0, tmax);

YtL : map ('yd, tL);

wxdraw2d (xlabel = "t", ylabel = "Y[t]", xrange = [0, tmax],yrange = [0,20],
    title = " double real root, a = 0.5, b = 4*a/(1+a)^2, IaG0 = 1", 
    points (tL, YtL), line_width = 1, 
    color = red, key = "Ye = 9", explicit(9, t, 0, tmax))$

In this example with a = 0.5 < 1 we see a steady progression towards Ye = 9, with each
step smaller than the previous step.

Details, starting Y[0] above Ye (Y[0] = 15, Y[1] =13):

MAlist (0.5, bmin(0.5), 1, 15, 13, 10)$
makeLists (%)$
YtL;

tmax : 10$
wxdraw2d(xlabel = "t", xrange = [0, tmax], yrange = [-2, 25], 
    title = "a = 0.5, b = 4*a/(1+a)^2, IaG0 = 1, Y[0] = 15, Y[1] = 13",
    point_type = filled_circle,  key = "Yt",points(tL,YtL),color = red,key = "Ct",
    points (tL, CtL), color = brown,key = "Idt",points(tL,IdtL),
     color = magenta, line_width = 1,key = sconcat ( "Ye = ", Ye), explicit(Ye, t, 0,tmax),
     color = black, key = "",  explicit (0, t, 0, tmax))$

Details, starting Y[0] below Ye:

MAlist (0.5, bmin(0.5), 1, 3, 5, 10)$
makeLists (%)$
YtL;

tmax : 10$
wxdraw2d(xlabel = "t", xrange = [0, tmax], yrange = [-3, 20], key_pos = top_left,
    title = "a = 0.5, b = 4*a/(1+a)^2, IaG0 = 1, Y[0] = 3, Y[1] = 5",
    point_type = filled_circle,  key = "Yt",points(tL,YtL),color = red,key = "Ct",
    points (tL, CtL), color = brown,key = "Idt",points(tL,IdtL),
     color = magenta, line_width = 1,key = sconcat( "Ye = ", Ye), explicit(Ye, t, 0, tmax),
     color = black, key = "",  explicit (0, t, 0, tmax))$

### Case 2B:  a > 1, r > 1, divergence

MAlist (1.5, bmin(1.5), 1, 10, 13, 10)$
makeLists (%)$
YtL;

tmax : 10$
wxdraw2d(xlabel = "t", xrange = [0, tmax], yrange = [-20, 250], key_pos = top_left,
    title = "a = 1.5, b = 4*a/(1+a)^2, IaG0 = 1, Y[0] = 10, Y[1] = 13",
    point_type = filled_circle,  key = "Yt",points(tL,YtL),color = red,key = "Ct",
    points (tL, CtL), color = brown,key = "Idt",points(tL,IdtL),
     color = magenta, line_width = 1,key = sconcat ("Ye = ", Ye), explicit(Ye, t, 0, tmax),
     color = black, key = "",  explicit (0, t, 0, tmax))$

Details, with a declining GDP at the start:

MAlist (1.5, bmin(1.5), 1, 35, 32, 10)$
makeLists (%)$
YtL;

tmax : 10$
wxdraw2d(xlabel = "t", xrange = [0, tmax], yrange = [-200, 60], key_pos = bottom_left,
    title = "a = 1.5, b = 4*a/(1+a)^2, IaG0 = 1, Y[0] = 35, Y[1] = 32",
    point_type = filled_circle,  key = "Yt",points(tL,YtL),color = red,key = "Ct",
    points (tL, CtL), color = brown,key = "Idt",points(tL,IdtL),
     color = magenta, line_width = 1,key = sconcat ( "Ye = ", Ye), explicit(Ye, t, 0, tmax),
     color = black, key = "",  explicit (0, t, 0, tmax))$

### Case 3: Complex Roots: b < 4*a/(1+a)^2

### Case 3A: b < min(4*a/(1+a)^2, 1/a), sqrt(a*b) < 1,convergence

Let's make a plot to show the region in the (a,b) plane which produces this case.

wxdraw2d(xlabel = "a",ylabel = "b",key = "4*a/(1+a)^2",
    explicit (4*a/(1+a)^2,a,1,5),color = red,key = "1/a",
    explicit (1/a,a,0.5, 1), filled_func = true, fill_color = dark_gray,
    key = "complex roots, k<1",  explicit (4*a/(1+a)^2,a,0,1), 
   key = "", explicit(1/a,a,1,5),filled_func = false,
    color = blue, explicit(4*a/(1+a)^2,a,0,1),color = red,explicit(1/a,a,1,5),
    color = black,line_width = 1, explicit (1,a,0,5),parametric(1,yy,yy,0,1))$

This plot shows that for 0 < a < 1, we need 0 < b < 4*a/(1+a)^2 
and for a > 1 we need 0 < b < 1/a, to have converging k < 1 solutions arising from
complex root cases.

The analytic solution for complex roots has the form k^t * (sinusoid), where
k = sqrt (a*b), and for convergence we need k < 1.

For given values of a, survey the value of k = sqrt(a*b) for
  0 <  b < min (4*a/(1+a)^2, 1/a)
using surveyComplex1(a):

surveyComplex1(%a) := 
block ([blist, bstart,bend,delb],       
        bstart :  1e-10,
        if %a <= 0 then print("%a must be greater than zero")
        else (
           if %a < 1 then bend : 4*%a/(1+%a)^2 - 1e-10
            else bend : 1/%a - 1e-10,
           delb : (bend - bstart)/10,
   /*        print("bstart = ",bstart,"bend = ", bend," delb = ",delb), */
           blist : makelist (bstart + n*delb,n, 0,10),      
           print (" [b,  k]"),
          for %b in blist do print ( [%b, sqrt (%a*%b)] )),
       done)$

fpprintprec:16$
surveyComplex1(0.5);

Experimenting with surveyComplex1(a) shows that with
      0 < b < min(4*a/(1+a)^2, 1/a)
leads to k < 1 and a convergent solution, since the complementary solution has the
form 
      k^t*[A1*cos(th*t) + A2*sin(th*t)] 
with real adjustable constants A1 and A2.

As a complex root example showing convergence we choose a = 0.5, b = 0.5 (see plot
above), IaG0 = 1, Y[0] = 6, Y[1] = 5, in which k = sqrt(a*b) = 1/2 < 1 and Ye = 2.

fpprintprec : 5$
yindef : yMA(0.5,0.5, 1, t);

ys(t) := ''yindef;

solns : solve([ys(0) = 6, ys(1) = 5]), numer;

ydef : at (yindef, solns[1]);

yd(t) := ''ydef;

tL : makelist (t, t, 0, 8);

YtL : map ('yd, tL);

tmax : 8$
wxdraw2d (xlabel = "t", ylabel = "Y[t]", xrange = [0, tmax],yrange = [0, 8],
    title = "Complex Roots: a = 0.5, b = 0.5, IaG0 = 1, Y[0] = 6, Y[1] = 5 ", 
    points (tL, YtL), line_width = 1, color = red, key = "Ye = 2", explicit(2, t, 0, tmax))$

We clearly have convergence, with some oscillation once the curve crosses Ye = 2.

Looking at longer times:

map ('yd, makelist(t,t,0,20));

Details in this scenario with declining GDP at the start:

MAlist (0.5, 0.5, 1, 6, 5, 8)$
makeLists (%)$
YtL;

tmax : 8$
wxdraw2d(xlabel = "t", xrange = [0, tmax], yrange = [-1, 8], 
    title = " Complex roots: a = 1.5, b = 0.5, IaG0 = 1, Y[0] = 6, Y[1] = 5",
    point_type = filled_circle,  key = "Yt",points(tL,YtL),color = red,key = "Ct",
    points (tL, CtL), color = brown,key = "Idt",points(tL,IdtL),
     color = magenta, line_width = 1,key = sconcat ( "Ye = ", Ye), explicit(Ye, t, 0, tmax),
     color = black, key = "",  explicit (0, t, 0, tmax))$

Details with a rising GDP at the start:

MAlist (0.5, 0.5, 1, 0.5, 0.7, 8)$
makeLists (%)$
YtL;

tmax : 8$
wxdraw2d(xlabel = "t", xrange = [0, tmax], yrange = [-0.25, 3.5], key_pos = top_left,
    title = "Complex roots: a = 1.5, b = 0.5, IaG0 = 1, Y[0] = 0.5, Y[1] = 0.7",
    point_type = filled_circle,  key = "Yt",points(tL,YtL),color = red,key = "Ct",
    points (tL, CtL), color = brown,key = "Idt",points(tL,IdtL),
     color = magenta, line_width = 1,key = sconcat ( "Ye = ", Ye), 
     explicit(Ye, t, 0, tmax),  color = black, key = "",  explicit (0, t, 0, tmax))$

With a rising GDP at the start, we again see convergence to Ye.

For a convergent complex root case in which a > 1, we take a = 1.5 = 3/2, b = 0.5 < 1/a,
Y[0] = 6, Y[1] = 5, IaG0 = 1, Ye = 1/(1-b) = 2.

MAlist (1.5, 0.5, 1, 6, 5, 30)$
makeLists (%)$
YtL;

tmax : 30$
wxdraw2d (xlabel = "t", ylabel = "Y[t]", xrange = [0, tmax],yrange = [-1, 8],
    title = "Complex Roots: a = 1.5, b = 0.5, IaG0 = 1, Y[0] = 6, Y[1] = 5 ", 
    points (tL, YtL), line_width = 1, color = red, 
    key = sconcat ( "Ye = ", Ye), explicit(Ye, t, 0, tmax), color = black, key = "",
    explicit (0, t, 0, tmax))$

### Borderline Case: Eternal Oscillation

 Consider a borderline case in which a = 2, b = 0.5, IaG0 = 10, Y[0] = 10, Y[1] = 25.
We call this "borderline" because a > 1 and b = 1/a, on the border between convergence
and divergence. Since a*b = 1, we have k = 1, k^t = 1 for t = 0, 1, 2,.... and we get
a purely sinusoidal expression. We are taking IaG0 = 10 just for some variety.

MAlist (2, 0.5, 10, 10, 25, 40)$
makeLists (%)$
YtL;

tmax : 40$
wxdraw2d(xlabel = "t", xrange = [0, tmax], yrange = [- 10, 80], 
    title = "Complex roots: a = 2, b = 1/2, k = 1, IaG0 = 10, Y[0] = 10, Y[1] = 25",
    point_type = filled_circle,  key = "Yt",  points(tL,YtL), 
    color = magenta, line_width = 1,key = sconcat ( "Ye = ", Ye),  
     explicit(Ye, t, 0, tmax),  color = black, key = "",  explicit (0, t, 0, tmax))$

tmax : 40$
wxdraw2d(xlabel = "t", xrange = [0, tmax], yrange = [-30, 80], 
    title = "Complex roots: a = 2, b = 1/2, k = 1, IaG0 = 10, Y[0] = 10, Y[1] = 25",
    point_type = filled_circle,  key = "Yt",points(tL,YtL),color = red,key = "Ct",
    points (tL, CtL), color = brown,key = "Idt",points(tL,IdtL),
     color = magenta, line_width = 1,key = sconcat ( "Ye = ", Ye), 
     explicit(Ye, t, 0, tmax),  color = black, key = "",  explicit (0, t, 0, tmax))$

### Case 3B: a > 1, 1/a < b < 4*a/(1+a)^2, divergence

We draw a plot in the (a,b) plane showing this region of oscillating divergence 
where k  = sqrt(a*b) > 1.

wxdraw2d(xlabel = "a",ylabel = "b",key = "4*a/(1+a)^2",
    explicit (4*a/(1+a)^2,a,0,1),color = red,key = "1/a",
    explicit (1/a,a,0.5, 1),filled_func = true, fill_color = dark_gray,
    key = "complex roots, k > 1",   explicit (4*a/(1+a)^2,a,1,5), 
    fill_color = light_gray, key = "", explicit(1/a,a,1,5), filled_func = false,
    color = blue, explicit (4*a/(1+a)^2,a,1,5), color = red, explicit(1/a,a,1,5),
    color = black, line_width = 1,explicit(1,a,0,5),parametric(1,yy,yy,0,1))$

We need a different function to survey values of k = sqrt(a*b) for a divergent set (a,b).

surveyComplex2(%a) := 
block ([blist, bstart,bend,delb],       
        if %a <= 1 then print("%a must be greater than 1")
        else (
             bstart : 1/%a + 1e-10,
             bend : 4*%a/(1+%a)^2 - 1e-10,           
            delb : (bend - bstart)/10,
   /*        print("bstart = ",bstart,"bend = ", bend," delb = ",delb), */
           blist : makelist (bstart + n*delb,n, 0,10),      
           print (" [b,  k]"),
          for %b in blist do print ( [%b, sqrt (%a*%b)] )),
       done)$

We try this out for a = 1.5.

fpprintprec : 16$
surveyComplex2(1.5);

Experimenting with surveyComplex2(a) shows that with
       a > 1 and  1/a < b < 4*a/(1+a)^2 
one has a complex roots case in which k > 1, hence a divergent solution, since the 
complementary solution has the form 
      k^t*[A1*cos(th*t) + A2*sin(th*t)] 
with real adjustable constants A1 and A2.

As a complex root example showing divergence we choose a = 2, b = 0.75 (see plot
above), k = 1.225,  IaG0 = 1, Y[0] = 6, Y[1] = 5.

fpprintprec : 5$
MAlist (2, 0.75, 1, 6, 5, 30)$
makeLists (%)$
YtL;

tmax : 30$
wxdraw2d (xlabel = "t", ylabel = "Y[t]", xrange = [0, tmax],yrange = [-400, 1300],
    title = " Complex: a = 2, b = 0.75, k = 1.225 ", key_pos = top_left,
    points (tL, YtL), line_width = 1, color = red, key = sconcat ( "Ye = ", Ye), 
    explicit(Ye, t, 0, tmax), color = black,  key = "", explicit (0, t, 0, tmax))$

This diverging solution blows right past Ye = 4, with each step larger than the previous
step.

tmax : 30$
wxdraw2d(xlabel = "t", xrange = [0, tmax], yrange = [-400, 1300], key_pos = top_left,
    title = "Complex roots: a = 2, b = 0.75, k = 1.2, IaG0 = 1, Y[0] = 6, Y[1] = 5",
    point_type = filled_circle,  key = "Yt",points(tL,YtL),color = red,key = "Ct",
    points (tL, CtL), color = brown,key = "Idt",points(tL,IdtL),
     color = magenta, line_width = 1,key = sconcat ( "Ye = ", Ye),  
     explicit(Ye, t, 0, tmax),  color = black, key = "",  explicit (0, t, 0, tmax))$

which shows explosive instability.

### Summary Plot : Stable and Unstable Regions

We draw a summary plot somewhat similar to Fig. 18.2 in Chiang and Wainwright
(except that we use different colors rather than different fill styles to distinguish different
regions).

We draw two curves: b vs a for which a*b = 1, and b vs a for which b = 4*a/(1+a)^2.
In Maxima's wxdraw2d we use filled_func = true and fill_color = light_green, for example,
to paint a region beneath a function with light green color. Regions can then be painted
over by changing the choice of color (including white or gray if needed) to end up with
different colors for different regions.

The following two calls to solve allow one to design straight lines of the form
y = al*x + be which can then be used to relate labels to curves.

solve ([1.1=3.6*al + be,0.75 = 3*al + be], [al, be]), numer;

solve ([1.65=al + be, 1.3 = 0.8*al + be], [al, be]), numer;

kill(a,b)$
wxdraw2d (xlabel = "a (acceleration)",ylabel = "b (propensity to consume)", 
    yrange = [0,2], 
      color = dark_green, explicit (1/a, a, 0.5, 1), key = "Stable",
    fill_color = light_green, filled_func = true,  explicit (1,a,0,1), 
    fill_color = light_blue, key = "Unstable",explicit (1, a, 1, 5), 
    fill_color = light_salmon, key = "Explosively Unstable",  explicit (4*a/(1+a)^2, a,1, 5),
    fill_color = dark_gray, key = "Damped Stepped Fluctuation",
      explicit (4*a/(1+a)^2, a, 0, 1), key = "",  explicit (1/a, a, 1, 5),
    filled_func = false, color = black, line_width = 1,explicit (1, a, 0, 5),
    parametric (1, yy, yy, 0, 1),
    label (["b = 1/a", 1, 1.7]), label ([ "b = 4*a/(1+a)^2", 4, 1.2]),
explicit (1.75*x -0.1,x, 0.8, 1),    explicit (0.583*x -1, x, 3.05, 3.65),
    line_width = 2, color = blue, explicit (1/a, a, 1, 5),
    color = red, explicit (4*a/(1+a)^2, a, 0, 5))$
    

The stable green region has 0 < a < 1 and b > 4*a/(1+a)^2 > 0, or a*b > 4*a^2/(1+a)^2.

### Economic Instability with the Acceleration Principle

Quoting from Ch. 17, "The Business Cycle", p. 405, Economics, by Paul A. Samuelson, 
1st. ed., 1948, McGraw-Hill:

"The Acceleration Principle:
According to this law, society's needed stock of capital, whether inventory or equipment,
depends primarily upon the level of income or production. Additions to the stock of capital,
or what we customarily call net investment [induced investment], will take place only when 
income is growing. As a result, a prosperity period may come to an end - not simply because
consumption sales have gone down  - but simply because sales have leveled off at a high 
level or have continued to grow, but at a lower rate than previously."

"A simplified arithmetical example will make this clear. Imagine a typical textile-manufacturing
firm whose stock of capital equipment is always kept equal to about ten times the value of its
yearly sales of cloth. Thus, when its sales have remained at 6 million dollars per year for 
some time, its balance sheet will show 60 million dollars of capital equipment, consisting of
perhaps 20 machines of different ages with one wearing out each year and being replaced.
Because replacement balances depreciation, there is no "net" investment or saving being
done by the corporation. "Gross" investment is taking place at the rate of 3 million dollars per
year, representing the yearly replacement of one machine. (The other 3 million dollars of 
sales may be assumed to go for wages and dividends.) [Gross investment (G.I.) is the
sum of net investment (N.I.) plus replacement.]"

Let us call the first three years (sales = 6 million/year) "Phase 1", the fourth, fifth, and sixth 
years (sales = 9, 12, 15 million/yr)  "Phase 2", the seventh year "Phase 3", and the eighth 
year "Phase 4".

"... suppose that, in the fourth year, sales [of cloth]  rise by 50 per cent - from 6 to 
9 million dollars. Then the number of machines must also rise by 50 per cent, or from 20 to 
30 machines. In the fourth year, 11 machines must be bought, 10 new ones in addition to the
replacement of the worn-out one."

"Sales have gone up by 50 per cent. How much has machine production gone up? By
1000 per cent! [ the fractional increase = (11 - 1)/1,  times 100%] "

"It is this accelerated effect of a change in consumption [of cloth] on [net] investment levels
that gives the acceleration principle its name."

"If sales continue to rise in [each of] the fifth and sixth years by 3 million dollars [per year], 
then we shall continue to have 11 new machines ordered in every year [of Phase 2]."

"So far the acceleration principle has given us no trouble. On the contrary, it has given us
a tremendous increase in investment spending as a result of a moderate increase in
consumption sales. But now we are riding a tiger. Consumption has to keep increasing in
order for investment to stand still! If consumption should stop growing at so rapid a rate - 
if it should level off in the seventh year even at the high level of 15 million dollars per year -
then net investment will fall away to zero, and gross investment will fall back to 1 machine
per year (see the table). In other words, a drop of zero per cent in sales has resulted in a
90 per cent drop in gross investment and a 100 per cent drop in net investment!"

Δ(N.I.) = (0 - 30)/30 times 100% = - 100%
Δ(G.I.) = (3 - 33)/33 times 100% = - 91%

"The Lord giveth and the Lord taketh away. The acceleration principle is a two-edged sword.
If sales should drop below 15 million dollars per year, gross investment would drop away to
nothing; in fact, the firm would want to disinvest (N.I. < 0) by selling off some of its machinery
[or not replacing a 20 yr old machine no longer needed] on the used-equipment market. It is
now clear that a depression can set in just because consumption has stopped growing so
rapidly, even if it has not dropped off absolutely, but only leveled off at a high level."

"Needless to say, the curtailment of production in the machine-producing industries will 
cause them to shut down, will curtail their income and spending on food and clothing, and 
lead to still further 'multiplier' changes in spending. This might ultimately cause textile sales 
to stop growing altogether, or even to decline. This will cause a further accelerated drop in 
net investment. Thus, we may be in a vicious circle whereby the acceleration principle and 
the multiplier interact so as to produce a cumulative deflationary (or inflationary) spiral."

"It is easy to see that in the acceleration principle we have a powerful factor making for
economic instability."
....
"If business sales go up and down, the acceleration principle intensifies their fluctation. It
induces net investment on the upswing but causes about the same amount of net
disinvestment on the downswing. In the long run, if the system is growing because of 
population increase or higher real incomes, then the acceleration principle works primarily
as a stimulating factor."
...
"The economic system is more or less without a steering wheel. There is no mechanism or
automatic governor which keeps purchasing power right at the full-employment level, neither
too low nor too high. On the contrary, even if businessmen and workers are doing their best
and acting intelligently and unselfishly, the system may still be in the throes of inflation or
deflation, depending upon the chance circumstances of the complex interaction between
investment and saving.....This is the challenge to fiscal policy."

Let's make a table using the Maxima matrix function, which will worry about formatting.
Each list represents one row of the matrix, and each must have the same number of 
elements.

matrix (["time", "Yearly Sales", "Stock of Capital", "N.I.","G.I."],
    ["======","======","======","======","======"],
    ["Phase 1"," ", " "," "," "],
    ["First year", 6,60,0,3],
    ["Second year",6,60,0,3],
    ["Third year",6,60,0,3],
    ["Phase 2"," ", " "," "," "],
    ["Fourth year",9, 90, 30,33],
    ["Fifth year",12,120,30,33],
    ["Sixth year",15,150,30,33],
    ["Phase 3"," ", " "," "," "],
    ["Seventh year",15,150,0,3],
    ["Phase 4"," ", " "," "," "],
    ["Eighth year", 14.7, 147, -3,-3]);