Maple_P3S2.mw

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 Maple 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

diff(P(t), t) = r*P(t) .

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

> DE := diff(P(t),t) = 0.04*P(t);
principal := dsolve( {DE, P(0)=1000} );

DE := diff(P(t), t) = 0.4e-1*P(t)

principal := P(t) = 1000*exp(1/25*t)

> soln := fsolve( rhs(principal)=1500, {t} );
check, eval(principal,soln);

soln := {t = 10.13662770}

check, P(10.13662770) = 1500.000000

Answer: It will take about 10.14 years.

Let's plot the graph of P(t) including a point corresponding to P(t) = 1500.

> plot( [rhs(principal),[[10.14,1500]]], t=0..12, 0..1600,
      style=[line,point], color=[red,black], ytickmarks=4);

[Plot]

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.

> doubling_time = evalf(ln(2)/0.04)*years;

doubling_time = 17.32867952*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.

> DE2 := diff(y(t),t) = 0.09*y(t) - 1200;

DE2 := diff(y(t), t) = 0.9e-1*y(t)-1200

> debt := dsolve( {DE2, y(0)=2500} );

debt := y(t) = 40000/3-32500/3*exp(9/100*t)

The debt is liquidated when y(t) = 0.

> debt_free := fsolve( rhs(debt)=0, {t} );
check, eval(debt,debt_free);

debt_free := {t = 2.307104053}

check, y(2.307104053) = -0.1e-4

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

> subs(debt_free,12*t)*months;

27.68524864*months

> plot( rhs(debt), t=0..2.5, ytickmarks=3);

[Plot]

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

> total_payments = subs(debt_free,1200*t)*dollars;

total_payments = 2768.524864*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.

> DE3 := diff(x(t),t) = 0.02*x(t) - 1200;
savings := dsolve( {DE3, x(0)=35000} );

DE3 := diff(x(t), t) = 0.2e-1*x(t)-1200

savings := x(t) = 60000-25000*exp(1/50*t)

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

> eval(savings,debt_free);

x(2.307104053) = 33819.42022

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*exp(Float(9, -2)/12)-100 .

Remember that time is in years so one month is 1/12 years.

> d[2] := 2500*exp(0.09/12)-100;

d[2] := 2418.820488

Study the following for..do loop until you are convinced that it generates d[n] for n = 3, ... , 28. Output is suppressed.

> for n from 3 to 28
 do

  d[n] := d[n-1]*exp(0.09/12)-100:

 end do:

unassign('n'):

The next output shows that at the beginning of the 27th month, Juanita has paid $2600 and owes $178.22.

> d[27];

178.2163853

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

> d[28];

79.5580330

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

diff(T(t), 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

diff(T(t), t) = k*(A(t)-T(t))+H(t) .

If this equation is given to dsolve, Maple outputs the solution as a textbook integral. (See Ledder's Variation of Parameters formula, Chapter 4, page 260.)

> diff(T(t),t) = k*(A(t) - T(t)) + H(t);
dsolve( % );

diff(T(t), t) = k*(A(t)-T(t))+H(t)

T(t) = (Int(exp(k*t)*(k*A(t)+H(t)), t)+_C1)*exp(-k*t)

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.

> diff(T(t),t) = k*(A - T(t));
dsolve( % );

diff(T(t), t) = k*(A-T(t))

T(t) = A+exp(-k*t)*_C1

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.)

> DE_coffee := diff(T(t),t) = k*(66 - T(t));
coffee_temp := dsolve( {DE_coffee, T(0)=150} );

DE_coffee := diff(T(t), t) = k*(66-T(t))

coffee_temp := T(t) = 66+84*exp(-k*t)

The value of k can be found by substituting t = 5 and T(5) = 110 into coffee_temp and solving for k.

> subs(t=5,T(5)=110,coffee_temp);
k := solve( %, k );

110 = 66+84*exp(-5*k)

k := -1/5*ln(11/21)

The next calculation finds the time constant. (The k variable is now assigned.)

> evalf(1/k);

7.732431100

This tells us that if the coffee sits for 7.7 minutes, the difference (84 degrees) between T(0) and the ambient temperature will decrease to about 1/3 of its value. That is, the difference will be about 30 degrees: T(7.7) is roughly 96. We will find out when T(t) = 100 using fsolve.

> cools_to_100_in, fsolve( rhs(coffee_temp)=100, t)*minutes;

cools_to_100_in, 6.993645823*minutes

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

> plot( [rhs(coffee_temp),66,[[0,100],[7,100],[7,0]]], t=0..15, 0..160,
      linestyle=[SOLID$2,DASH], color=[red,blue,black],

      tickmarks=[3,4]);

[Plot]

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 unassign k.

> k := 'k':
DE_house := diff(T(t),t) = k*(45 + 6*t - T(t));

house_temp := dsolve( {DE_house, T(0)=45} );

DE_house := diff(T(t), t) = k*(45+6*t-T(t))

house_temp := T(t) = 45-6/k+6*t+6*exp(-k*t)/k

The value of k is obtained, as before, by substituting the extra temperature data into this equaton. Once this is done, Maple has a special function called LambertW that can provide an exact value for k.

> k_equation := subs(t=2,T(2)=50,house_temp);
k := solve( k_equation, k);

k_equation := 50 = 57-6/k+6*exp(-2*k)/k

k := 1/2*LambertW(-12/7*exp((-12)/7))+6/7


Let's turn
k into floating point form and check that it satisfies k_equation.

> k := evalf(k);
k_equation;

k := .5978788369

50 = 49.99999999

It would also be wise to verify that there is only one solution to k_equation. You are asked to do that in the exercises.

The next plot shows the ambient temperature (red) and the house temperature (blue) from 7 AM to 12 noon.

> plot( [rhs(house_temp), 45 + 6*t], t=0..5, 40..75, color=[blue,red]);

[Plot]

The temperature inside the house at noon is 65.5 degrees.

> eval(house_temp,t=5);

T(5) = 65.46948615

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*Pi*t/2)

