7.1 Modeling with Differential Equations

Mathematica script by Chris Parrish,
   cparrish@sewanee.edu

Sources and references for some of these problems include
    James Stewart, "Calculus: Concepts and Contexts," Second Edition, Brooks/Cole, 2001
    Deborah Hughes-Hallett, Andrew M. Gleason, et. al., "Calculus," Second Edition, John Wiley & Sons, 1998
    Robert Fraga, ed., "Calculus Problems for a New Century," The Mathematical Association of America, 1993

How a Person Learns

Hughes-Hallett, Gleason, et al, Text Example, page 478

In[138]:=

<<Graphics`Colors`

Let dy/dt = f[t,y]. Here is the data we are working with.

In[139]:=

Clear[f,t,y,k]

f[t_,y_] := 100 - y;  (* dy/dt = f[t,y]; t in weeks, y is % of task learned *)

t0 = 0.0;
y0 = 0.0;

deltaT = 0.2;      (* 1 day = 1/5th week *)

m0 = f[t0,y0];     (* m0 = slope of solution curve at t = t0 *)

We need the % of the task that has been learned by the end of the first day.

In[145]:=

y1 = y0 + f[t0,y0] deltaT

Out[145]=

20.

By the end of the second day, ...

In[146]:=

t1 = t0 + deltaT;

y2 = y1 + f[t1,y1] deltaT

Out[147]=

36.

There is a pattern here. On the k+1 st day, ...

In[148]:=

t[0]   = t0;
t[k+1] = t0 + k deltaT;

y[0]   = y0;
y[k+1] = y[k] + f[t[k],y[k]] deltaT;

We can ask Mathematica to do that sort of incremental evaluation,
and accumulate the series of numbers in a variable named "result."

In[152]:=

y = y0; t = a = t0;
result = {{t,y}};
Do[y = y + f[t,y] deltaT;
   t = a + k deltaT;
   AppendTo[result,{5 t,y}],  (* t is in weeks, but we want the result to be in days *)
   {k,20}];

TableForm[result,
          TableHeadings -> {{"% of Task Learned"},{"time (days)","%"}},
          TableSpacing -> {0,3}]

Out[155]//TableForm=

time (days) %
% of Task Learned 0. 0.
1. 20.
2. 36.
3. 48.8
4. 59.04
5. 67.232
6. 73.7856
7. 79.0285
8. 83.2228
9. 86.5782
10. 89.2626
11. 91.4101
12. 93.1281
13. 94.5024
14. 95.602
15. 96.4816
16. 97.1853
17. 97.7482
18. 98.1986
19. 98.5588
20. 98.8471

Let's put this data in a picture.

In[156]:=

ListPlot[result,
         PlotStyle -> {CornflowerBlue,PointSize[0.02]},
         PlotLabel -> "% of Task Learned",
         AxesLabel -> {"time (days)","%"}];

[Graphics:HTMLFiles/7.1_modelingWithDEs_4.gif]

Logistic Equation

Hughes-Hallett, Gleason, et al, Exercise 9.1.6, page 482

In[157]:=

Clear[p, t] <br /> p[t_] := 1/(1 + ^(-t)) ; <br /> Plot[p[t], {t, -5, 5}, Plot ... lotLabel"Logistic Equation", AxesLabel {"t", None}] ;

[Graphics:HTMLFiles/7.1_modelingWithDEs_6.gif]

In[160]:=

p '[t]

Out[160]=

^(-t)/(1 + ^(-t))^2

In[161]:=

p[t] ( 1 - p[t]) //Simplify

Out[161]=

^t/(1 + ^t)^2

They don't look the same. Are they?
Try harder, Mathematica.

In[162]:=

p '[t] - p[t] (1 - p[t])//Simplify

Out[162]=

0

Yup. They are same -- so p[t] is a solution of the Logistic Equation.

Hanging Cable

Hughes-Hallett, Gleason, et al, Exercise 9.1.8a, page 482

In[163]:=

Clear[y, x]  y[x_] := (^x + ^(-x))/2 Plot[y[x], {x, -5, 5}, &# ... 71;PlotLabel"Hanging Cable", AxesLabel {"x", None}] ;

[Graphics:HTMLFiles/7.1_modelingWithDEs_13.gif]

Now calculate.

In[166]:=

y''[x] - (1 + y '[x]^2)^(1/2)

Out[166]=

1/2 (^(-x) + ^x) - (1 + 1/4 (-^(-x) + ^x)^2)^(1/2)

They don't look the same. Are they?
Try harder, Mathematica.

In[167]:=

y''[x] - (1 + y '[x]^2)^(1/2)//Simplify

Out[167]=

1/2 (^(-x) + ^x - (2 + ^(-2 x) + ^(2 x))^(1/2))

This is about as far as Mathematica can take it. Can you take it farther?

Hughes-Hallett, Gleason, et al, Exercise 9.1.8b, page 482

In[168]:=

Clear[y, x, k, a] <br /> y[x_] := (^(a x) + ^(-a x))/(2a) ;

Let's calculate the two sides of the equation separately.

In[170]:=

y''[x]

Out[170]=

(a^2 ^(-a x) + a^2 ^(a x))/(2 a)

In[171]:=

(1 + y '[x]^2)^(1/2)//Simplify

Out[171]=

1/2 (2 + ^(-2 a x) + ^(2 a x))^(1/2)

If y[x] is a solution of the more general differential equation,
the quotient of these two results should give us k.

In[172]:=

y''[x]/(1 + y '[x]^2)^(1/2)//Simplify

Out[172]=

(a ^(-a x) (1 + ^(2 a x)))/(2 + ^(-2 a x) + ^(2 a x))^(1/2)

If you are given k, what's a?

Families of Curves

Hughes-Hallett, Gleason, et al, Exercise 9.1.8a, page 482

In[173]:=

Clear[y, x, k]  y[x_] := x ^(k x) ;

In[175]:=

y'[x]

Out[175]=

^(k x) + ^(k x) k x

In[176]:=

Table[{y[x]/x,
       x Log[y[x]]/x,
       y[x]/x (1 + Log[y[x]/x]),
       (y[x] Log[y[x]])/(x Log[x])}]

Out[176]=

{^(k x), Log[^(k x) x], ^(k x) (1 + Log[^(k x)]), (^(k x) Log[^(k x) x])/Log[x]}

Will the REAL dy/dx please come forward!

In[177]:=

Clear[y, x, k, m]     (* Our families ! *)<br /> y[x_] := x ^(k x) ; y[x_] := x^p ; y[x_] := ^(k x) ; y[x_] := m x ;


Created by Mathematica  (April 22, 2004)