Part III. First Order Ordinary Differential Equations
Section 2. Working with Solutions: Modeling
A mathematical model often includes differential equations describing how one variable changes in relation to  another. Examples in this section show how Mathematica can be used to help in the analysis of differential equations  arising in the modeling process.

Exponential models

When a bank compounds interest continuously the principal,  P, as a function of time t (in years) grows  exponentially. The governing differential equation, given in terms of an annual interest rate  r, is

P' = r P

Question: How long will it take $1000 to grow to $1500 at an annual interest rate of 4%?

In[1]:=

DE = P'[t] == 0.04*P[t]
soln = DSolve[ {DE, P[0]==1000}, P, t ]

Out[1]=

P^′[t] 0.04 P[t]

Out[2]=

{{PFunction[{t}, 1000 ^(0.04 t)]}}

Make the solution into the function named g and solve g(t) = 15000.

In[3]:=

g := soln[[1,1,2]]
timesoln = FindRoot[ g[t]==1500, {t,1} ]

Out[4]=

{t10.1366}

It will take about 10.14 years. Here is plot of the principal and the point when P(t) = 1500.

In[11]:=

Show[  ListPlot[ {{10.14,1500}}, PlotStyle->PointSize[0.03],
                  DisplayFunction->Identity ],
       Plot[ g[t], {t,0,12}, DisplayFunction->Identity ],
       PlotRange->{0,1600}, DisplayFunction->$DisplayFunction]

[Graphics:HTMLFiles/Math_P3S2_4.gif]

The exponent on e in the solution formula tells us that the principal will triple in a little over 25 years (the  multiplier for t = 25 would be e, which is 2.78). Recall that the doubling time is ln(2)/r. See Ledder, Chapter 1. See below.

In[12]:=

Log[2]/0.04*years

Out[12]=

17.3287 years

A more interesting problem

Juanita's credit card bill is $2500, she pays an annual 9% interest rate, compounded continuously. She decides  to pay it off by transferring money from her savings account continuously at the rate of $100 per month.  Assuming she has $35,000 in savings, which accumulates interest at the annual rate of 2%, also compounded  continuously, how long will it take? What will her total payments be? How much money will she have in her  savings account when the debt is settled?

Let y(t) be the amount that Juanita owes the creditor at time  t (in years), y(0) = 2500. Since she reduces debt at  a constant yearly rate of $1200, the differential equation for y(t) looks like this.

In[13]:=

DE2 = y'[t] == 0.09*y[t] - 1200
soln = DSolve[ {DE2, y[0]==2500}, y, t]

Out[13]=

y^′[t]  -1200 + 0.09 y[t]

Out[14]=

{{yFunction[{t}, 2.71828^(-0.09 t) (-10833.3 7.38906^(0.09 t) + 13333.3 ^(0.09 t))]}}

Name the solution function g once more. The debt is liquidated when g(t) = 0.

In[16]:=

g := soln[[1,1,2]]
debtfree = FindRoot[ g[t]==0, {t,1}]

Out[17]=

{t2.3071}

Check the solution.

In[18]:=

g[t]/.debtfree

Out[18]=

0.

It will take Juanita about 2.3 years to pay off the debt. This is between 27 and 28 months.

In[19]:=

(12*t/.debtfree)*months

Out[19]=

27.6852 months

In[25]:=

Plot[ g[t], {t,0,2.5}, AspectRatio->1/3 ]

[Graphics:HTMLFiles/Math_P3S2_12.gif]

Her total payment can be calculated as the product of the annual payment rate, 1200, and the time it took to  settle the debt.

In[26]:=

(1200*t/.debtfree)*dollars

Out[26]=

2768.52 dollars

To calculate the balance in Juanita's savings account, let  x(t) denote the principal at time  t, x(0) = 35000. The  differential equation for x(t), and its solution are entered below.

In[31]:=

DE3 = x'[t] == 0.02*x[t] - 1200
savings = DSolve[ {DE3, x[0]==35000}, x, t ]