> H := t -> 3 - 3*cos(Pi*t);
plot( H(t), t=0..5);

H := proc (t) options operator, arrow; 3-3*cos(Pi*t) end proc

[Plot]

We assume that the value of k stays about the same, k = 0.6. The exact solution formula is rather messy so the dsolve output is suppressed. The solution formula, evaluated to 3 digits, is displayed with the name soln.

> k := 0.6:
DE_heated_house := diff(T(t),t) = k*(45 + 6*t - T(t)) + H(t);

heated_house_temp := dsolve( {DE_heated_house, T(0)=45} ):

DE_heated_house := diff(T(t), t) = 30.0+3.6*t-.6*T(t)-3*cos(Pi*t)

> soln := evalf[3](heated_house_temp);

soln := T(t) = 5.0*exp(-.600*t)+40.0+6.00*t-.176*cos(3.14*t)-.925*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 sinuoidal componant that oscillates with the same period as the furnace.

> plot( [rhs(soln), 45 + 6*t], t=0..5, 40..75, color=[blue,red]);

[Plot]

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

> eval(soln,t=5);

T(5) = 70.41756381

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

> plot( [H(t), -.176*cos(3.14*t)-.925*sin(3.14*t)], t=0..5,
      color=[red,blue]);

[Plot]

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.

Orthogonal trajectories (Ledder, Section 2.2, Exercises)

Two families of curves are said to be orthogonal when they meet at right angles at all points of intersection. For example, the family of circles, centered at the origin, is orthogonal to the family of lines, passing throuth the origin. See below where the contourplot procedure (plots package) is used to create the circles. The plot procedure creates lines and all the curves are shown using the display procedure. The optional equation scaling=constrained is needed to guarantee that orthogonal curves appear to meet a right angles.

> with(plots):
Circles := contourplot( x^2 + y^2, x=-3..3, y=-3..3, color=red):

Lines := plot( [0.4*x, 0.7*x, 1.2*x, 2*x, 4*x, -0.4*x, -0.7*x, -1.2*x,

                -2*x, -4*x], x=-3..3, y=-3..3, color=blue):

display( Circles, Lines, view=[-3..3,-3..3], scaling=constrained);

[Plot]

Orthogonal trajectories can be found using diffferential equations. For example, circle equations,

x^2+y^2 = c

are the solution formulas for the differential equation

y' = -x/y

(differentiate the circle equation implicitly with respect to x and solve for y' ). Since the orthogonal trajectories meet the circles with perpendicular tangents they must satisfy the differential equation

y' = y/x

Why? You ask?

Answer: Because the slopes of perpendicular lines are negative reciprocals of one another.

Solving the last equation shows that the orthogonal trajectories are straight lines through the origin, as claimed.

> dsolve( diff(y(x),x) = y(x)/x );

y(x) = _C1*x

Another Example. Obtain the trajectories that are orthogonal to the family of hyperbolas determined by the equation

x^2-y^2 = c .

Solution. Differentiate implicitly with respect to x to see that the hyperbolas are solutions to the differential equation y' = x/y. Therefore, the orthogonal trajectories solve the equation y' = -y/x. These solutions are displayed below.

> soln := dsolve( diff(y(x),x) = -y(x)/x );

soln := y(x) = _C1/x

Thus, the orthogonal trajectories are the level curves for the function f (x, y) = x y . See the picture below.

> Hyper := contourplot( x^2 - y^2, x=-3..3, y=-3..3, color=red):
Orth := contourplot( x*y, x=-3..3, y=-3..3, color=blue):

display( Hyper, Orth, view=[-3..3,-3..3], scaling=constrained);

[Plot]

Go to: Next Section Previous Section