Slope Fields

Enter the following line to load the necessary packages.

<< Graphics`Graphics`
<< Graphics`PlotField`

Part 1. The Slope Field Concept

Enter your answers to questions in Part 1 in the cell below:

Answer:

Define the function f as the right-hand side of the differential equation (dP/dt) = t - P.

f[t_,  P_] :=  t - P;

The next command draws the slope (direction) field.  Mathematica does not have a direction field command, but we can use PlotVectorField.  A direction field is just a special case of a vector field.

field  =  
        PlotVectorField[{1, f[t, P]}, {t, -6, 6}, {P, -6, 6}, 
          ScaleFunction -> (0.5 &), HeadWidth -> 0,  
          HeadLength -> 0, Axes -> True, PlotPoints -> 20, 
          AxesLabel -> {"t", "P"}];

We add solution curves by solving the differential equation under the specified initial conditions and then plotting them with the direction field.  For efficiency, we will use Mathematica's NDSolve to find a numerical solution instead of an exact analytical solution.  Note that we could not write P'[t] == f[t, P], since Mathematica needs to be explictly told that P is a function of t.  Here is the same slope field with three solution curves.

initial = {P[0] == 0, P[0] == -1, 
          P[0] == -6};
solutions  =  
        Table[NDSolve[{P'[t] == t - P[t], initial[[k]]},  
            P,  {t,  -6,  6}], {k, Length[initial]}];
solutionsPlot  =  
        Table[Plot[
            Evaluate[P[t] /. solutions[[i]]],  {t,  -6,  6},  
            PlotRange -> {-6, 6},  
            PlotStyle -> RGBColor[1, 0, 0], 
            DisplayFunction -> Identity], {i,  
            Length[solutions]}];
Show[field,  solutionsPlot,  PlotRange -> {-6, 6}, 
        DisplayFunction -> $DisplayFunction];

Part 2. Slope Fields for Natural Growth

Enter the value of k.

k  =  0.1;

Construct the slope field, and assign it a name for later use:

field  =  
        PlotVectorField[{1, k*P}, {t, 0, 20}, {P, 0, 1000}, 
          ScaleFunction -> (100 &),  HeadWidth -> 0, 
          Axes -> True, AspectRatio -> 1, 
          AxesLabel -> {"t", "P"}];

Answer to question 2:

Redraw the slope field with three solution curves.  Change value1, value2, and value3 in the list initial to numbers of your choice.

initial = {P[0] ==  value1, P[0] ==  value2, 
          P[0] ==  value3};
intcurves  =  
        Table[NDSolve[{P'[t] == k*P[t], initial[[i]]},  
            P,  {t,  0, 20}], {i, Length[initial]}];
intcurvePlot  =  
        Table[Plot[Evaluate[P[t] /. intcurves[[i]]],  {t,  0,  20}, 
            PlotStyle -> RGBColor[0, 0, 1],  
            DisplayFunction -> Identity], {i,  3 }];
Show[field,  intcurvePlot,  
        DisplayFunction -> $DisplayFunction,  
        PlotRange -> {0, 1000}];

Enter your formulas for three solution curves.  Name the starting populations P1, P2, P3, and the plots of these curves solution1, solution2, solution3, respectively.

P1[t_] :=  ?;
solution1  =  
        Plot[P1[t], {t, 0, 20},  PlotStyle -> RGBColor[0.5, 0.5, 0], 
          DisplayFunction -> Identity];
P2[t_] :=  ?;
solution1  =  
        Plot[P2[t], {t, 0, 20},  PlotStyle -> RGBColor[0.2, 1, 0.8], 
          DisplayFunction -> Identity];
P3[t_] :=  ?;
solution1  =  
        Plot[P3[t], {t, 0, 20},  PlotStyle -> RGBColor[0, 0, 1], 
          DisplayFunction -> Identity];

Display the field and solution plots together:

Show[field, solution1, solution2, solution3,  
        PlotRange -> {0, 1000}, 
        DisplayFunction -> $DisplayFunction];

Part 3. Slope Fields for Radioactive Decay

The Meyer-von Schweidler data and a plot of the data points:

activity  =  {{0.2, 35.0}, {2.2, 25.0}, {4.0, 22.1}, {5, 17.9}, {6, 
            16.8}, {8, 13.7}, {11, 12.4}, {12, 10.3}, {15, 7.5}, {18, 
            4.9}, {26, 4.0}, {33, 2.4}, {39, 1.4}, {45, 
            1.1}};
actplot  =  
        ListPlot[activity, 
          PlotStyle -> {RGBColor[0, 0, 1], 
              PointSize[0.015]}];

Enter (Subscript[P, 0]) and b:

P0  = 28;
b  =  0.92466;

Find k such that P = (Subscript[P, 0])(b^t) = (Subscript[P, 0])(e^(k t)).

k  =  ?;

Answer to the question in step 2:

Draw the slope field.  Then add the solution starting at (Subscript[P, 0]).

field  =  
        PlotVectorField[{1, k*P}, {t, 0, 45}, {P, 0, 35}, 
          ScaleFunction -> (0.5 &), HeadWidth -> 0,  
          HeadLength -> 0,  Axes -> True, 
          AxesLabel -> {"t", "P"}];

solution  =  
        NDSolve[{P'[t] == k*P[t],  P[0] ==  P0},  
          P,  {t,  0, 45}];
solutionPlot  =  
        Plot[Evaluate[P[t] /. solution],  {t, 0, 45},  
          PlotStyle -> RGBColor[1, 0, 0], 
          DisplayFunction -> Identity];
Show[field,  solutionPlot,  
        DisplayFunction -> $DisplayFunction];

Display the data plot along with the field plot and solution.

Show[actplot, field,  solutionPlot];

Answer to the question in step 3:

Part 4. Slope Fields for Limited Population Growth

Enter the constants and function definition for the limited growth model.  Fill in the right-hand side of dP/dt at the "?" below.  Remember that this should be done in terms of P and not P[t].

M  =  1000;
c  =  0.000098;
rhs  =  ?;

Construct the direction field:

field  =  
        PlotVectorField[{1, rhs}, {t, 0, 100}, {P, 0, 1500}, 
          ScaleFunction -> (200 &),  Axes -> True, 
          AspectRatio -> 1, 
          AxesLabel -> {"t", "P"}];

Answer to the question in step 2:

Answer to the question in step 3:

We repeat here the commands from the Limited Population Growth module for constructing an approximate solution, using repeated applications of "rise = slope * run." You may vary n if you wish.

n  = 50;
Delta  = 100/n;
Clear[t, p, slope];
\(\(Do[
	
        t[k]\  = \ k*Delta, 
	{k, 0, 
          n}
];\)\)
p[0] = 111;
slope[0]  =  c*p[0]*(M - p[0]);
\(\(Do[
	
        p[k]\  = \ p[k - 1] + slope[k - 1]*Delta; \ 
	
        slope[k]\  = \ c*p[k]*\((M - p[k])\), 
	{k, 1, 
          n}
];\)\)

popData  =  Table[{t[k], p[k]}, {k, 0, n}];
popGraph  =  
        ListPlot[popData, 
          PlotStyle -> {RGBColor[0, 0, 1], 
              PointSize[0.015]}];

Now, we display the approximate solution on the direction field.

Show[popGraph, field];

Redraw the field with an "automatic" solution, and draw again with the approximate solution.

solution  =  
        NDSolve[{P'[t] == c*P[t]*(M - P[t]),  P[0] == 111},  
          P,  {t, 0, 100}];
solutionPlot  =  
        Plot[Evaluate[P[t] /. solution],  {t,  0,  100},  
          DisplayFunction -> Identity,  
          PlotStyle -> RGBColor[1, 0, 0]];
Show[popGraph, field,  solutionPlot,  
        DisplayFunction -> $DisplayFunction];

Answer to question 6:

Part 5. Summary