Chapter 9	Introduction to Differential Equations

9.1  Solving Differential Equations

Students should read Section 9.1 of Rogawski's Calculus [1] for a detailed discussion of the material presented in this section.

An ordinary differential equation is an equation that involves an unknown function, its derivatives, and an independent variable. Differential equations are useful for modeling many physical phenomena some of which are discussed in the next section.

Given a differential equation, our objective is to find all functions that satisfy it.  Mathematica's command for solving a differential equation is DSolve[eqn,y[x],x] where eqn is the differential equation to be solved and y[x] is the unknown function that depends on the independent variable x.

If the differential equation has initial conditions, we use braces { } to group the equation as well as the initial conditions (separated by commas): DSolve[{eqn,cond1,cond2,...,condn},y[x],x], where cond1, cond2,...,condn are initial conditions.

9.1.1. Separation of Variables

As discussed in your textbook, there is a special class of first-order differential equations that can be solved by hand using the method of separation of variables.  Mathematica can help in applying this method but of course it can solve the differential equation outright.  This makes Mathematica useful for verifying solutions obtained by other methods or for solving more complicated differential equations.  Since your textbook focuses on solving differential equations by hand, we will primarily discuss how to solve them using Mathematica.

Example 9.1.  Solve the given differential equation and plot the graph of the solutions.

a)  y'=2(4-y), y(0)=1		b)  (Sqrt[1-(x^2)])y'=x y		 c) y ((d y)/(d x)) +5x =0

Solution:

a) This is an initial value problem. Let us first solve this differential equation by hand using the method of separation of variables:
				(dy/dx)=2(4-y)	   ⟹ 	(dy/(y-4))=-2dx
							∫(dy/(y-4))=-∫2dx
							log|y-4|= -2x+C
							|y-4|= (e^(-2x+C))=(e^C)(e^(-2x))
							y-4=±(e^C)(e^(-2x))=C (e^(-2x)) (±(e^C) replaced by C)
This shows that the general solution is given by

y=C (e^(-2x))+4

It remains to determine the value of the constant C using the initial condition y(0)=1 (recall from your textbook that each value of C corresponds to a particular solution):

1=y(0)=C (e^(-2·0))+4=C+4

Thus, C=-3 and the unique solution is

y=-3 (e^(-2x))+4

Next, let us confirm this solution using Mathematica.  Recall that when entering a differential equation in Mathematica, we write y[x] instead of y to make explicit the dependence on x.

