Math 230, Spring 1999
Mathematica assignment 2
(due Monday, 2/8/99)

The following problems are taken from pages 59 to 67  of the Differential Equations with Mathematica book.  Please put some effort into presenting your work well.  That is, include
comments where appropriate, and  use the Format menus to distinguish text from headings
from actual computation.  It might work best just to do all the work and then go back and insert comments, format things right, etc.  Also, remember that when the instructions don't explicitly require the use of Mathematica, you're free to use pencil and paper to do some of the work--whatever gets the job done most effectively.

Problem 2.

a. (2 points)

soln[x_]  =  
      First[y[x] /. 
          DSolve[{y'[x]  +  y[x]/x == 2,  y[1] == c}, y[x], 
            x]]

b. (2 points)

Plot[
        Evaluate[Table[soln[x], {c,  .8, 1.2,  .1}]], {x,  .75, 
          1.25}];

c. (2 points)

Table[soln[10^n], {n, -2, 1}, {c,  .8, 1.2,  .1}]

% // TableForm

The columns correspond to c  =   .8,   .9,  1,  1.1, 
     1.2 going from left\(TextForm
    \`to\ right . \ \ The\ rows\ correspond\ to\ x\  = \ \(10\^n = \  .01\), 
    \  .1, \ 1, \ 10\ going\)from top to bottom .

d. (4 points)

Plot[
        Evaluate[Table[soln[x], {c,  .8, 1.2,  .1}]], {x, 0, 
          2.5}];

All solutions are asymptotic to the same line (looks like y=x)
as x goes to infinity--regardless of initial condition.  Solutions
with initial value greater than one tend to infinity as x goes
to zero.  Those with initial value less than one tend to negative 
infinity.  The solution with initial value one tends to zero as
x goes to zero.

\(\(






\)\)

Problem 4.

a. (2 points)

\(Clear[y]
\)
\(DSolve[\(y'\)[
          x]\  == \ \((x\  - \ Exp[\(-x\)])\)/\((y[x] + 
              Exp[y[x]])\), 
				\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ y[x], x]\)

b. (2 points)

I'll cut and paste the left hand side and give it a convenient
name  (it's important to remove the [x]'s from the y[x]'s):

lhs[x_, y_]  =  
      Simplify[ 
        (1/2) (E^(-x)) (-2 + 2 (E^(x + y)) - (E^x) (x^2) + 
              (E^x) (y^2))]

ContourPlot[lhs[x, y], {x, -1, 3}, {y, -2, 2}, 
        ContourShading -> False];

c. (3 points)

plot1  = 
        ContourPlot[lhs[x, y], {x, -1, 3}, {y, -2, 2}, 
          Contours  ->  {lhs[1.5,  .5]},  
          ContourShading -> False];

The top curve is the graph of the solution since it (evidently)
passes through the point (1.5,.5).

d. (3 points)

I'll use the previous picture to take a guess at the correct 
value of y and then use FindRoot to improve the guess:

FindRoot[lhs[0, y]  ==  lhs[1.5,  .5], {y,  .5}]
FindRoot[lhs[1, y]  ==  lhs[1.5,  .5], {y,  .5}]
FindRoot[lhs[1.8, y]  ==  lhs[1.5,  .5], {y,  .5}]
FindRoot[lhs[2.1, y]  ==  lhs[1.5,  .5], {y,  .5}]

That is:
	   y(0)     = .318386...
               y(1)     = .235633...
	   y(1.8) = .682187...
               y(2.1) = .866212...

\(\(Show[plot1, 
        Graphics[{PointSize[0.03], Point[{0,  .318386}], 
            Point[{1,  .235633}], 
		\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Point[{1.8,  .682187}], \ 
            Point[{2.1,  .866212}]}]];\)\)

Problem 10.

a. (4 points)

Before I plot the vector field, I'd like to check the graph of
the right side of the ODE as a function of y.  This will give me
some idea of what vertical range to use.  Since the equation is
autonomous, the horizontal range is unimportant.

Plot[3  Sin[y] + y - 2, {y, -10, 10}];

looks like I should expect equilibrium solutions at about y=1/2,
y=4 and y=5.  Therefore, I'll choose [-1,6] as a range for y.

Needs["Graphics`PlotField`"]

\(\(PlotVectorField[{1, 3  Sin[y]\  + \ y\  - \ 2}, {t, 0, 
          3}, {y, \(-1\), 6}, Axes -> \ True, 
	
        ScaleFunction -> \((1 &)\), Ticks -> None, 
        Frame -> True, AspectRatio -> 1, \ 
        PlotPoints\  -> \ 15];\)\)

To find the equilibrium solutions:

FindRoot[3  Sin[y] + y - 2 == 0, {y,  .5}]
FindRoot[3  Sin[y] + y - 2 == 0, {y, 4}]
\(FindRoot[3  Sin[y] + y - 2 == 0, {y, 5}]
\)

From this and the two plots above, it is apparent that the ODE has
two unstable equilibrium solutions:  y= 4.9295 and y=.517049
and one stable equilibrium solution: y= 3.77452.

b. (4 points)

This ODE is not autonomous, so I cannot proceed as I did in the last
example.  I'll have to get along with just the direction field.  
Note that the right side vanishes (i.e. the direction field is
horizontal) when y=0 and y=x.

\(\(PlotVectorField[{1, y^2\  - x*y}, {x, \(-4\), 4}, {y, \(-4\), 4}, 
        Axes -> \ True, 
	ScaleFunction -> \((1 &)\), 
        Ticks -> None, Frame -> True, AspectRatio -> 1, \ 
        PlotPoints\  -> \ 15];\)\)

Looks like y=0 is a constant solution of the ODE--an observation
that is easily verified by plugging into the ODE.  All solutions
y[x] starting below y= 0 appear to converge to y=0 as x goes to
infinity.  Solutions y[x] starting above  y= 0 appear to 
converge to either infinity or 0.  Solutions exhibiting these two
types of behavior appear to be separated by a curve asymptotic to
the line y=x as x goes to infinity and y=0 as x goes to negative
infinity.   At any rate it is clear from the right side of the 
differential equation that the dividing curve must lie at or 
above the line y=x, because the right side is negative when 
y < x (and y > 0) so that solutions lying between y=0 and y=x will
slope downward as x becomes large (and positive).

c. (4 points)

DSolve[{y'[x] ==  y[x]^2  -  x *y[x], y[0]  ==  c}, y[x], 
      x]

What the *&^&* is this Erf thing?

Plot[Erf[x/Sqrt[2]], {x, -10, 10}];

Doesn' t really matter,  
      I guess .   The point is that Erf[x/(Sqrt[2])]
increases from  - 
            1 to 1 as x goes from negative infinity to 
infinity .  Therefore,  
          if c = y[0] is less than (Sqrt[2/π]) then it' s clear 
\(that\ y[x]\ is\ defined\ for\ all\ x\ \((that\ is, \ 
          the\ denominator\ never
vanishes)\) 
        and\ that\ y[x]\ converges\ to\ zero\ as\ x\ goes\ to\ infinity\)
\(\((because\ the\ numerator\ of\ y[
                  x]\ definitely\ goes\ to\ zero\ while\ the

                denominator\ doesn' t)\) . \ \ If\ y[0]\  > \@\(2/π\), \ 
      then\ at\ some\ finite\)
\(x\ the\ denominator\ vanishes\ and\ the\ solution\ goes\ to\ infinity\ along\)
\(a\ vertical\ asymptote\ \((this\ is\ different\ than\ what\ I\ predicted\ from\ the\ 
direction\ field)\) . \ \ It' 
        s\ harder\ to\ tell\ just\ what\ happens\ if\)
y[0]  =  (Sqrt[2/π]),  
      so I plot that particular solution:

Plot[-((2 (Sqrt[2/π]) (E^(-((x^2)/2))))/(-2 + (Sqrt[2/π]) (Sqrt[2 π]) Erf[(x/(Sqrt[2]))])), {x, 0, 7}];

This looks pretty close to a straight line, so I compare with  the line y=x that looked important in the direction field:

Plot[{-((2 (Sqrt[2/π]) (E^(-((x^2)/2))))/(-2 + (Sqrt[2/π]) (Sqrt[2 π]) Erf[(x/(Sqrt[2]))])), x}, {x, 0, 7}];

Sure enough, the "dividing" solution is asymptotic to the line y=x.

Problem 12.

a. (3 points)

Needs["Graphics`PlotField`"]

\(PlotVectorField[{1, \(-t\)*x^2}, {t, \(-3\), 3}, {x, \(-3\), 3}, 
      Axes -> \ True, 
	ScaleFunction -> \((1 &)\), 
      Ticks -> None, Frame -> True, AspectRatio -> 1, \ 
      PlotPoints\  -> \ 19]\)

Looks like x = 0 is a constant solution.Moreover, it looks like 
solutions x[t] that start above x=0 will converge to zero as t goes
to either positive or negative infinity,whereas solutions x[t] 
that start below x=0 will converge to negative infinity.

b. (3 points)

I'm going to get the solution to the ODE in a slightly 
different form than is asked for in the book.

soln  =  
      First[x[t] /. 
          DSolve[{x'[t] == -t*x[t]^2, x[0]  ==  c}, x[t], 
            t]]

Plot[
        Evaluate[Table[soln, {c, -3, 3, 1}]], {t, -5, 
          5}];

So in fact, my predictions from the direction field weren't quite
right.  It's true that any solution satisfying x[0] >0 will 
converge  to 0 as t goes to positive or negative infinity.  
However, any negative solution appears to have to vertical asymptotes and
to go to negative infinity along both of them--i.e. solutions that start below 
x-axis go to negative infinity at some *finite* positive and negative 
values of t.

c. (4 points)

Five types of solution curves:
1. x[t] < 0, defined for all t in some finite interval symmetric about the 
    y-axis, and x[t] tends to negative infinity at each end of the interval.
2. x[t] = 0 for all t
3. x[t] > 0, defined for all t and tending to zero as t goes to ±∞.
4. x[t] > 0 everywhere decreasing, and defined for t in an interval of the form [a,∞)
      where a > 0; asymptote at t = a.
5. x[t] > 0 everywhere, increasing and defined for t in an interval of the form (-∞,-a]
      where a > 0; asymptote at t = -a.

d. (2 points)

This `Show' command is a little tricky--took some playing around to get a good
picture.

\(\(\(\ \ \)\(plot1\  = \ 
        PlotVectorField[{1, \(-t\)\ *x^2}, {t, \(-3\), 3}, {x, \(-4\), 4}, \ 
          PlotPoints\  -> \ 17, 
	\ \ \ ScaleFunction -> \((1 &)\), 
          Axes -> \ True, Ticks -> None, Frame -> True, 
          AspectRatio -> 1];\)\)\)

\(\(\(	\ \)\(plot2\  = \ 
        Plot[Evaluate[Table[soln, {c, \(-3\), 3, 1}]], {t, \(-3\), 
            3}, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ PlotRange -> {\(-4\), 4}, 
          AspectRatio -> 1];\)\)\)

Show[plot1, plot2, PlotRange -> {-4, 4}, Frame -> False, 
        Ticks -> None];

Problem 16.

a. (3 points)

\(pop[t_, a_, b_, c_]\  = \ 
      Expand[First[
          p[t] /. DSolve[
				\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {\(p'\)[t]\  == \ a\ p[t] + b, \ 
                p[0]\  == \ 
                  c}, 
			\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ p[
                t], t]]]\)

The parameter c is simply the population of Nevada in 1950--i.e. the initial
value.  The parameter a determines the contribution to population growth
coming from births - deaths , and b is the contribution to growth 
coming from (people moving in) - (people moving out).

b (3 points)

\(Solve[{pop[0, a, b, c] == 285.3, 
        pop[10, a, b, c] == 
          488.7, 
		\ \ \ \ \ \ \ \ \ \ \ \ pop[20, a, b, c] == 
          800.5}, {a, b, c}]\)

Only the first in this list of solutions gives real values for a,b,
and c.  I'll use that solution.

c. (4 points)

population[t_]  =  pop[t, 0.042718,  4.1162, 285.3]

Now Icheck the population in 1950 and 1990 predicted by 
this model:

population[-10]

population[30]

These are at least close to the true values.  To predict when the population reaches three
million, I'll guess that it'll be fifteen years after 1990 (since the population nearly
tripled in the twenty years between 1970 and 1990).

FindRoot[population[t] == 3000, {t, 45}]

So I'm predicting a population of three million by 2009
Given the real data for 1990 and 1950 this model seems to be 
projecting too much population growth.  I'd probably guess that
it'll be a few years after 2009 when the population really
does reach three million.  Based on the formula above, it looks like I'd
have to decrease both a and b to get a more accurate population model.

Plot[population[t], {t, -10, 90}];

If the model is right, Nevada will get pretty crowded in time!

Problem 17.

a.  (2 points)

Mathematica is really overkill for this part...

Solve[(α - 1) x -  x^3  ==  0, x]

That is, we get three real roots when α>1, and only one real
root when α<1.  When α=1, we also get one real root, but with
multiplicity three.

b. (2 points)

\(\(Table[
        PlotVectorField[
		{1, \((α - 1)\) x - \ x^3}, {t, \(-2\), 
            2}, {x, \(-2\), 2}, Axes -> \ True, 
	
          ScaleFunction -> \((1 &)\), Ticks -> None, 
          Frame -> True, AspectRatio -> 1, \ 
          PlotPoints\  -> \ 15], {α, \(-2\), 0}];\)\)

The picture is more or less the same in all three cases:  x=0 is
a stable equilibrium solution.

c.  (2 points)

\(PlotVectorField[
		{1, \(-x^3\)}, {t, \(-2\), 2}, {x, \(-2\), 2}, 
      Axes -> \ True, 
	ScaleFunction -> \((1 &)\), 
      Ticks -> None, Frame -> True, AspectRatio -> 1, \ 
      PlotPoints\  -> \ 15]\)

x=0 is  still a (barely) stable equilibrium.

d. (3 points)

\(PlotVectorField[
		{1,  .5  x\  - \ x^3}, {t, \(-2\), 
        2}, {x, \(-2\), 2}, Axes -> \ True, 
	
      ScaleFunction -> \((1 &)\), Ticks -> None, 
      Frame -> True, AspectRatio -> 1, \ 
      PlotPoints\  -> \ 17]\)

Now x = 0 is an unstable equilibrium,  
    but there are two stable equilibrium\(TextForm\`solutions\ at\ x\  = \ \(±\@ .5\)\)

\(PlotVectorField[
		{1, x\  - \ x^3}, {t, \(-2\), 2}, {x, \(-2\), 2}, 
      Axes -> \ True, 
	ScaleFunction -> \((1 &)\), 
      Ticks -> None, Frame -> True, AspectRatio -> 1, \ 
      PlotPoints\  -> \ 17]\)

This is similar to the last case, except that the stable 
equilibria have migrated to x=±1.

e.  ( 3 points)

When α<0, we had a single stable equilibrium solution at x=0.
When α increased past 0, we suddenly had two stable equilibrium
that started just above and below x=0 and moved further from
x=0 as α increased.