Out[31]=

x^′[t]  -1200 + 0.02 x[t]

Out[32]=

{{xFunction[{t}, 2.71828^(-0.02 t) (-25000. 7.38906^(0.02 t) + 60000. ^(0.02 t))]}}

The savings balance at time of payoff can be calculated as follows.

In[34]:=

x[t]/.savings/.debtfree

Out[34]=

{33819.4}

Monthly payments

It is unlikely that any bank would agree to make continuous payouts for Juanita. On the other hand, many  banks offer the service of regularly monthly payments. Let's compare the total payout calculated above with  the payout if Juanita has her bank send $100 to arrive to her creditor on day 30 of month 1 and $100 on day 30 of every month thereafter. Assume that each month has 30 days and a year has 360.
Let  d[n] be Juanita's credit card debt on the first day of month n. Then  d[1] = 2500 and d[2], her debt on the first day of month 2, is given  by the following formula:

d[2] = 2500 e^0.09/12- 100

Remember that time is in years so one month is 1/12 years.
The function d can be defined recursively in Mathematica as follows.

In[48]:=

d[1] = 2500;
d[n_] := d[n-1]*Exp[0.09/12] - 100

At the beginning of the 27th month, Juanita has paid $2600 and owes d[27] dollars, $178.22. See below.

In[47]:=

d[27]

Out[47]=

178.216

One more month goes by, the 27th payment of $100 is made, and Juanita still owes the following amount.

In[50]:=

d[28]

Out[50]=

79.5581

Assuming she pays it all off at that time, her total payout is $2779.56. This is about $10 more than it would be if the payouts were made continuously. See the exercises for this section.  

Heating and cooling

Isaac Newton proposed a simple model for the temperature T of an object that cools or warms while immersed  in a medium of ambient temperature A (in the absence of convective effects). If  T(t) denotes the object's  temperature at time t and A(t) is the ambient temperature at the same time, then Newton's law of cooling is the  assumption that

T' = k ( A(t) - T(t) )

where k is a positive constant. An additional heating or cooling term can be added by changing the equation to

T' = k ( A(t) - T(t) ) + H(t)