sola=DSolve[{y'[x]==2(4-y[x])},y[x],x]

This solution agrees with the solution obtained earlier by hand (the arbitrary constant C[1] is the same to the constant C).  We can visualize the behavior of the particular solutions by plotting some of their graphs for different values of C[1].  First, let us define the general solution to be g(x,c), where c=C[1] as follows (see Section 1.2.3 to learn how to extract elements from lists):

Clear[g,x,c]
g[x_,c_]=sola[[1,1,2]]/.C[1]->c

We then plot the one-parameter family of solution curves by combining the graphs of g(x,c) for c=-5,-4,...,5.

plotgeneralsolution=Plot[Table[g[x,c],{c,-5,5}],{x,-2,2},PlotRange->{-20,20},ImageSize->250]

Can you explain how the graph of g(x,c) varies as c varies?  Which c value corresponds to the top graph?

Next, to find the unique particular solution satisfying the given initial condition y(0)=1, we solve the equation g(0,c)=1 for c:

Solve[g[0,c]==1,c]

Thus, our unique solution is y=-3 (e^(-2x))+4.  This agrees with the solution we obtained earlier by hand.  Of course, Mathematica can solve for the unique solution on its own, bypassing the algebraic steps involved:

sola=DSolve[{y'[x]==2(4-y[x]), y[0]==1},y[x],x]

However, this unique solution does not appear to be the same as the one we obtained earlier.  To remedy this, let us extract solution from the output and define it as y=f(x):

f[x_]= sola[[1,1,2]]

We then apply the Expand command to simplify f(x):

Expand[f[x]]

Thus, the unique solution obtained by Mathematica is the same as the one obtained by hand.  Here is the plot of the unique solution:

plotuniquesolution=Plot[f[x], {x,-2,2},PlotRange->{-20,20},PlotStyle->Thickness[0.01],ImageSize->250]

Lastly, we combine plots of the general solution and the unique solution to show where the latter (bold graph) is situated in the former:

Show[{plotgeneralsolution,plotuniquesolution},ImageSize->250]

b)  From this point on, we shall skip using the method of separation of variables, which we leave for the reader to employ, and proceed directly to solving all differential equations using Mathematica as in part a) above.

solb=DSolve[(Sqrt[1-(x^2)])y'[x]==x y[x],y[x],x]

Again, we can visualize the behavior of these particular solutions by plotting graphs of some particular solutions corresponding to different values of C[1].  As before, we define the general solution to be g(x,c), where c=C[1].

Clear[g,x,c]
g[x_,c_]=solb[[1,1,2]]/.C[1]->c

We then make a combined plot of the graphs of g(x,c) for c=-5,-4,...,5.

Plot[Table[g[x,c],{c,-5,5}],{x,-2,2},ImageSize->250]

c) We again use Mathematica to directly obtain the solution:

Clear[y]
solde=DSolve[ y[x] y'[x]-5 x==0,y[x],x]

Observe that the two solutions, which we denote by f(x,c) and g(x,c), differ only in sign:

f[x_,c_]=solde[[2,1,2]]/.C[1]->c
g[x_,c_]=solde[[1,1,2]]/.C[1]->c

The following two plots show the graphs of f(x,c) and g(x,c) corresponding to c=-50, -40,...,0,..., 40,50.

Plot[Table[f[x,c],{c,-50,50,10}],{x,-5,5},PlotRange->{0,10},ImageSize->250]

Plot[Table[g[x,c],{c,-50,50,10}],{x,-5,5},PlotRange->{-10,0},ImageSize->250]

Observe that the two solutions y=-(Sqrt[5(x^2)+2c]) and y=(Sqrt[5(x^2)+2c])can be represented by a single equation:

(y^2)-5(x^2)=2c

which describes a family of hyperbolas. Here is a contour plot of this equation.  Observe that it nothing more than a combination of the two plots above as to be expected.

ContourPlot[(y^2)-5(x^2),{x, -5,5}, {y,-10,10}, Frame->False, Axes->True,ContourShading->False,Contours->10,ImageSize->250]

Exercises

In Exercises 1 through 8, solve the given differential equations.  If initial conditions are also given, then plot the unique solution.  If not, then make a combine plot of several particular solutions by choosing various values of the arbitrary constant.  Then describe the graphs and explain how they vary as the arbitrary constant varies.

1.  (1+(x^2))y'=(x^(2 ))y; y(0)=2		2.  y'+3(x^4)(y^2)=0; y(0)=1
3.  y'+(y^2)=-1; y(0)=-1			4.  y'+3 y=sin x; y(0)=0

5.  y'=-2x y (bell-shaped curves)		6.  16y y'+9x=0
7.  y'- y=(y^2)					8. 2x y y'-(y^2)+(x^2)=0

9. Consider the differential equation

(3+2 y)y'=2-(e^x),  y(0)=a

a) Solve the equation.
b) Plot the graphs for values of a=-2,-1, 0, 1, 2.
c) Plot the graphs for the values of a=-.5, -.1, .1,.5.
NOTE: For parts b) and c), make sure to use a sufficiently large interval for x.

10. Consider the differential equation

y=x y (b-y)/(4+x),    y(0)=a

a) Solve the equation.
b) Plot the graphs for values of a=-2,-1, 0, 1, 2  and b=-2, -1, 0, 1, 2.
c) Plot the graphs for the values of a=-.5, -.1, .1,.5 and b=-.5, -.1, .1,.5
d) Show that the limit as x-> ∞ of the solution does not depend on a. Does the limit depend on b? If so, how?

11. Suppose a skydiver falls from rest toward the earth and assume that the air resistance caused by his open parachute is  proportional to the square of his velocity v with proportionality constant k (we neglect air resistance due to the skydiver himself).  A model for describing the skydiver's velocity after his parachute opens is then given by the differential equation

v'=-(k/m)((v^2)-((m g)/k))

where m is the mass of the skydiver and g=9.8 meters/(sec^2) is his acceleration due to gravity.

a) Solve the equation assuming an initial velocity v(0)=(Subscript[v, 0]).
b) Suppose that for a particular skydiver m=70 kg and k=30 kg/meter.  Solve the equation again using these values and plot the particular solutions for the following values of (Subscript[v, 0]): 0, 2, ..., 10.
c) What is the skydiver's limiting (terminal) velocity as t->∞ for each of the particular solutions in part b)?  Does it depend on (Subscript[v, 0])?
d) Find a formula for the terminal velocity in terms of m, g, and k.

12. Recall that the first-order linear differential equation y'+y=0 has solution y=C (e^(-x)).  Solve the following higher-order generalizations of this equation:

