# 2nd Order Differential Equations 
Dowling18A.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 ("Econ2.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. 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 slightly changed some of the symbols used in particular problems.

Ted Woollett
https://web.csulb.edu/~woollett/
Oct. 26, 2021

## References

Fundamental Methods of Mathematical Economics, Fourth Ed., A. C. Chiang, and
K. Wainwright, McGraw-Hill, 2005.

## Fast Track, ode2(eqn,y,t), check (y, t)

Fast track linear 2nd order ode with constant coefficients.
Standard form:  y'' + A y' + B y(t) = C.
For ease of checking solution, write eqn with C brought over to the left hand side with a
minus sign:
eqn :  'diff (y,t,2) + A*'diff(y,t) + B*y - C;
Our Maxima function in Econ2.mac is:
Lode2(y,t,type,A,B,C), with type = real, complex, or double, depends on A and B:

The characteristic roots are solutions of r^2 + A*r + B = 0.
distinct real roots if A^2 > 4*B
complex roots if A^2 < 4*B
double real root if A^2 = 4*B.

eqn : 'diff(y, t, 2) + 9*'diff (y,t) + 14*y - 7;

For numerical coefficients A, B, Maxima's ode2 returns a solution without asking questions.
If A and/or B are symbolic, ode2 asks questions (see examples later).

Msoln : ode2 (eqn,y,t);

Solution assuming values of y and y' at t = 0:

Mpsoln : ic2(Msoln,t=0,y=10,'diff(y,t) = -2);

Using our Lode2 function (Econ2.mac) you must say what kind of roots you expect.

Lsoln : Lsoln : Lode2(y,t,real,9,14,7);

Lpsoln : ic2(Lsoln,t=0,y=10,'diff(y,t) = -2);

A solution returned by the Maxima function ode2 ( we used ode2 for first order ordinary
differential equations in Dowling16.wxmx) should never need to be checked for accuracy
(since why use Maxima if it is not always reliable), but just to show off Maxima tricks, we will 
define check(y,t) here, which should produce zero (0) if y is replaced by an actual solution 
as a function of t (given our way of defining eqn).

If you have defined eqn with the constant C term brought over to the lhs with a minus
sign (as above), this function should produce 0.

check(y, t) := expand (ev (eqn,diff));

ys : rhs (Lpsoln);

check(ys,t);

What if you made the wrong type assumption for a numerical case?

soln : Lode2(y,t,complex,9,14,7);

## Numerical Initial Value Problems, ode2, ic2

### Dowling Prob. 18.8

Given the 2nd order ODE: y'' + 9 y' + 14 y = 7, find the particular (definite) solution for which
y(0) = -5/2 and y'(0) = 31. Comment on stability.

This example illustrates the required order of arguments for the Maxima functions ode2 & ic2.
Notice the use of the single quote mark (') in front of diff, to prevent immediate action by
the Maxima function diff, in the definition of eqn.

eqn : 'diff(y, t, 2) + 9*'diff (y,t) + 14*y - 7;
gsoln : ode2 (eqn, y, t);
psoln : ic2 (gsoln, t = 0, y = -5/2, 'diff (y,t) = 31 );

grind(%)$

check (rhs (psoln),t);

### yc, yp

If the ODE has the form y'' + A y' + B y = C, the "complimentary solution" yc is
yc = A1 exp (r1*t) + A2 exp (r2*t), in which r1, r2 are the roots of the quadratic equation
r^2 + A r + B = 0. This equation comes from using for y(t) the trial solution
  y = exp(r*t) in the "reduced equation" gotten by setting C = 0:
      y'' + A y' + B y = 0, whence
  (r^2 + A*r + B)*exp(r*t) = 0, but exp (r*t) # 0, so we must have r^2 + A*r + B = 0 for this
to work.

For this example, these roots are:

solns : solve (r^2 + 9*r + 14, r);

The general solution (with the two arbitrary adjustable constants %k1 and %k2) has the 
structure of ygen = yc + yp, in which yp is the "particular solution" yp = C/B, which we can 
call the "equilibrium solution" (hopefully) since if we assume we have waited so long that 
both the first and second derivatives of y(t) equal zero, the given ODE reduces to B*y = C, 
or y(t) = C/B = independent of time. In our example, C/B = 7/14 = 1/2.

Since both "characteristic roots" r1 and r2 are negative numbers, and exp(-|r|*t) --> 0 as
t >>1/|r|, the general solution will approach the particular solution yp = 1/2 if we wait long 
enough (t >> 1/7 = 0.143). "much greater" (>>)  means, say, ten or twenty times, 
so t > 1.4 as a rough estimate.

ys : rhs (psoln);

wxdraw2d (xlabel = "t", ylabel = "y",key_pos = bottom_right,
    key = "y(0) = -5/2, dy/dt(0) = 31", explicit (ys, t, 0, 2), 
    color = red, key = "ye = 1/2", explicit (1/2, t, 0, 2), 
    color = black, key = "", line_width = 1, explicit (0, t, 0, 2))$

### Characteristic Roots Notation: g, h

Maxima returns two solutions of the equation r^2 + A r + B = 0:

solve (r^2 + A*r + B, r);

We adopt a notation based on Dowling's treatment:
  r1, r2 = g  +/- i h,  or in Maxima notation, g +/- %i*h, where %i = sqrt(-1), (%i)^2 = -1,
and g = - A/2, h  = sqrt(4*B - A^2)/2, connected with our chosen notation
for a general linear 2nd order ODE:  y'' + A y' + B y(t) = C.

For this example, A = 9, B = 14, so g = -9/2, and
h = sqrt(4*B - A^2)/2 = sqrt(56 - 81)/2 = %i*sqrt(81 - 56)/2 = %i*(5/2),
so r1 = g + %i*h = -9/2 + (%i)^2*5/2 = -9/2 - 5/2 = -14/2 = -7,
and r2 = g - %i*h = -9/2 + 5/2 = -4/2 = -2.

Since both r1 and r2 are negative (real) numbers, the solution y(t) will approach a
stable solution for "large" positive values of t.

### croots (A,B)

Let croots (A,B) be a Maxima function which returns a list of characteristic roots
in the form [r = r1, r = r2].

croots(%A,%B) := 
block ([gg,hh],
    gg : -%A/2,
    hh : sqrt (4*%B - %A^2)/2,
    [r = gg + %i*hh, r = gg - %i*hh])$

croots (9,14);

We can get a list of the right hand sides of this list using map('rhs, alist).

map ('rhs, %);

We can define r1 and r2 via a list as follows:

[r1, r2] : map('rhs, croots(9,14));

Try interrogating this list:

[r1,r2];

Try interrogating just r1:

r1;

### r1r2 (A, B)

Let's define a second Maxima function that uses this mapping to define r1 and r2.

r1r2(AA, BB) := map ('rhs, croots (AA, BB))$

To check that this function works, we first remove Maxima's knowledge of our prior
assignment of values to the symbols r1 and r2, by using the Maxima function kill.

kill(r1,r2)$

Now Maxima should not know any value for r1 and r2.

[r1,r2];

Now use our new Maxima function r1r2 (A, B) to define r1 and r2:

[r1,r2] : r1r2 (9, 14);

r1;

### Dowling Prob 18.9  

Find the definite integral of the 2nd order ODE
 y'' - 4 y' - 5 y = 35
such that y(0) = 5 and y'(0) = 6.

The characteristic roots r1,r2 of the characteristic equation
r^2 + A r + B = 0, or r^2 - 4 r - 5 = 0, are given by

solve (r^2 - 4*r - 5, r);

Or, we can use our Maxima function r1r2 (A, B):

[r1,r2] : r1r2 (-4, -5);

Since we have one positive root, that will lead to explosive increase with time, and 
NOT approach  ye = 35/(-5) = -7 for "large" times.

Reverting to the more general approach, we use the Maxima function ode2.
The symbol eqn is defined here using 'diff (...), putting a single quote (') mark in front of
the Maxima function diff, which prevents diff from executing immediately.
gsoln is defined to be the general solution with two arbitrary constants %k1 and %k2,
and dsoln is the "definite" solution corresponding to the t = 0 values of y and y'.

eqn : 'diff(y, t, 2) - 4*'diff (y,t) - 5*y = 35;
gsoln : ode2 (eqn, y, t);
dsoln : ic2 (gsoln, t = 0, y = 5, 'diff (y,t) = 6 );

Let ys be a Maxima expression, depending on the time t, which is the desired
solution for y(t).

ys : rhs (dsoln);

wxdraw2d (xlabel = "t", ylabel = "y", yrange = [-10, 35],
    key_pos = top_left,
    key = "y(0) = 5, dy/dt(0) = 6", explicit (ys, t, 0, 0.5),
    color = red, key = "ye = -7",explicit (-7, t, 0, 0.5),
    color = black, key = "", line_width = 1, explicit (0,t,0,0.5))$

which shows explosive (exponential) growth.

## y''(t) + A y'(t) + B y(t) = C,  General Solution from ode2

Consider the "2nd ODE with constant coefficients"
     d²y/dt² + A dy/dt + B y(t) = C,
in which we assume A, B, and C are each real constants.

We define eqn as an expression (which equals zero for later convenience) in a form
which ode2 can understand. The single quote mark (') in front of diff prevents diff from
going through the actions normally taken by the diff function (technically, the single
quote mark "prevents evaluation", producing a "noun form" rather than a "verb form").

eqn : 'diff (y,t,2) + A*'diff(y,t) + B*y - C;

A "general solution" y(t) = yc(t) + yp.

The "particular solution" yp  is proportional to C and may or may not depend on t.

If C = 0, then yp = 0.

If C # 0 then 

  1.) if B # 0, then yp = C/B, which is gotten by assuming y has ceased changing with
time in the sense that both the first and second derivative of y wrt t are zero. In that case
the given ODE reduces to B*y = C.

2.) if B = 0, and A # 0 then yp = C*t/A, which is gotten by assuming y'' = 0, in which
case the ODE reduces to dy/dt = C/A, and the indefinite integral is y = C*t/A.

3.) if B = 0 and A = 0 then yp = C*t^2/2, since the ODE reduces to y'' = C, the first
indefinite integral is y'  = C*t, and the indefinite integral of y' is C*t^2/2.

The "complementary solution" yc = A1*exp(t1*t) + A2*exp(r2*t), in which r1 and r2 are
the solutions of the quadratic equation r^2 + A*r + B = 0, and A1 and A2 are constants
which can be chosen to fit initial conditions (for example).

solns : solve (r^2 + A*r + B, r);

CASE A^2 > 4*B

If A^2 > 4*B, r1 and r2 are real roots and depending on the signs of r1 and/or r2 there is a 
possibility that the complementary solution can grow exponentially (vide infra). 

Let s = sqrt(A^2 - 4*B), then r1 = -(s + A)/2, and if (s + A) < 0 [which is equivalent to s < - A] 
then r1 > 0 and exp(r1*t) grows exponentially. If  (s + A) > 0 [ equivalent to s > - A], then 
r1 < 0 and exp(r1*t) decays exponentially to zero.
--------------------------------------------------------------------------
CASE   4*B > A^2

If, on the other hand, 4*B > A^2, then r1 and r2 are complex roots.

Let p = sqrt(4*B - A^2) and use sqrt (-1) = %i and sqrt (-x) = %i*sqrt(x). Then we get 
r1 = -A/2 - sqrt(- (4*B - A^2))/2 = -A/2 -%i*sqrt(4*B - A^2)/2 = -A/2 -%i*p/2.

Similarly, r2 = -A/2 +%i*p/2. 

Then yc = exp(-A*t/2)*(A1*exp(-%i*p*t/2) + A2*exp(%i*p*t/2)).

Use exp(%i*θ) = cos(θ) + %i*sin(θ), and exp(-%i*θ) = cos(θ) - %i*sin(θ), to replace the
complex exponentials with trig functions, replace A1 and A2 by a different pair of constants
k1 and k2 to get exp(-A*t/2) times a linear combination of trig functions with the
argument p*t/2. 

If A > 0, the overall factor exp(-A*t/2) drives the whole term to zero for t >> 1/A.
If A < 0, the factor exp(-A*t/2) shows explosive growth. 
If A = 0, yc is an eternally oscillating function of time within fixed limits.

### Case:   A^2  <  4*B,   Complex Roots,   Lode2(y,t,type,A,B,C)

This is the simpler case in which the characteristic roots r1 and r2 are complex, and the
nature of the complementary solution depends crucially on the sign of A. The Maxima
function ode2 asks us immediately about the sign of (4*B - A^2).

gsoln : ode2 (eqn,y,t);

In this case in which (4*B - A^2) > 0, we get a solution which approaches C/B as long
as A > 0, so that the exponential factor exp( -A*t/2) drives the oscillating sine and cosine
term to zero for t >> 1/A.

Let ypos be an expression depending on t which is a time dependent solution
depending on the constant real values of A, B, and C.

ypos : rhs(gsoln);

Here we make sure (OCD syndrome) this solution satisfies the given ODE.

check (ypos, t);

An alternative approach (in Maxima) is to design a Maxima function which requires us to
decide (ahead of time) on the relative sizes of A,B,C (or express them symbolically)
and ask for the specific cases: complex roots (sines and cosines) if A^2 < 4*B, distinct
real roots (if A^2 > 4*B), or double real roots (if A^2 = 4*B). For the arg 'type' in Lode2
then use either complex, real, or double.

We showed an example of using Lode2 already in the "Fast Track" section at the top.
You need to load Econ2.mac first (which we did at the top of this worksheet).

Lsoln : Lode2 (y, t, complex, A, B, C);

### Example: y'' + A y' + y = 0, A^2 < 4

We assume here that C = 0, B = 1, and %k1 = %k2 = 1.
Since C = 0, yp = 0. Also, since we are assuming (above) that  4*B > A^2, the
roots are complex and we get sines and cosines modulated by an exponential factor.
B = 1, together with 4*B > A^2, implies that 0 < A < 2 or -2 < A < 0.

We first use our Maxima function Lode2 and specialize to C = 0, B = 1.

Lode2(y,t,complex,A,1,0);

A global flag dofactor is set to true in Econ2.mac

dofactor;

If we change that setting of dofactor here, we get an unfactored argument of the
square root function sqrt.

dofactor : false$
Lode2(y,t,complex,A,1,0);

grind(%)$

Next we use Maxima's solution from ode2 (gsoln and ypos above) and specialize to the
case C = 0, B = 1, %k1 = %k2 = 1. The name "ypos" reflects our assumption that 
4*B > A^2 in this section, corresponding to the question ode2 asks.

We can use the at (...) function to define yex for the case C = 0, B = 1, %k1 = %k2 = 1.

yex : at (ypos, [C = 0, B = 1, %k1 = 1,%k2 = 1]);

In economics, we are interested (usually) in positive times (the future). What is the value of
yex at time t = 0? Again we can use the at function.

at(yex, t = 0);

So if we vary the parameter A, this solution always starts at y(0) = 1.

We now define a Maxima function depending on (t, A) using TWO single quotes (') in front
of yex (to force evaluation of the expression yex):

ys (t, A) := ''yex;

The t = 0 value of y(t), for any real value of A (satisfying A^2 < 4) can then be found using
this function ys (t, A).

ys(0,A);

The Maxima function ys can now also be used to show plots for different values of A.

Positive values of A, respecting (4*B - A^2) > 0, or A^2 < 4 here, show ys(t) converging
to C/B  = 0, as seen in this plot. The overall factor exp (-A*t/2), which starts out with unit
value at t = 0, approaches zero rapidly for "large" positive values of t, (t >> 1/A).

wxdraw2d (title =  " y’’ + A y’ + y = 0, 0 < A < 2", xlabel = "t",ylabel = "y", 
   key = "A = 0.5", explicit (ys(t, 0.5), t, 0, 20),
   key = "A = 0.25",  color = brown, explicit (ys(t,0.25),t,0,20), 
   key = "A = 0.125", color = red, explicit ( ys (t, 0.125),t, 0, 20),
   color = black, key = "", line_width = 1, explicit (0,t,0,20))$

Larger positive values of A lead to faster convergence to the equilibrium value y = 0.

Next we look at negative values of A. In this case the overall factor exp (-A*t/2) increases
explosively away from unit value (at t = 0) for positive times t.
Negative values of A, respecting (4*B - A^2) > 0, or A^2 < 4 here, show ys(t) diverging
(in absolute value) away from y(0) = 1.

tmax : 15;
wxdraw2d (title =  " y’’ + A y’ + y = 0,  -2 < A < 0", 
    xlabel = "t",ylabel = "y", yrange = [-20,20],  key_pos = bottom_left,
    key = "A = -0.5", explicit (ys(t, -0.5), t, 0, tmax),
   key = "A = -0.25",  color = brown, explicit (ys(t,-0.25),t,0,tmax), 
   key = "A = -0.125", color = red, explicit ( ys (t, -0.125),t, 0, tmax),
   color = dark_cyan, key = "y(0) = 1", explicit (1,t,0,0.75),
    color = black, key = "", line_width = 1, explicit (0,t,0,tmax))$

Each case of negative A shows divergence from y = 0. The greater the magnitude of A, the
faster the divergence.

### Example: y'' +  y' + B y = 0,  B > 1/4

A second example of a linear second order ODE with constant coefficients and term, which
has complex characteristic roots is y'' + y' + B y = 0, as long as B > 1/4.

Lode2 (y,t,complex,1,B,0);

Again, for practice using results from the chain starting with ode2, we redefine yex by
specializing ypos (defined above) for this case, including unit values for %k1 and %k2.

yex : at (ypos, [C = 0, A = 1, %k1 = 1,%k2 = 1]);

yex;

This is clearly damped for positive t >> 1.
We again use two single quote marks (') in front of the Maxima expression yex to define
a Maxima function ysf (t, B).

ysf (t,B) := ''yex;

ysf(t, B);

The initial value (t = 0) for any B > 1/4 is found from this function using ys(0,B):

ysf(0,B);

tmax : 10$
wxdraw2d ( xlabel = "t", ylabel = "y",yrange = [-1, 1.5],
    title =  " y’’ + y’ + B*y = 0,  B > 1/4",    
   key = "B = 0.3", explicit (ysf(t, 0.3), t, 0, tmax),
key = "B = 0.5", color = red, explicit (ysf(t, 0.5), t, 0, tmax),
    key = "B = 1", color = dark_cyan, explicit (ysf(t, 1), t, 0, tmax),
    key = "B = 5", color = brown, explicit (ysf(t, 5), t, 0, tmax),
    color = black,key = "", line_width = 1, explicit (0,t,0,tmax))$

Large positive values of B produce clear damped oscillations, small values of B > 1/4 are
also damped, but without much oscillation.

### Case:   A^2  > 4*B,   Distinct Real Roots

This is the more complicated case, in which the characteristic roots r1 and r2 are real.
Here (A^2 - 4*B) > 0 and the solution is expressed in terms of sqrt(A^2 - 4*B).

Lode2 (y, t, real, A, B, C);

gsoln : ode2 (eqn,y,t);

yneg : rhs(gsoln);

check (yneg, t);

### Example: y'' + A y' + y = 0, A^2 > 4

Lode2 (y,t,real,A,1,0);

C = 0 implies that yp = 0, and y = yc.
B = 1, together with A^2 > 4*B, implies that  A > 2 or  A < -2.
We also assume %k1 = %k2 = 1.

yex : at (yneg, [C = 0, B = 1, %k1 = 1,%k2 = 1]);

at (yex, t = 0);

Define a Maxima function ysf(t,A) which can be used for A either positive or negative,
as long as A^2 > 4.

kill(ysf)$
ysf (t, A) := ''yex;

We first consider the case that A > 2, so A is positive.

Since sqrt(A^2 - 4) < A , the first exponential terms has the form exp(-r*t) with r > 0, t > 0,
which rapidly becomes vanishingly small as t gets "large".  The second exponential term
also has the form exp (-r*t) with both r and t positive and goes to zero even faster than the
first term. For longer times the first term dominates and is zero for t >> 1/(A - sqrt(A^2 - 4)),

Initial value of y at t = 0, for any A > 2:

ysf(0,A);

wxdraw2d (xlabel = "t", ylabel = "y", yrange = [0.2, 2],
    title =  " y’’ + A y’ + y = 0,  A > 2", 
    key = "A = 3", explicit (ysf (t, 3), t, 0, 3),
    color = brown, key = "A = 6", explicit (ysf (t, 6), t, 0, 3),
     color = red, key = "A = 12",explicit (ysf (t, 12), t, 0, 3))$

We see convergence to C/B = 0 as long as A > 2, with faster convergence for values of A 
closer to 2. ysf(t) gets close to C/B for t >> 1/ (A - sqrt(A^2 - 4)).

We next consider negative values of A, respecting A < -2 (so that A^2 > 4).
The first term of ysf (t,A) is exp(r*t) with r = (sqrt(A^2-4) - A) = (sqrt(A^2-4) + |A|) > 0.
The second term of ysf(t,A) is exp(r*t) with r = -sqrt(A^2-4) - A = (|A| - sqrt(A^2-4))> 0.
Both terms grow explosively.

tmax : 0.5;
wxdraw2d (xlabel = "t", ylabel = "y", yrange = [0, 20], 
    title =  " y’’ + A y’ + y = 0,  A < - 2 ", 
    key_pos = top_left,
    key = "A = -3", explicit (ysf(t, -3), t, 0, tmax),    
    color = brown, key = "A = -6", explicit (ysf (t, -6), t, 0, tmax),
     color = red, key = "A = - 9",explicit (ysf(t, -9), t, 0, tmax),
    color = dark_cyan, key = "y(0) = 2", explicit (2, t, 0, tmax))$

We see exponential growth in y(t) away from y(0) = 2, with faster growth occurring with more
negative values of A < -2 (greater in magnitude).

### Example: y'' + y' + B y = 0, B < 1/4

This is a second example of distinct real roots.

C = 0 implies that yp = 0, and y = yc.
A = 1, together with A^2 > 4*B, implies that 1 > 4*B, or   B < 1/4.
We also assume %k1 = %k2 = 1.

Lode2 (y, t, real, 1, B, 0);

yex : at (yneg, [C = 0, A = 1, %k1 = 1,%k2 = 1]);

at (yex, t = 0);

kill(ysf)$
ysf (t, B) := ''yex;

ys(0,B);

tmax : 10$
wxdraw2d (xlabel = "t", ylabel = "y", yrange = [0, 2.1], 
    key_pos = bottom_left,
    title =  " y’’ + y’  + B y = 0,  B < 1/4", 
    key = "B = 0.1", explicit (ysf (t, 0.1), t, 0, tmax),
    color = dark_cyan, key = "B = -0.1", explicit (ysf (t, -0.1), t, 0, tmax),
    color = brown, key = "B = -0.2", explicit (ysf (t, -0.2), t, 0, tmax),
     color = red, key = "B = -0.3",explicit (ysf (t, -0.3), t, 0, tmax))$

The case B = 0.1 illustrates the result that if 0 < B < 1/4, ysf(t) has no minimum value.
The other three cases have B negative, and we see that ysf(t) has a minimum value for 
some t.

### Case:   4 B = A^2, Double Real Root,  Examples A > 0, A < 0

When 4B = A^2, the solution has both linear and exponential factors, which sometimes 
leads to a curve with a point of inflection. With B = A^2/4, the characteristic equation
becomes
  r^2 + A r  + A^2/4 = 0, or (r + A/2)^2 = 0, or r = -A/2.

factor (r^2 + A*r + A^2/4);

solve (r^2 + A*r + A^2/4, r);

Lode2(y,t,double,A,A^2/4,C);

gsoln : ode2 (eqn, y, t);

yzero : rhs(gsoln);

yex : at (yzero, [C = 0, %k1 = 1, %k2 = 1]);

kill(ysf)$
ysf (t, A) := ''yex;

tmax : 10$
wxdraw2d (xlabel = "t", ylabel = "y", yrange = [0, 1.4],
    title =  " y’’ + A y’ + (A^2/4) y = 0,  A > 0 ", 
    key = "A = 1", explicit (ysf (t, 1), t, 0, tmax),
    color = brown, key = "A = 2", explicit (ysf (t, 2), t, 0, tmax),
     color = red, key = "A = 4",explicit (ysf (t, 4), t, 0, tmax))$

A numerical example in which A = -4, and we vary the relative
sizes of the constants %k1, %k2:

yex : at (yzero, [C = 0, A = -4]);

kill(ysf)$
ysf(t, %k1,%k2) := ''yex;

We now consider positive and negative combinations of %k1 and %k2 over the t interval [0,1].

wxdraw2d (xlabel = "t", ylabel = "y", yrange = [-15, 15],
    title =  " y’’ - 4 y’ + 4y = 0 ", 
    key_pos = top_left,
    key = "k1=1, k2 = 1", explicit (ysf (t,1, 1), t, 0, 1),
    color = brown, key = "k1=1, k2 = -1", explicit (ysf (t,1, -1), t, 0, 1),
     color = red, key = "k1= -1, k2 = 1", explicit (ysf (t, -1, 1), t, 0, 1),
      color = magenta, key = "k1= -1, k2 = -1", explicit (ysf (t, -1, -1), t, 0, 1),
    color = black, key = "", line_width = 1, explicit (0, t, 0, 1))$

### Case:  B = 0, Dowling 18.4

Consider the equation 
   y'' + 7 y' = 28
for which A = 7, B = 0, C = 28, A^2 = 49 > 4*B = 0.

eqn : 'diff(y,t,2) + 7*'diff(y,t) -28;
gsoln : ode2(eqn,y,t);

If A or B is zero, Lode2 automatically calls ode2.

Lode2(y,t,real,7,0,28);

expand(%);

Naturally, one can combine %k1 and the term -4/7 into some other adjustable constant, 
say A1, and write this as y =  A1 + 4 t + A2 exp(- 7*t),  as in Dowling.

## Market Model with Price Expectations

We adapt Sec. 16.4 in Chiang and Wainwright.

In prior chapters the quantity demanded Qd and the quantity Qs supplied were assumed to
depend only on the current price P. We can relax that assumption by recognising that
both buyers and sellers often are influenced by the "price trend" prevailing at a given
time t. The current price trend can induce expectations about the price level in the future, 
and affect either (or both) Qd and Qs.

In the continuous-time context, price trend information comes mainly from the current value
of both the current rate of change of the price, P' =  dP/dt, and the current rate of change of 
dP/dt, P'' = d²P(t)/dt².

A continuous time market model including the influence of the current price trend (in this
context) would  require a model of the form:
    Qd = D[P(t), P'(t), P''(t)],
    Qs = S[P(t), P'(t), P''(t)].

If we specialize to the linear version of these functions D and S, expressed in terms of
constant coefficients b,m,n,d,u,w and terms a,g:
   Qd = a - b P + m P' + n P'', with a,b > 0,
   Qs = -g + d P + u P' + w P'',  with g,d > 0.

The price trend parameters, whose signs have not been restricted, depend on the future
price expectations of the buyers (m and n) and the sellers (u and w).

If m > 0, a rising price (P' >0) will result in an increase of Qd - buyers expect the price to
continue to rise significantly and opt to make purchases at the current time t, while the
price is affordable. The magnitude of m, in both sign cases, would reflect the strength of 
the corresponding expectations.

If m < 0, a rising price (P' > 0) will result in a decrease of Qd - buyers expect a quick
reversal of the price rise behavior and decide to stop current purchases and wait for a
lower price level to appear.

The inclusion of the parameter n makes the buyers' behavior depend also on the rate of
change of dP/dt, and both parameters m and n include the effects of price speculation on
the behavior of the market price in this model.

The signs and magnitudes of the supply side parameters u and w likewise would depend
on sellers' expectations of future price levels.

### A Simplified Model

Assume here that we only need to include price trend assumptions in the demand function
expression Qd:
    Qd = a - b P + m P' + n P'', with a,b > 0,
    Qs = -g + d P,  with g,d > 0
Assume the market is cleared at every point of time (continuous time model):
    Qd = Qs, or Qd - Qs = 0.

As a help for algebraic manipulation let P1 stand for P' and P2 stand for P''.

kill(a,b,m,n,g,d)$
Qd : a - b*P + m*P1 + n*P2;
Qs : - g + d*P;
aeqn : Qd - Qs;

We want to arrive at the 2nd order ODE form
   P'' + A P' + B = C, or  P'' + A P' + B - C = 0
which we adopted in prior sections. So divide by the parameter n.

factorsum(aeqn);

%/n;

Dividing by n we get the form 
  P'' + A P' + B P - C = 0
with
A = m/n,      B = - (b + d)/n,  and    C = -(a + g)/n.

When P'' = 0 and P' = 0, we reach "intertemporal equilibrium" Pe = C/B (perhaps only
for a moment, if this is an unstable equilibrium),
      
Pe = (a+g)/(b+d) > 0.

We get two distinct real roots if A^2 > 4*B, or if (m/n)^2 > - 4*(b+d)/n.

If n > 0, -4*(b+d)/n is a negative number (since b and d are positive), and we certainly get
two distinct real roots. What about stability?

Without assuming a sign for n, Maxima's solve function gets the roots r1 and r2 
in the form:

assume(b > 0, d > 0)$
[r1, r2] : map ('rhs, solve (r^2 + (m/n)*r - (b+d)/n = 0,r) );

If we assume n > 0, the expression under the radical is both positive and greater than m^2;   
[4*(b+d)*n + m^2 ] > m^2, so its square root is greater than the magnitude of m, so r1 < 0 
and r2 > 0, no matter the sign of m. So n > 0 only leads to distinct real roots.

With one negative root and one positive root, there is no stable intertemporal equilibrium 
possible (except for the case that the initial conditions exclude the effect of the positive root).

If we assume n < 0, then all three cases become possible:
   real roots if m^2 >  4*(b+d)*|n|, in which case 0 < sqrt(4*(b+d)*n + m^2)) < |m|


   complex roots if m^2 < 4*(b+d)*|n|.
   

We get complex roots r1 and r2 if A^2 < 4*B, or if  (m/n)^2 < -4*(b+d)/n, or (multiplying 
both sides by n^2) if m^2 < -4*(b+d)*n, or (m^2 + 4*(b+d)*n) < 0,  which is not possible 
 if n is positive; we need n < 0 to get complex roots.

### Distinct Real Roots Case:    (m/n)^2 > -4*(b+d)/n

Distinct real roots occur if A^2 > 4*B (from our previous discussion), or
if (m/n)^2 > -4*(b+d)/n, which is ALWAYS true if n > 0, and MAY be true if n < 0.
In this case P(t) = A1*exp(r1*t) + A2*exp(r2*t) + Pe, with the signs and magnitudes of
the real roots r1 and r2 dependent of the values of A and B.

The distinct real roots r1, r2 are the roots of the equation  r^2 + A*r + B = 0.

solve (r^2 + (m/n)*r -(b+d)/n = 0,r);

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

yreal : %k1*exp(r1*t) + %k2*exp(r2*t) + (a+g)/(b+d);

### Example 1

Consider the numerical example
  Qd = 42 - 4 P - 4 P' + P''
  Qs = -6 + 8 P
with the initial conditions P(0) = 6, P'(0) = 4.

case1 : [a=42,b=4,m=-4,n=1,g=6,d=8]$

at ([r1, r2], case1);

yex1 : at (yreal, case1);

ic2(P = yex1,t = 0, P= 6, 'diff(P,t) = 4);

The existence of the positive root r2 = 6 means that the intertemporal equilibrium Pe = 4
is dynamically unstable.

y1 : rhs(%);

at (y1, t = 0);

wxdraw2d ( xlabel = "t", ylabel = "P", yrange = [0, 100],
    title = "P = exp(6*t) + exp(-2*t) + 4",
    explicit (y1,t,0,1),color = red, line_width = 1,
    key_pos = top_left, key = "P(0) = 6", explicit (6,t,0,0.2))$

### Complex Roots Case:  (m/n)^2 < -4*(b+d)/n

Assume the demand and supply functions are
  Qd = 40 - 2 P - 2 P' - P'' = a - b P + m P' + n P'', m/n = (-2)/(-1) = 2, A^2 = 4.
  Qs = -5 + 3 P = -g + d P, 
b = 2, d = 3, b+d = 5, n = -1, B = -(b+d)/n = -5/(-1) = 5, A^2 < B => complex roots,
Easiest to use ode2 directly.
Equate Qd and Qs and divide by n to get
     P'' + (m/n) P' - (b + d) P/n = - (a + g)/n
         P'' + 2 P' + 5 P = 45.

eqn : 'diff (P,t,2) + 2*'diff(P,t) + 5*P = 45;
gsoln : ode2(eqn, P, t);

For t >> 1 (say 5 or 10), P ~ 9.

### Example 2

Assume P(0) = 12, P'(0) = 1.

dsoln : ic2 (gsoln, t = 0, P = 12, 'diff (P,t) = 1);

y2 : rhs (dsoln);

at (y2, t = 0);

wxdraw2d (xlabel = "t", ylabel = "P", yrange = [8,12.5],
    explicit (y2, t, 0, 5), line_width = 1, color = red, key = "Pe = 9",
    explicit (9, t, 0, 5))$

### Including Supply Side Price Trend Expectations

Assume
  Qd = 9 - P + P' + 3 P'',
  Qs = -1 + 4 P - P' + 5 P'',
Equating Qd and Qs and writing in form P'' + A P' + B P = C, 
we get P'' - P' + (5/2) P = 5, A = -1, B = 5/2, C = 5, Pe = C/B = 2, A^2 = 1, 4*B = 10,
A^2 < 4*B =>  complex roots.

eqn: 'diff(P,t,2) - 'diff (P,t) + (5/2)*P = 5;
gsoln : ode2 (eqn, P, t);

We see a divergent price P(t) predicted by this model.
The intertemporal equilibrium  Pe = 2 is dynamically unstable.
Assume P(0) = 4 and P'(0) = 4.

dsoln : ic2(gsoln, t = 0, P = 4, 'diff (P,t) = 4);

y3 : rhs (dsoln);

at (y3, t = 0);

wxdraw2d (xlabel = "t", ylabel = "P", explicit (y3, t, 0, 2),
    color = black, line_width = 1, explicit (0, t, 0, 2))$

solve (diff (y3,t) = 0, t);

An extremum of P(t) thus occurs at tan(3*t/2) = 2, or t = (2/3)*arcTan(2),

tmax : 2.0*atan(2.0)/3;

Pmax : at (y3, t = tmax);

wxdraw2d (xlabel = "t", ylabel = "P", yrange = [-2, 6],
    explicit (y3, t, 0, 2), color = black, line_width = 1, 
    explicit (0, t, 0, 2),color = red, key = "tmax = 0.738",
    parametric (tmax, yy, yy, -2, Pmax))$

## Dowling Prob 18.36

In many markets supply and demand are influenced by current prices and price trends
 (i.e., whether prices are rising or falling and whether they are rising or falling at an 
increasing or decreasing rate). The economist, therefore, needs to know the current 
price P(t), the first derivative dP(t)/dt, and the second derivative d²P(t)/dt². 

Assume the quantity supplied Qs(t) at time t depends on price P(t) subject to
   Qs = c1 + w1 P(t) + u1 P'(t) + v1 P''(t) with c1, w1, u1, v1 constants,
and assume the quantity demanded Qd(t) at time t depends on price subject to
   Qd = c2 + w2 P(t) + u2 P'(t) + v2 P''(t) with c2, w2, u2, v2 constants.

Assuming price clears the market at each time t, Qs(t) = Qd(t), and we get a second order
ODE for P(t). Let  c = c1 - c2, w = w1 - w2,  u = u1 - u2, v = v1 - v2. Then we get, after
dividing by v # 0:

   P''(t) + (u/v) P'(t) + (w/v) P(t) = - c/v.

Since c = c1 - c2 and w = w1 - w2, where under ordinary supply conditions c1 < 0 & w1 > 0,
and under ordinary demand conditions c2 > 0 and w2 < 0, c is negative and w is positive,
so -c/w is positive and is the predicted possible eqilibrium price Pe found by setting
P''(t) = 0, P'(t) = 0, leaving (w/v) P(t) = - c/v, or P(t) = -c/w = positive constant = Pe.     

Here we ask ode2 to find the general solution. In answer to the question:
Is 4*v*w - u^2 positive, negative or zero, we first answer positive, or simply p.

kill(P,u,v,w,c)$
eqn : 'diff (P,t,2) + (u/v)*'diff (P,t) + (w/v)*P = -c/v;
gsoln : ode2 (eqn, P, t);

Likewise, you can get other solutions by answering negative or zero.

gsoln : ode2 (eqn, P, t);

gsoln : ode2 (eqn, P, t);

## Differential Equations with a Variable Term

We adapt the discussion in Chiang and Wainwright, Sec. 16.6.
The following examples show the use of ode2 to find a general solution.
These are all second order ODE's.

### Example 1

General solution of
y'' + 5 y' = 3 y = 6 t^2 - t -1.

eqn : 'diff (y,t,2) + 5*'diff (y,t) + 3*y = 6*t^2 - t - 1;
gsoln : ode2 (eqn, y, t);

### Example 2

General solution of y'' + 5 y' = 6 t^2 - t - 1.

eqn : 'diff (y,t,2) + 5*'diff (y,t) = 6*t^2 - t - 1;
gsoln : ode2 (eqn, y, t);

### Example 3

General solution of y'' + 3 y' - 4 y = 2 exp(- 4*t).

eqn : 'diff (y,t,2) + 3*'diff (y,t) - 4*y = 2*exp (- 4*t);
gsoln : ode2 (eqn, y, t);

## Higher Order ODE's With desolve (eqn, [y(t)])

ode2 only can only solve first or second order ODE's. If we try a 3rd order ODE, the 
eqn is returned and the return value is 'false'.

eqn : 'diff (y,t,3) + 2*'diff (y,t,2) + 'diff (y,t) + 2*y = 8;
gsoln : ode2 (eqn,y,t);

We can instead use desolve, which uses Laplace transform methods.
The syntax required in defining eqn is slightly different, since we must write the
dependent variable as y(t) instead of just y. The following illustrates desolve syntax.

eqn : 'diff (y(t),t,3) + 2*'diff (y(t),t,2) + 'diff (y(t),t) + 2*y(t) = 8;
gsoln : desolve (eqn, [y(t)]);

grind(gsoln)$

Use the atvalue function to set initial conditions and get a definite solution.

atvalue ('diff (y(t),t,2), t = 0, 1);
atvalue ('diff (y(t), t), t = 0, 2);
atvalue (y(t), t = 0, 4);
dsoln : desolve (eqn, [y(t)]);