If this equation is given to DSolve, Mathematica outputs the solution as a textbook integral. (See Ledder's Variation of  Parameters formula, Chapter 4, page 260.)  The following input asks for the solution in form T[t] -> ... .

In[57]:=

T'[t] == k*(A[t] - T[t]) + H[t]
DSolve[ %, T[t], t ]

Out[57]=

T^′[t] H[t] + k (A[t] - T[t])

Out[58]=

{{T[t] ^(-k t) C[1] + ^(-k t) ∫_1^t^(k K$535) (k A[K$535] + H[K$535]) K$535}}

The strange looking integral. Mathematica has made the solution formula in terms of a definite integral and has chosen the symbol K$535 for the variable of integration. Send an e-mail to Mathematica and ask why.

Constant ambient temperature, no heat

If the ambient temperature is constant and there is no additional heating or cooling, the temperature T(t) will  move towards the ambient temperature A exponentially. The number 1/k, called the "time constant", tells us  how fast T(t) approaches A.

In[59]:=

T'[t] == k*(A - T[t])
DSolve[ %, T[t], t ]

Out[59]=

T^′[t] k (A - T[t])

Out[60]=

{{T[t] A + ^(-k t) C[1]}}

Observe that for every increase in t by 1/k the difference between T(t) and A will decrease by a factor of 1/e  (about 1/3).This is useful for graphing purposes and making rough estimates. In a typical problem one knows  the initial temperature T(0), the ambient temperature A, and the object's temperature at some later time, T(t0).
Example. Joe is served a cup of coffee at 10 AM. Its temperature is 150 degrees F. Five minutes later the coffee  temperature has dropped to 110 degrees. If the room temperature is 66 degrees, when will the coffee  temperature be 100 degrees?

Start with the differential equation and the solution to the IVP obtained using the initial condition T(0) = 150.  (10 AM is taken as time 0.)  

In[72]:=

DEcoffee = T'[t] == k*(66 - T[t])
soln = DSolve[ {DEcoffee, T[0]==150}, T, t ]

Out[72]=

T^′[t] k (66 - T[t])

Out[73]=

{{TFunction[{t}, 6 ^(-k t) (14 + 11 ^(k t))]}}

Name the solution function g. The temperature at t = 5 is 110. This is used to find k.

In[74]:=

g := soln[[1,1,2]]
ksoln = Solve[ g[5] == 110, k ]

Solve :: ifun : Inverse functions are being used by Solve, so some solutions may not be found; use Reduce for complete solution information.  More…

Out[75]=

{{k1/5 Log[21/11]}}

The time constant is calculated first. The natural entry shown below seems to fail, N is not applied.

In[76]:=

N[1/k]/.ksoln

Out[76]=

{5/Log[21/11]}

The next entry illustrates an alternative way to apply N. Follow the calculation of 1/k with //N.

In[77]:=

1/k/.ksoln//N

Out[77]=

{7.73243}

This tells us that if the coffee sits for 7.7 minutes, the difference (84 degrees) between g(0) and the ambient  temperature will decrease to about 1/3 of its value. That is, the difference will be about 30 degrees:  g(7.7) is  roughly 96. We will find out when g(t) = 100 using FindRoot but first let's give k the value described by ksoln.

In[81]:=

k = k/.ksoln[[1,1]]

Out[81]=

1/5 Log[21/11]

You might like to see the formula for g[t] now that k has its exact value.

In[82]:=

g[t]

Out[82]=

2 (11/7)^(t/5) 3^(1 - t/5) (14 + 11^(1 - t/5) 21^(t/5))

When is the coffee at 100 degrees? At about 7 minutes after 10.

In[83]:=

FindRoot[ g[t]==100, {t,5} ]

Out[83]=

{t6.99365}

The next plot shows the graph of the temperature, g(t), the ambient temperature, 66, and dashed lines indicating when the coffee  reaches 100 degrees.

In[96]:=

Show[ ListPlot[ {{0,100},{6.99,100},{6.99,0}}, PlotJoined->True,
                 PlotStyle->Dashing[{0.02,0.02}]],
      Plot[ {g[t],66}, {t,0,15}], PlotRange->{0,160} ]

[Graphics:HTMLFiles/Math_P3S2_33.gif]

Variable ambient temperature, no heat

The next example considers the case where the ambient temperature changes linearly. The solution formula  shows that the temperature of the immersed object also approaches a straight line (same slope).
Example. Overnight the temperature in the house dropped to 45 degrees. At sunrise (7 AM) the outdoor  temperature was also 45 degrees and increased linearly to 75 degrees at 12 noon. What is the temperature  inside the house at 12 noon if it had increased to 50 degrees at 9 AM? The furnace in the house is turned off.
Take 7 AM to be time t = 0. The ambient temperature is given by the function  

A(t) = 45 + 6 t

t measured in hours. The temperature T(t) inside the house at time t satisfies T(0) = 45 and T(2) = 50. Here is  the differential equation and the solution to the IVP,  T(0) = 45. We first clear the k variable.

In[151]:=

Clear[k]
DEhouse = T'[t] == k*(45 + 6*t - T[t])
soln = DSolve[ {DEhouse, T[0]==45}, T, t ]

Out[152]=

T^′[t] k (45 + 6 t - T[t])

Out[153]=

{{TFunction[{t}, (3 ^(-k t) (2 - 2 ^(k t) + 15 ^(k t) k + 2 ^(k t) k t))/k]}}

Name the solution function h and find the value of k as we did above.

In[154]:=

h := soln[[1,1,2]]
ksoln = Solve[ h[2]==50, k ]

InverseFunction :: ifun : Inverse functions are being used. Values may be lost for multivalued inverses.  More…

Solve :: ifun : Inverse functions are being used by Solve, so some solutions may not be found; use Reduce for complete solution information.  More…

Out[155]=

{{k1/14 (12 + 7 ProductLog[-12/(7 ^(12/7))])}}

Let's assign k the value found above, but in floating point form, and check that h(2) = 50.

In[156]:=

k = N[ksoln[[1,1,2]] ]
h[2]

Out[156]=

0.597879

Out[157]=

50.

It would also be wise to verify that there is only one k value that works (see the warning above). You are asked to do that in the  exercises.
The next plot shows the ambient temperature (dashed) and the house temperature from 7 AM to 12 noon.

In[158]:=

Plot[ {45+6*t,h[t]}, {t,0,5}, PlotRange->{40,75},
      PlotStyle->{Dashing[{0.02,0.02}],Dashing[{0,0}]}]

[Graphics:HTMLFiles/Math_P3S2_41.gif]

Out[158]=

⁃Graphics⁃

The temperature inside the house at noon is 65.5 degrees.

In[113]:=

h[5]

Out[113]=

65.4695

Variable ambient temperature, variable heat

Finally, suppose the same house, on the same day, has a furnace that cycles on and off providing just enough heat to cause the temperature to increase at the rate of 6 degrees per hour at its peak output. It cycles on and  off sinusoidally over a 2 hour period:

H(t) = 3 - 3 cos((2π t)/2)

In[116]:=

H[t_] := 3 - 3*Cos[2*Pi*t/2]
Plot[ H[t], {t,0,5}, AspectRatio->1/3 ]

[Graphics:HTMLFiles/Math_P3S2_45.gif]

We assume that the value of k stays about the same,  k = 0.6. The number is entered exactly (no decimal) in order to have a nice looking solution formula.

In[135]:=

k = 6/10;
DEheatedhouse = T'[t] == k*(45+6*t-T[t])+H[t]
soln = DSolve[ {DEheatedhouse, T[0]==45}, T, t]

Out[136]=

T^′[t] 3 - 3 Cos[π t] + 3/5 (45 + 6 t - T[t])

Out[137]=

{{TFunction[{t}, 1/(9 + 25 π^2) (^(-3 t/5) (90 + 360 ^(3 t/5) + 1 ... /5) π^2 t - 45 ^(3 t/5) Cos[π t] - 75 ^(3 t/5) π Sin[π t]))]}}