a) y''+2y'+y=0
b) y'''+3y''+3y'+y=0
c) (y^((4)))+4y'''+6y''+4y'+y=0
d) Do you recognize the coefficients involved in the differential equations above?  What would be the next differential equation (of order 5) that follows this pattern?  Solve this differential equation to verify that its solution follows that same pattern exhibited in parts a) through c).

9.2  Models of the Form y'=k(y-b)

Students should read Section 9.2 of Rogawski's Calculus [1] for a detailed discussion of the material presented in this section.

NOTE: The differential equations we encounter in this section can be solved by the method of separation of variables and is discussed in the text. We leave it to the reader to solve the examples in this section by hand to verify the solutions obtained using Mathematica.

9.2.1. Bacteria Growth

The growth of bacteria in a culture is known to be proportional to the amount of the bacteria present at time t.  Suppose the initial amount of the bacteria is (Subscript[y, 0]) and the amount at time t is y(t). Then the above physical law is modeled by the differential equation

y'=k y,           y(0)=(Subscript[y, 0])

where k is the proportionality (growth) constant. Such a model exhibits exponential growth as can be seen from its solution below:

Clear[k]
DSolve[{y'[x]==k*y[x],y[0]==y0},y[x],x]

NOTE: Since the bacteria is growing in number, y(t) is increasing and hence y'(t)>0. Thus, k must be a positive number.

Example 9.2.  Suppose the amount of bacteria in a culture was 200 at time t=0.  It was found that there were 450 bacteria after 2 minutes.
a) Find the amount of the bacteria at any time t.   
b) At what time will the number of bacteria exceed 10,000?

Solution:

a) First, note that y(0)=200 and y(2)=450. We solve the differential equation y'=k y with the former as the initial condition:

Clear[y,t,k]
solde=DSolve[{y'[t]==k y[t],y[0]==200},y[t],t]

f[t_]=solde[[1,1,2]]

To find the value of k we solve f(2)=450 for k.

solk=Solve[f[2]==450,k]

N[%]

Thus, the proportionality constant is k=(1/2)l n(9/4)≈0.405465. Substituting this value into y(t), we see that the amount of bacteria at a given time t is

y(t)=200 (e^(0.405465 t))

b) To find the amount of time it takes for the bacteria to exceed 10,000, we solve

k=solk[[1,1,2]]
Solve[f[t]==10000,t]

We can approximate this value for t by

N[%]

Thus, it takes about 9.64824 minutes for the bacteria to reach 10,000. To visually see this, we plot the graphs of the solution y(t)=200 (e^(0.405465 t)) (blue curve) and y=10000 (red line) on the same axes.

Plot[{f[t],10000}, {t, 0, 15},PlotStyle->{Blue,Red},ImageSize->250]

NOTE: The solution y(t)=200 (e^(0.405465 t)) is only approximate since we approximated the growth constant k.  By using the exact value for k=(1/2)l n(9/4)=l n(3/2), we can derive the exact solution:

y(t)=200 (e^(k t))=200 (e^( t ln(3/2)))=200 (e^( ((ln((3/2)))^t)))=200((((3/2)))^t)

This agrees with the answer obtained by Mathematica:

f[t]

9.2.2. Radioactive Decay

The differential equation y'=k y is also used to model the amount of a radioactive substance whose rate of decay is proportional to the amount present.  However, in this case we note that the proportionality constant k<0. (Explain this!)

Example 9.3.  Carbon dating is a method used to determine the age of a fossil based on the amount of radioactive Carbon-14 in it compared to the amount normally found in the living environment. Suppose that a bone fossil contains 5% of the amount of Carbon-14 normally found in living animals. If the half-life of Carbon-14 is 5600 years, estimate the age of the bone.

Solution: Let y(t) be the amount of Carbon-14 in the bone and let (Subscript[y, 0]) be the initial amount of Carbon-14. Then the differential equation we need to solve is

Clear[k,y,y0]
solde=DSolve[{y'[t]==k y[t],y[0]==y0} , y[t],t]

Thus, the solution to the differential equation is y(t)=(Subscript[y, 0])(e^(k t)).  The half-life of Carbon-14 is 5600 implies that y(5600)=(1/2)(Subscript[y, 0]). We solve this equation for k:

y[t_]=solde[[1,1,2]]
solk=Solve[y[5600]==(1/2)y0,k]

N[%]

Thus, k=-0.000123776. To find the age of the bone, we solve y(t) = 0.05(Subscript[y, 0]) (5% of the initial amount) for t.

k=-((Log[2])/5600);
Solve[y[t]==0.05y0,t]

Thus, the bone is about 24,203 years old. Observe that it not necessary to know the original amount (Subscript[y, 0]) of Carbon-14 in the bone.

9.2.3. Annuity

An annuity is an  investment in which a principal amount of money is placed in a bank account that earns interest at an annual rate (compounded continuously) and the money is withdrawn at a regular interval. The differential equation that models an annuity is given by the annuity equation (rate of change = growth due to interest - withdrawal rate):

P'(t)=r P(t)-W =r(P(t)-(W/r))

where P(t) is the balance in the annuity at time t, r is the interest rate, and W is the rate (dollars per year) at which money is withdrawn continuously.

Example 9.4.  Find the general solution of the annuity equation for P(t) and then use it to calculate the following:
a) Assume r=6% and W=$ 6000 per year and P(0)=$ 50000.  Find P(t) and determine if and when the annuity runs out of money.
b) Assume r=6% and W=$ 6000 per year and P(0)=$ 100000.  Find P(t) and determine if and when the annuity runs out of money.
c) Assume r=6% and W=$ 12000 per year.  If we want the annuity to run out of money after 20 years, how much should be invested now?

Solution: We solve

DSolve[P'[t]==r(P[t]-(W/r)),P[t], t]

Thus, the general solution is P(t)=W/r+c (e^(r t)).

a) We set r=0.06, W=6000, and solve the initial value problem:

Clear[r,W,P]
r=0.06;
W=6000;
solde=DSolve[{P'[t]==r(P[t]-(W/r)),P[0]==50000},P[t], t]

We then define P(t) to be the solution above and plot it to see when the money will run out.

P[t_]=solde[[1,1,2]]
Plot[P[t], {t, 0,15},ImageSize->250]

As the graph indicates, the money runs out after approximately 11.5 years. We can confirm this by solving P(t)=0:

NSolve[P[t]==0,t]

b) We repeat the procedure in part a) with the obvious modifications:

Clear[r,W,P]
r=0.06;
W=6000;
solde=DSolve[{P'[t]==r(P[t]-(W/r)),P[0]==100000},P[t], t];
P[t_]=solde[[1,1,2]]
Plot[P[t], {t, 0, 80},ImageSize->250]

Observe that the balance P(t)=100,000 remains constant (can you explain why?) and thus the account will never run out of money.  What happens if we invest $100,001? $99,999?

c) In this case, we have r=0.06 and W=10000 per year.  The general solution is then given by

Clear[r,W,P,c]
r=0.06;
W=12000;
dsol=DSolve[{P'[t]==r(P[t]-(W/r)),P[0]==c},P[t], t];
P[t_]=dsol[[1,1,2]]

To determine the principal amount that will make the account run out of money in 20 years, we solve P(20)=0 for c:

NSolve[P[20]==0,c]

Thus, we need to invest $139,761.00 now.

9.2.4. Newton's Law of Cooling

Newton's Law of Cooling states that the rate of change in the temperature of an object is proportional to the difference between its temperature and that of the surrounding environment (known as the ambient temperature).  If A is the ambient temperature and T(t) is the temperature of the object at time t, then the differential equation that models this law is

T'(t)=-k(T(t)-A),   T(0)=(Subscript[T, 0])

where (Subscript[T, 0]) is the initial temperature of the object and k is a positive proportionality constant.

Example 9.6.  The temperature in an oven is (350^°) F when the oven is turned off. After 15 minutes, the temperature is (250^°) F.  Assume the temperature in the house is (70^°) F.

a) Find the temperature of the oven at any time t.
b) At what time will the temperature become (75^°) F? 
c) What will the temperature be in the limit as t-> ∞?
d) Does your answer in c) conform with your physical intuition?

Solution:

a) The ambient temperature here is room temperature.  Hence, A=70.  The initial temperature is (Subscript[T, 0])=350.  Newton's Law of Cooling then gives the model

T'(t)=-k(T(t)-70)),   T(0)=350

We solve this equation to get