The next input defined the solution function as HH and displays a simplified formula.

Note that Simplify is applied to N[...] by entering N[...]//Simplify.

In[140]:=

HH := soln[[1,1,2]]
N[HH[t],3]//Simplify

Out[141]=

40.0 + 5.18 ^(-0.60 t) + 6.00 t - 0.176 Cos[3.14 t] - 0.921 Sin[3.14 t]

The solution has an exponentially decaying term, a linear term increasing at the same rate as the ambient  temperature, and a sinusiodal component that oscillates with the same period as the furnace.

In[144]:=

Plot[ {45+6*t,HH[t]}, {t,0,5}, PlotRange->{40,75},
      PlotStyle->{Dashing[{0.02,0.02}],Dashing[{0,0}]}]

[Graphics:HTMLFiles/Math_P3S2_49.gif]

At 12 noon the house has warmed up to 70.4 degrees.

In[146]:=

HH[5]//N

Out[146]=

70.4337

The plot below shows the input heat function  H(t) and the sinusoidal component of the room temperature.

In[150]:=

Plot[ {H[t],-0.176 Cos[3.14 t]-0.921 Sin[3.14 t]}, {t,0,5},
      PlotStyle->{RGBColor[1,0,0],RGBColor[0,0,1]}, AspectRatio->1/3]

[Graphics:HTMLFiles/Math_P3S2_51.gif]

This picture confirms that the house temperature cycles up and down with the same period as the furnace.  However, the house temperature lags behind the heater by almost exactly a quarter cycle so that the room  temperature is increasing at its maximum rate at the same time that the furnace is putting out heat at its  maximum rate.

Go the Next Section Previous Section


Created by Mathematica  (September 13, 2004)