Clear[T,k]
sol=DSolve[{T'[t]==-k (T[t]-70),T[0]==350},T[t],t]

T[t_]=sol[[1,1,2]]

Thus, the solution is T(t)=70(e^(-k t))(4+(e^(k t))) or T(t)=70+280 (e^(-k t)). To find the value of k, we solve T(15)=250 for k:

solk=Solve[T[15]==250,k]

k=solk[[1,1,2]]

N[%]

Thus, k=((ln(14/9))/15)=0.0294555.  Hence, the temperature of the oven at any time t is given by

T(t)=20+280(e^(-0.029455516818602614`t))

b) We solve T(t)=75 for t:

NSolve[T[t]==75, t]

Thus, the temperature will be (75^°) F after about two hours and 17 minutes.

c) We make a plot of the solution:

Plot[T[t], {t, 0, 100},AxesOrigin->{0,0},ImageSize->250]

To find the limiting temperature, we evaluate

Limit[T[t], t->Infinity]

d) Since heat flows from a region of higher temperature to a region of lower temperature, it is intuitively clear that the oven will cool down to the room (ambient) temperature.  Hence, the limit should be (70^°) F as expected.

Exercises

1. Mass of bacteria in a culture grow at a rate proportional to its size. Suppose the culture contains 200 cells intially and there are 800 cells after 3 hours.

a) Find the formula for the number of cells in the culture at time t.
b) Find the amount of bacteria after 2 hours.
c) At what time will the bacteria exceed 10,000 cells?

2. A mummy excavated from an archaelogical site in Egypt is found to contain 20% of Carbon-14 normally found in living humans.  Use carbon dating to estimate the age of the mummy.

3. Plutonium-239 is a highly radioactive element generated from waste in nuclear power plants and has a half-life of approximately 24,000 years.  How many years would it take for Plutonium-239 to decay to a safe level of 1/1000 its original amount?

4. Solve the following using the annuity differential eqaution discussed in this section.
a) Assume r=6% and W=$500 per year and P(0)=$5,000. Find P(t) and determine when the annuity runs out of money.
b) Assume r=6% and W=$500 per year and P(0)=$9,000. Find P(t) and determine when the annuity runs out of money.
c) Assume r=6% and W=$20,000 per year.  If we want the annuity to run out after 40 years, how much should we invest now?

5. Suppose a retired worker wants to invest in an annuity that will pay out $10,000 per year.
a) Assuming the annuity has an interest rate of 5%, find the minimum principal amount that should be invested so that the annuity never runs out of money.
b) Assuming the principal amount of money invested is $250,000, find the minimum interest rate that the annuity should bear so that it never runs out of money.

6. A hot metal rod is placed in a water bath whose temperature is  40° F. The rod cools from 300° F to 200° F in 1 minute. How long will take the rood to cool down to 150° F?  45° F?

9.3  Numerical Methods Using Slope Fields

Students should read Section 9.3 of Rogawski's Calculus [1] for a detailed discussion of the material presented in this section.

9.3.1.  Slope Fields

Consider a differential equation in the form

y'=f(x,y)

Since y' represents the slope of the line tangent to the graph of the solution y, we can think of f(x,y) as the slope of the same tangent line at the point (x,y), which we indicate by drawing a segment of it at the point of tangency.  The set of all such line segments (normalized to have the same length) is called the slope (or direction) field of the differential equation. Note that the slope field gives a graphical approximation to the solution.  It enables us to draw or visualize the graph of the unique solution of the equation passing through a given point. We will illustrate this in an upcoming example.

To plot the slope field of the differential equation y'=f(x,y) along the intervals (a,b) and (c,d) on the x- and y-axis, respectively, we use the command VectorPlot[{1, f[x, y]}, {x, a, b}, {y, c, d}], where slope is represented as a two-dimensional vector (1,f(x,y)) with the change in x normalized to equal 1.

NOTE: The command VectorPlot replaces the command VectorFieldPlot, which is obsolete in version 7 of Mathematica.

Example 9.9.  Consider the differential equation y'=(x^2)-2y, y(0)=-1.

a)  Draw the slope fields for the differential equation.
b)  Solve the differential equation.
c)  Plot both the slope field and the solution on the same axes.
d)  Redo parts b) and c) for the same equation but with initial condition given by y(a)=b.  Choose various values for a and b.

Solution:

a) Here, f(x,y)=(x^2)-2y.  We use the VectorPlot command to plot the corresponding slope field:

f[x_,y_]:=(x^2)-2y

plot1=VectorPlot[{1,f[x,y]},{x, -5,5},{y,-10,10},Axes->True,Frame->False,VectorScale->{Tiny,Tiny,None},ImageSize->250]

b) We use the DSolve command to find the exact solution of the differential equation.

Clear[y,x,g]
sol=DSolve[{y'[x]==f[x,y[x]],y[0]==-1},y[x],x]
g[x_]=sol[[1,1,2]]

c) We now plot the slope field together with the solution above passing through the point (0,-1):

plot2=Plot[g[x],{x,-5,5}, PlotRange->{-10,10},ImageSize->250]

Show[plot1,plot2,Graphics[{PointSize[Large],Point[{0,-1}]}],ImageSize->250]

d) We can show several graphs of solution curves (called integral curves) together with the corresponding slope field. Here is an example of how this can be done.

Clear[y,x,h,a,b]
sola=DSolve[{y'[x]==f[x,y[x]],y[a]==b},y[x],x];
h[x_,a_,b_]=Simplify[sola[[1,1,2]]]

plot3=Plot[Evaluate[Table[h[x,a,b],{a,-3,3,2},{b, -3,3,2}]],{x,-5,5},PlotRange->{-10,10}];
Show[plot1,plot3,ImageSize->250]

9.3.2.  Euler's Method

The simplest numerical method for solving a first order differential equation is Euler's Method. This method approximates the solution by moving along tangent lines described by the slope field of the differential equation.  Here is a brief description.

Let y=ϕ(x) be the solution of the differential equation

y'=f(x,y), y((Subscript[x, 0]))=(Subscript[y, 0])

Then the equation of the line tangent to the graph of y=φ(x) at x=(Subscript[x, 0]) is given by

y=φ'((Subscript[x, 0]))(x-(Subscript[x, 0]))+φ((Subscript[x, 0]))

But when x=(Subscript[x, 0]), we have φ((Subscript[t, 0]))=(Subscript[y, 0]) and φ'((Subscript[x, 0]))=f((Subscript[x, 0]),(Subscript[y, 0])). Thus, when x is close to (Subscript[x, 0]), φ(x) can be approximated by

y≈f((Subscript[x, 0]),(Subscript[y, 0]))(x-(Subscript[x, 0]))+(Subscript[y, 0])

We now choose h>0 to be a small positive number, called the step size, and define (Subscript[x, 1])=(Subscript[x, 0])+h. Then φ((Subscript[x, 1])) is approximately equal to

(Subscript[y, 1])=(Subscript[y, 0])+f((Subscript[x, 0]),(Subscript[y, 0]))((Subscript[x, 1])-(Subscript[x, 0]))

or

(Subscript[y, 1])=(Subscript[y, 0])+h f((Subscript[x, 0]),(Subscript[y, 0]))

We repeat the above argument at the point ((Subscript[x, 1]),(Subscript[y, 1])) to get an approximation of φ((Subscript[x, 2])), where (Subscript[x, 2])=(Subscript[x, 1])+h=(Subscript[x, 0])+2h:

(Subscript[y, 2])=(Subscript[y, 1])+h f((Subscript[x, 1]),(Subscript[y, 1]))

Proceeding in this manner, we obtain Euler's Method:

(Subscript[y, n+1])=(Subscript[y, n])+h f((Subscript[x, n]),(Subscript[y, n]))  for n=0, 1, 2, 3, ....

where φ((Subscript[x, n]))≈(Subscript[y, n]).

If the approximated solution is calculated over an interval [a,b] and the step size h is specified, then the number of iterations (or steps) required is given by m=(b-a)/h, where (Subscript[x, 0])=a and (Subscript[x, n])=(Subscript[x, 0])+n h.

Here is a Mathematica program called Euler for evaluating Euler's Method in m steps (the option SetPrecision sets the precision of our calculations to 10 digits).

Clear[f,x,y,x0,y0,h, m]
Euler[f_,h_,m_]:=Module[{n},
Do[
y[n+1]=SetPrecision[N[y[n]+h*f[x[n],y[n]]],10];
x[n+1]=x[n]+h,
{n,0,m}]]

Example 9.7.  Use the Euler program to construct a table of solution values for the differential equation y'=(x^2)+2y, y(0)=1 with a step size of h=0.1 and for m=10 steps.

Solution: Here f(x,y)=(x^2)+2y, (Subscript[x, 0])=0, (Subscript[y, 0])=1.

f[x_,y_]:=(x^2)+2y
m=10;
x0=0;
y0=1;
h=0.1;
x[0]=x0;
y[0]=y0;
Euler[f,h,m]
TableForm[Table[{n,x[n],y[n]},{n,1,m}],TableHeadings->{{},{"n ","(Subscript[x, n])","(Subscript[y, n])" }}]

To see how accurate the above approximation is, we solve the differential equation for the exact solution and plot both the approximate and exact solutions on the same axes.

Clear[z,t,exact]
exact=DSolve[{z'[t]==f[t,z[t]],z[x0]==y0},z[t],t];
z[t_]=z[t]/.exact[[1]]

Clear[plot1,plot2]
plot1=Plot[z[t],{t,0,1} ];
plot2=ListPlot[Table[{x[n],y[n]},{n,0,m}],PlotStyle->PointSize[0.01] ];
Show[{plot1,plot2},ImageSize->250]

Observe that the approximations become less accurate as we move away from the initial point (0,1). This is typical of numerical methods such as Euler's Method.  However, we can increase the accuracy of our approximation by decreasing the step size.  For example, we recompute the solution using h=0.05 (this increases the number of steps to m=20):

h=0.05;
m=20;
Euler[f,h,m]
TableForm[Table[{n,x[n],y[n]},{n,1,m}],TableHeadings->{{},{"n ","(Subscript[x, n])","(Subscript[y, n])" }}]

The following plot of the two numerical solutions corresponding to h=0.1 (small blue dots) and h=0.05 (large red dots) clearly shows that the latter is more accurate in comparison to the exact solution (curve):

plot3=ListPlot[Table[{x[n],y[n]},{n,0,m}],PlotStyle->{PointSize[0.01],Red,PointSize[0.015] }];
Show[{plot1,plot2,plot3},ImageSize->250]

Here is a modification of the Euler progam that allows the user to input the endpoints a and b directly (instead of the step size h) and m.

Clear[f,x,y,h,a,b,m]
EulerEndpt[f_,a_,b_,m_]:=Module[{n,h},h=N[((b-a)/m)];
Do[
y[n+1]=SetPrecision[N[y[n]+h*f[x[n],y[n]]],10];
x[n+1]=x[n]+h,
{n,0,m}]]

Example 9.8.  For the differential equation y'=(x^2)+2y, y(0)=1, approximate its solution over the interval [0,2] using m=10 steps.

Solution: Again, we have f(x,y)=(x^2)+y, (Subscript[x, 0])=0, (Subscript[y, 0])=1.  However, we now input the interval [a,b]=[0,2] into EulerEndPt.

f[x_,y_]:=(x^2)+y
m=10;
x0=0;
y0=1;
x[0]=x0;
y[0]=y0;
a=0;
b=2;

EulerEndpt[f,0,2,m]
TableForm[Table[{n,x[n],y[n]},{n,1,m}],
TableHeadings->{{},{"n ","(Subscript[x, n])","(Subscript[y, n])" }}]

This time we numerically compare the approximate solution with the exact solution:

Clear[z,t,exact]
exact=DSolve[{z'[t]==f[t,z[t]],z[x0]==y0},z[t],t];
z[t_]=z[t]/.exact[[1]]
TableForm[Table[{n,x[n],y[n],N[z[x0+n(b-a)/m]]},{n,1,m}],
TableHeadings->{{},{"n ","(Subscript[x, n])","(Subscript[y, n])","y(n)" }}]

Exercises

In Exercises 1 through 4, plot the slope field of the given differential equations:
1.  y'=(x^2)+(y^2)		2.  y'=(t^2)y		3.  y'=sin(x+y)		4.  y'=x (e^(-y))

5. Consider the differential equation  y'=3y-2(y^2).
a) Draw the slope field for the differential equation.
b) Solve the differential equation.
c) Assume y(0)=2.  Plot the graph of the solution for this case and also the slope field on the same axes.  Discuss the behavior of the solution as x->∞.
d)  Redo part c) for the same differential equation but with initial condition given by y(a)=b (choose various values for a and b).

In Exercises 6 through 9, use Euler's Method to find a numerical solution to the given initial value problem along the stated interval and using the prescribed number of steps.  Also, find their exact solutions and compare the results.
6.  y'=(x^2)-y,  y(0)=1; [0,1]; m=10	7.  y'=(1-(x^2)) cos y,  y(1)=0; [1,2]; m=20	8.  y'-(y^2)=x,  y(0)=1; [0,5]; m=50		9.  y'=-3x/2+ln((x^2)+y),  y(1)=1; [1,2]; m=100

10. Heun's method is a numerical method that improves on Euler's method.  It uses the approximation from Euler's method as an auxiliary value (called a predictor), which we denote by (Subscript[y, n+1]^( * )):

(Subscript[y, n+1]^( * ))=(Subscript[y, n])+h f((Subscript[x, n]),(Subscript[y, n]))

The actual approximation (called the corrector) is then computed as the mean of (Subscript[y, n+1]^( * )) (based on the slope of the tangent line at ((Subscript[x, n]),(Subscript[y, n]))) and (Subscript[y, n+1]^( ** ))=(Subscript[y, n])+h f((Subscript[x, n+1]),(Subscript[y, n+1]^( * )))  (based on the slope of the tangent line at ((Subscript[x, n+1]),(Subscript[y, n+1]^( * )))):

(Subscript[y, n+1])=(1/2)[(Subscript[y, n+1]^( * ))+(Subscript[y, n+1]^( ** ))] =(Subscript[y, n])+(1/2)[h f((Subscript[x, n]),(Subscript[y, n]))+h f((Subscript[x, n+1]),(Subscript[y, n+1]^( * )))]

a) Apply Heun's method to Exercise 6 and the compare the results obtained by Euler's method with the exact solution.  How much more accurate is Heun's method?
b) Redo part a) using m=20.  How much more accurate is the solution compared to that for m=10?

9.4  The Logistic Equation

Students should read Section 9.4 of Rogawski's Calculus [1] for a detailed discussion of the material presented in this section.

The differential equation

((d y)/(d t))=k y(1-(y/A))

is called the logistics equation. Here, k>0 and A is a constant called the carrying capacity.  This equation is useful for modeling the growth of a population where resources are limited and can only sustain a certain maximum population given by the carrying capacity.

Example 9.5.   The population p(t) of mosquito larvae growing in a tree hole increases according to the logistics equation with growth constant k=0.3 per day and carrying capacity A=1000.

a) Assuming that the initial population of the larvae is 50, find the population p(t) at any time t. 
b) After how many days will the larvae population exceed 500?
c) When does the larvae population reach 99% of the maximum capacity?

Solution:  
a) We use k=(3/10) and solve the corresponding differential equation in Mathematica:

Clear[y]
solde=DSolve[{y'[t]==(3/10) y[t](1-((y[t])/1000)), y[0]==50}, y[t],t]

NOTE:  Be careful with using a decimal approximation for k.  For example, try using k=0.3 and see what happens.

Next, for convenience we write the solution given in the previous output as

Clear[p,t]
p[t_]=solde[[1,1,2]]

Thus, the population of larvae at any time t is given by

p(t)=((1000 (e^(3 t/10)))/(19+(e^(3t/10))))

b) To find how long it takes for the larvae population to reach 500, we solve

NSolve[p[t]==500,t]

Thus, it takes about ten days for the larvae population to exceed 500.  Observe that p(10)≈513.887.

NOTE: We ignored the other two solutions in the previous output since they are complex-valued and not physically relevant.

c) We first plot the graph of p(t) to estimate the number of days required for the larvae population to reach 99% of the maximum capacity, that is, p(t)=999.

Plot[p[t], {t,0, 60},ImageSize->250]

It appears that the population reaches 999 larvae after t=30. We use the Table command to numerically confirm this.

TableForm[Table[{t,N[p[t],20]}, {t, 10, 50, 5}],TableHeadings->{{},{"Days     ", "Larvae Population"}}]

We can reasonably conclude from the table that the population reaches 999 larvae between 30 and 35 days. To obtain a more precise answer, we use Mathematica to solve p(t)=999 for t:

NSolve[p[t]==999,t]

Thus, the desired time is approximately t=33 days.

Exercises

1. A population of squirrels live in a forest with a carrying capacity of 3000. Assume logistic growth with growth constant k=0.8 per year.
a) Find the population of the squirrels at any time t assuming an initial population of P(0)=800.
b) How long will it take for the squirrel population to double?  Triple?

2. From 1960 to 2000, the world's population doubled from approximately 3 billion to 6 billion people.  Assuming that the human population follows logistic growth and that the earth has a carrying capacity of 100 billion people, determine the following:
a) Find the population at any time t.  What is the growth constant?
b) What will the population be in the year 2050?  How does this answer compare with projections from the United Nations?
c) Find the year in which the human population will reach half its carrying capacity?
d) When will the population grow the fastest, that is, the point on the graph where it changes from concave up to concave down (point of inflection)? What is the growth rate then?

3. In medicine, the logistics equation is used to model the growth of tumors in certain organs.  Suppose a patient is diagnosed with a tumor that has doubled in size to 2% of his liver when a year ago it only covered 1% of his liver.  How long will it take the tumor to grow to 10% of his liver when a transplant will most like by required to increase his chances of survival.

4. The current population of a herd of bison living inside a national park is 1000.  To ensure that the population does not reach more than 1500 bison in 50 years time and 2000 bison in 100 years time, what carrying capacity should the park maintain?  Assume that the bison population follows a logistics model.