[Go to the Notes on Diffy Qs home page] [PDF version]

[next] [prev] [prev-tail] [tail] [up]

Chapter 1
First order ODEs

1.1 Integrals as solutions

Note: 1 lecture (or less), §1.2 in EP

A first order ODE is an equation of the form

dy-= f(x,y),
dx

or just

 ′
y = f (x,y).

In general, there is no simple formula or procedure one can follow to find solutions. In the next few lectures we will look at special cases where solutions are not difficult to obtain. In this section, let us assume that f is a function of x alone, that is, the equation is

y′ = f (x).
(1.1)

We could just integrate (antidifferentiate) both sides with respect to x .

∫            ∫
   y′(x) dx =    f(x) dx + C,

that is

       ∫

y(x) =    f(x) dx + C.

This y(x)  is actually the general solution. So to solve (1.1), we find some antiderivative of f (x )  and then we add an arbitrary constant to get the general solution.

Now is a good time to discuss a point about calculus notation and terminology. Calculus textbooks muddy the waters by talking about integral as primarily the so-called indefinite integral. The indefinite integral is really the antiderivative (in fact the whole one-parameter family of antiderivatives). There really exists only one integral and that is the definite integral. The only reason for the indefinite integral notation is that we can always write an antiderivative as a (definite) integral. That is, by the fundamental theorem of calculus we can always write ∫
  f (x ) dx + C as

∫  x
    f (t) dt + C.
  x0

Hence the terminology to integrate when we may really mean to antidifferentiate. Integration is just one way to compute the antiderivative (and it is a way that always works, see the following examples). Integration is defined as the area under the graph, it only happens to also compute antiderivatives. For sake of consistency, we will keep using the indefinite integral notation when we want an antiderivative, and you should always think of the definite integral.

Example 1.1.1: Find the general solution of y′ = 3x2   .

Elementary calculus tells us that the general solution must be y = x3 + C . Let us check: y′ = 3x2   . We have gotten precisely our equation back.

Normally, we also have an initial condition such as y(x0) = y0   for some two numbers x0   and y0   (x0   is usually 0, but not always). We can then write the solution as a definite integral in a nice way. Suppose our problem is y′ = f(x)  , y(x ) = y
   0    0   . Then the solution is

      ∫
         x
y(x) =    f(s) ds + y0.
        x0
(1.2)

Let us check! We compute y′ = f (x)  (by fundamental theorem of calculus) and by Jupiter, y is a solution. Is it the one satisfying the initial condition? Well, y(x ) = ∫ x0f (x) dx + y = y
   0    x0           0   0   . It is!

Do note that the definite integral and the indefinite integral (antidifferentiation) are completely different beasts. The definite integral always evaluates to a number. Therefore, (1.2) is a formula we can plug into the calculator or a computer, and it will be happy to calculate specific values for us. We will easily be able to plot the solution and work with it just like with any other function. It is not so crucial to always find a closed form for the antiderivative.

Example 1.1.2: Solve

 ′   -x2
y = e   ,    y(0) = 1.

By the preceding discussion, the solution must be

       ∫ x
y(x) =    e-s2 ds + 1.
        0

Here is a good way to make fun of your friends taking second semester calculus. Tell them to find the closed form solution. Ha ha ha (bad math joke). It is not possible (in closed form). There is absolutely nothing wrong with writing the solution as a definite integral. This particular integral is in fact very important in statistics.

Using this method, we can also solve equations of the form

 ′
y =  f(y).

Let us write the equation in Leibniz notation.

dy-= f(y).
dx

Now we use the inverse function theorem from calculus to switch the roles of x and y to obtain

dx    1
dy-= f(y).

What we are doing seems like algebra with dx and dy . It is tempting to just do algebra with dx and dy as if they were numbers. And in this case it does work. Be careful, however, as this sort of hand-waving calculation can lead to trouble, especially when more than one independent variable is involved. At this point we can simply integrate

       ∫   1
x(y) =    ----dy + C.
          f(y)

Finally, we try to solve for y .

Example 1.1.3: Previously, we guessed y′ = ky (for some k > 0  ) has solution Cekx  . We can actually do it now. First we note that y = 0  is a solution. Henceforth, we assume y ⇔ 0  . We write

dx    1
--- = --.
 dy   ky

We integrate to obtain

          1
x(y) = x =--ln|y| + D,
          k

where D is an arbitrary constant. Now we solve for y (actually for |y| ).

     kx-kD    -kD kx
|y| = e    = e   e  .

If we replace e-kD  with an arbitrary constant C we can get rid of the absolute value bars (we can do this as D was arbitrary). In this way, we also incorporate the solution y = 0  . We get the same general solution as we guessed before, y = Cekx  .

Example 1.1.4: Find the general solution of y′ = y2   .

First we note that y = 0  is a solution. We can now assume that y ⇔ 0  . Write

dx   1
---= --.
dy   y2

We integrate to get

    -1-
x =  y + C.

We solve for y = -1-
    C-x  . So the general solution is

y = --1---    or    y = 0.
    C - x

Note the singularities of the solution. If for example C = 1  , then the solution “blows up” as we approach x = 1  . Generally, it is hard to tell from just looking at the equation itself how the solution is going to behave. The equation  ′    2
y  = y   is very nice and defined everywhere, but the solution is only defined on some interval (-∞, C )  or (C,∞ )  .

Classical problems leading to differential equations solvable by integration are problems dealing with velocity, acceleration and distance. You have surely seen these problems before in your calculus class.

Example 1.1.5: Suppose a car drives at a speed et∕2   meters per second, where t is time in seconds. How far did the car get in 2 seconds (starting at t = 0  )? How far in 10 seconds?

Let x denote the distance the car traveled. The equation is

x′ = et∕2.

We can just integrate this equation to get that

        t∕2
x(t) = 2e  + C.

We still need to figure out C . We know that when t = 0  , then x = 0  . That is, x(0) = 0  . So

0 = x(0) = 2e0∕2 + C = 2 + C.

Thus C  = -2  and

x(t) = et∕2 - 2.

Now we just plug in to get where the car is at 2 and at 10 seconds. We obtain

x(2) = 2e 2∕2 - 2 ≈ 3.44 meters,  x (10 ) = 2e10∕2 - 2 ≈ 294 meters.

Example 1.1.6: Suppose that the car accelerates at a rate of 2m 2
t ∕s    . At time t = 0  the car is at the 1 meter mark and is traveling at 10 m/s. Where is the car at time t = 10  .

Well this is actually a second order problem. If x is the distance traveled, then x′ is the velocity, and x′′ is the acceleration. The equation with initial conditions is

 ′′   2                   ′
x  = t,     x(0) = 1,    x (0) = 10.

Well, what if we call x′ = v and then we have the problem

v′ = t2,     v(0 ) = 10.

Once we solve for v , we can then integrate and find x .

Exercise 1.1.1: Solve for v and then solve for x . Then find x (1 0)  to answer the question.

1.1.1 Exercises

Exercise 1.1.2: Solve dy = x2 + x
dx for y(1) = 3  .

Exercise 1.1.3: Solve dy = sin 5x
dx for y(0) = 2  .

Exercise 1.1.4: Solve dy   -1--
dx = x2-1   for y(0) = 0  .

Exercise 1.1.5: Solve y′ = y3   for y(0) = 1  .

Exercise 1.1.6 (challenging): Solve  ′
y = (y - 1)(y + 1 )  for y(0) = 3  .

Exercise 1.1.7: Solve dy = -1-
dx   y+1   for y(0) = 0  .

Exercise 1.1.8: Solve  ′′
y  = sin x for y(0) = 0  .

Exercise 1.1.9: A spaceship is traveling at the speed 2t2 + 1km∕s   (t is time in seconds). It is pointing directly away from earth and at time t = 0  it is 1000 kilometers from earth. How far from earth is it one minute from time 0?

Exercise 1.1.10: Solve dxdt = sin(t2) + t , x(0) = 20  . It is OK to leave your answer as a definite integral.

1.2 Slope fields

Note: 1 lecture, §1.3 in EP

At this point it may be good to first try the Lab I and/or Project I from the IODE website: http://www.math.uiuc.edu/iode/.

As we said, the general first order equation we are studying looks like

y′ = f (x,y).

In general, we cannot simply solve these kinds of equations explicitly. It would be good if we could at least figure out the shape and behavior of the solutions, or if we could even find approximate solutions for any equation.

1.2.1 Slope fields

As you have seen in IODE Lab I (if you did it), the equation y′ = f(x,y)  gives you a slope at each point in the (x,y)  -plane. We can plot the slope at lots of points as a short line with this given slope. See Figure 1.1.


PIC

Figure 1.1: Slope field of  ′
y = xy .
PIC
Figure 1.2: Slope field of  ′
y =  xy with a graph of solutions satisfying y(0) = 0.2  , y(0) = 0  , and y (0) = -0.2  .


We call this picture the slope field of the equation. If we are given a specific initial condition y(x0) = y 0   , we can look at the location (x0,y 0)  and follow the slopes. See Figure 1.2.

By looking at the slope field we can get a lot of information about the behavior of solutions. For example, in Figure 1.2 we can see what the solutions do when the initial conditions are y(0) > 0  , y(0) = 0  and y(0) < 0  . Note that a small change in the initial condition causes quite different behavior. On the other hand, plotting a few solutions of the equation y′ = - y , we see that no matter what y(0)  is, all solutions tend to zero as x tends to infinity. See Figure 1.3.


PIC

Figure 1.3: Slope field of  ′
y = - y with a graph of a few solutions.


1.2.2 Existence and uniqueness

We wish to ask two fundamental questions about the problem

y′ = f(x,y),    y(x ) = y .
                   0    0

(i)
Does a solution exist?
(ii)
Is the solution unique (if it exists)?

What do you think is the answer? The answer seems to be yes to both does it not? Well, pretty much. But there are cases when the answer to either question can be no.

Since generally the equations we encounter in applications come from real life situations, it seems logical that a solution always exists. It also has to be unique if we believe our universe is deterministic. If the solution does not exist, or if it is not unique, we have probably not devised the correct model. Hence, it is good to know when things go wrong and why.

Example 1.2.1: Attempt to solve:

y′ = 1-,   y(0) = 0.
    x

Integrate to find the general solution y = ln|x| + C . Note that the solution does not exist at x = 0  . See Figure 1.4.


PIC

Figure 1.4: Slope field of y′ = 1∕x  .
PIC
Figure 1.5: Slope field of       ∘--
y′ = 2 |y| with two solutions satisfying y(0) = 0  .


Example 1.2.2: Solve:

 ′    ∘ --
y = 2   |y|,     y(0) = 0.

See Figure 1.5. Note that y = 0  is a solution. But also the function

       (
       ||{x2     if x ≥ 0,
y(x) = ||(   2
        - x    if x < 0.

It is actually hard to tell from the slope field that the solution will not be unique. Is there any hope? Of course there is. It turns out that the following theorem is true. It is known as Picard’s theorem* .

Theorem 1.2.1 (Picard’s theorem on existence and uniqueness). If f (x, y)  is continuous (as a function of two variables) and ∂f∂y  exists and is continuous near some (x0,y0)  , then a solution to

y′ = f(x,y),   y(x0) = y0,

exists (at least for some small interval of x ’s) and is unique.

Note that the problems  ′  1
y =  ∕x  , y(0) = 0  and  ′    ∘--
y = 2  |y| , y(0 ) = 0  do not satisfy the hypothesis of the theorem. Even if we can use the theorem, we ought to be careful about this existence business. It is quite possible that the solution only exists for a short while.

Example 1.2.3: For some constant A , solve:

 ′    2
y  = y ,    y(0) = A.

We know how to solve this equation. First assume that A ⇔ 0  , so y is not equal to zero at least for some x near 0. So  ′  1 2
x =  ∕y    , so     -1
x =  ∕y + C , so     -1-
y = C-x  . If y(0) = A , then     1
C =  ∕A  so

    ---1--
y = 1∕A - x.

If A =  0  , then y = 0  is a solution.

For example, when A = 1  the solution “blows up” at x = 1  . Hence, the solution does not exist for all x even if the equation is nice everywhere. The equation y′ = y2   certainly looks nice.

For the most of this course we will be interested in equations where existence and uniqueness holds, and in fact holds “globally” unlike for the equation y′ = y2   .

1.2.3 Exercises

Exercise 1.2.1: Sketch direction field for y′ = ex-y  . How do the solutions behave as x grows? Can you guess a particular solution by looking at the direction field?

Exercise 1.2.2: Sketch direction field for y′ = x2   .

Exercise 1.2.3: Sketch direction field for y′ = y2   .

Exercise 1.2.4: Is it possible to solve the equation  ′   xy--
y =  cosx  for y(0) = 1  ? Justify.

Exercise 1.2.5: Is it possible to solve the equation       √---
y′ = y |x| for y(0) = 0  ? Is the solution unique? Justify.

1.3 Separable equations

Note: 1 lecture, §1.4 in EP

When a differential equation is of the form y′ = f(x)  , we can just integrate:     ∫
y =  f (x )dx + C . Unfortunately this method no longer works for the general form of the equation y′ = f(x,y)  . Integrating both sides yields

    ∫
y =   f (x,y)dx + C.

Notice dependence on y in the integral.

1.3.1 Separable equations

Let us suppose that the equation is separable. That is, let us consider

y′ = f (x)g (y),

for some functions f(x)  and g(y)  . Let us write the equation in the Leibniz notation

dy-= f (x)g(y).
dx

Then we rewrite the equation as

-dy-
g(y) = f(x)dx.

Now both sides look like something we can integrate. We obtain

∫   dy    ∫
   ---- =   f (x)dx + C.
   g(y)

If we can find closed form expressions for these two integrals, we can, perhaps, solve for y .

Example 1.3.1: Take the equation

y′ = xy.

First note that y = 0  is a solution, so assume y ⇔ 0  from now on. Write the equation as dy= xy
dx , then

∫       ∫
   dy-
   y  =    x dx + C.

We compute the antiderivatives to get

        2
ln |y| = x-+ C.
       2

Or

|y| = ex22 +C = ex22 eC = De x22 ,

where D  > 0  is some constant. Because y = 0  is a solution and because of the absolute value we actually can write:

y = De x22 ,

for any number D (including zero or negative).

We check:

        x2    (   x2)
y′ = Dxe 2-= x De 2- = xy.

Yay!

We should be a little bit more careful with this method. You may be worried that we were integrating in two different variables. We seemed to be doing a different operation to each side. Let us work this method out more rigorously.

dy-= f (x)g (y)
dx

We rewrite the equation as follows. Note that y = y(x)  is a function of x and so is ddyx  !

-1--dy- = f(x)
g(y)dx

We integrate both sides with respect to x .

∫               ∫
    -1--dy-
    g(y) dx dx =    f (x ) dx + C.

We can use the change of variables formula.

∫    1       ∫
   ---- dy =    f(x) dx + C.
   g(y)

And we are done.

1.3.2 Implicit solutions

It is clear that we might sometimes get stuck even if we can do the integration. For example, take the separable equation

 ′    xy
y =  -2---.
     y + 1

We separate variables

 2          (     )
y--+-1 dy = y + 1- dy = x dx.
   y            y

Now we integrate to get

y2         x2
--+ ln|y| =-- + C.
2           2

Or maybe the easier looking expression (where D = 2C ):

 2            2
y + 2 ln |y| = x + D.

It is not easy to find the solution explicitly as it is hard to solve for y . We will, therefore, call this solution an implicit solution. It is easy to check that implicit solutions still satisfy the differential equation. In this case, we differentiate to get

  (      )
 ′      2
y  2y + --= 2x.
        y

It is simple to see that the differential equation holds. If you want to compute values for y you might have to be tricky. For example, you can graph x as a function of y , and then flip your paper. Computers are also good at some of these tricks, but you have to be careful.

We note above that the equation also has a solution y = 0  . In this case, it turns out that the general solution is  2            2
y  + 2ln|y| = x + C together with y = 0  . These outlying solutions such as y = 0  are sometimes called singular solutions.

1.3.3 Examples

Example 1.3.2: Solve x2y′ = 1 - x2 + y2 - x 2y2   , y(1 ) = 0  .

First factor the right hand side to obtain

x2y′ = (1 - x2)(1 + y2).

We separate variables, integrate and solve for y

    y′     1 - x2
  -----2 = ---2--,
  1 + y      x
  --y′--   1-
  1 + y2 = x2 - 1,

arctan(y) = -1-- x + C,
            x (          )
               -1-
       y = tan   x - x + C  .
Now solve for the initial condition, 0 = tan (- 2 + C )  to get C = 2  (or 2 + π , etc…). The solution we are seeking is, therefore,
       (          )
        -1-
y = tan x  - x + 2 .

Example 1.3.3: Suppose Bob made a cup of coffee, and the water was boiling (100 degrees Celsius) at time t = 0  minutes. Suppose Bob likes to drink his coffee at 70 degrees. Let the ambient (room) temperature be 26 degrees. Furthermore, suppose Bob measured the temperature of the coffee at 1 minute and found that it dropped to 95 degrees. When should Bob start drinking?

Let T be the temperature of coffee, let A be the ambient (room) temperature. Then for some k the temperature of coffee is:

dT
---= k(A - T).
dt

For our setup A = 26  , T(0) = 100  , T (1) = 95  . We separate variables and integrate (C and D will denote arbitrary constants)

--1---dT- = -k,
T - A  dt
ln(T - A) = -kt + C, (note that T - A > 0)
              -kt
    T - A = De   ,
        T = De-kt + A.
That is, T = 26 + De -kt  . We plug in the first condition 10 0 = T (0) = 2 6 + D and hence D  = 74  . Now we have             -kt
T = 26 + 74e  . We plug in                     -k
95 = T (1) = 26 + 74e  . Solving for k we get         95-26-
k = - ln 74 ≈ 0.07  . Now to solve for t which gives us 70 degrees. That is, we solve              -0.07t
70 = 26 + 74e  to get        70-26-
t = - ln0.7074 ≈ 7.43  minutes. So Bob can begin to drink the coffee at about 7 and a half minutes from the time Bob made it. Probably about the amount of time it took us to calculate how long it would take.

Example 1.3.4: Find the general solution to y′ = -xy2
     3   (including singular solutions).

First note that y = 0  is a solution (a singular solution). So assume that y ⇔ 0  and write

-3- ′
y2 y = x,
         2
   3-= x- + C,
   y    2
          3         6
   y = x2----- = x-2 +-2C-.
         ∕2 + C

1.3.4 Exercises

Exercise 1.3.1: Solve y′ = x∕y  .

Exercise 1.3.2: Solve y′ = x2y .

Exercise 1.3.3: Solve dx-= (x2 - 1)t
dt , for x(0) = 0  .

Exercise 1.3.4: Solve dx-
dt = x sin (t)  , for x(0) = 1  .

Exercise 1.3.5: Solve dy
---= xy + x + y + 1
dx  . Hint: Factor the right hand side.

Exercise 1.3.6: Solve xy′ = y + 2x2y , where y(1) = 1  .

Exercise 1.3.7: Solve        2
dy-   y-+-1-
dx =  x2 + 1  , for y(0) = 1  .

Exercise 1.3.8: Find an implicit solution for dy    x2 + 1
---=  -2----
dx    y + 1  , for y(0) = 1  .

Exercise 1.3.9: Find explicit solution for  ′    -y
y = xe  , y(0) = 1  .

Exercise 1.3.10: Find explicit solution for   ′   -y
xy = e  , for y(1) = 1  .

Exercise 1.3.11: Find explicit solution for y′ = ye-x2    , y(0) = 1  . It is alright to leave a definite integral in your answer.

Exercise 1.3.12: Suppose a cup of coffee is at 100 degrees Celsius at time t = 0  , it is at 70 degrees at t = 10  minutes, and it is at 50 degrees at t = 20  minutes. Compute the ambient temperature.

1.4 Linear equations and the integrating factor

Note: 1 lecture, §1.5 in EP

One of the most important types of equations we will learn how to solve are the so-called linear equations. In fact, the majority of the course will focus on linear equations. In this lecture we will focus on the first order linear equation. A first order equation is linear if we can put it into the following form:

y′ + p (x)y = f (x).
(1.3)

The word “linear” here means linear in y . The dependence on x can be more complicated.

Solutions of linear equations have nice properties. For example, the solution exists wherever p(x)  and f (x)  are defined, and has the same regularity (read: it is just as nice). But most importantly for us right now, there is a method for solving linear first order equations.

First we find a function r(x)  such that

     ′              d-[    ]
r(x)y + r(x)p(x)y = dx r(x)y .

Then we can multiply both sides of (1.3) by r(x)  to obtain

  [     ]
d--r(x)y = r(x)f(x).
dx

Now we integrate both sides. The right hand side does not depend on y and the left hand side is written as a derivative of a function. Afterwards, we can solve for y . The function r(x)  is called the integrating factor and the method is called the integrating factor method.

We are looking for a function r(x)  such that if we differentiate it, we get the same function back multiplied by p(x)  . That seems like a job for the exponential function! Let

       ∫
r(x) = e p(x)dx

We compute:

              ′
             y + p(x)y = f(x),
 ∫ p(x)dx ′  ∫ p(x)dx         ∫ p(x)dx
e     y + e      p(x)y = e     f (x),
          d--[ ∫ p(x)dx ]   ∫ p(x)dx
          dx  e     y  = e     f (x),
               ∫         ∫   ∫
              e  p(x)dxy =    e p(x)dxf(x)dx + C,
                                 (∫                   )
                          -∫ p(x)dx    ∫ p(x)dx
                     y = e           e     f (x)dx + C  .

Of course, to get a closed form formula for y we need to be able to find a closed form formula for the integrals appearing above.

Example 1.4.1: Solve

y′ + 2xy = ex-x2,    y(0) = - 1.

First note that p(x) = 2x and            2
f (x ) = ex-x    . The integrating factor is         ∫         2
r(x) = e p(x)dx = ex    . We multiply both sides of the equation by r(x)  to get

 x2 ′     x2    x-x2 x2
e  y + 2xe  y = e   e ,
     d  [x2 ]   x
     ---e  y = e .
     dx
We integrate
  2
ex y = ex + C,
        x-x2     -x2
   y = e   + Ce   .
Next, we solve for the initial condition - 1 = y(0) = 1 + C , so C  = -2  . The solution is
       2      2
y = ex-x - 2e-x .

Note that we do not care which antiderivative we take when computing  ∫
e p(x)dx  . You can always add a constant of integration, but those constants will not matter in the end.

Exercise 1.4.1: Try it! Add a constant of integration to the integral in the integrating factor and show that the solution you get in the end is the same as what we got above.

An advice: Do not try to remember the formula itself, that is way too hard. It is easier to remember the process and repeat it.

Since we cannot always evaluate the integrals in closed form, it is useful to know how to write the solution in definite integral form. A definite integral is something that you can plug into a computer or a calculator. Suppose we are given

y′ + p(x)y = f(x),   y(x0) = y0.

Look at the solution and write the integrals as definite integrals.

|----------∫-----(∫--x-∫---------------)---|
|         -xx0p(s)ds       tx0p(s)ds             |
| y(x) = e          x e      f (t)dt + y0 . |
---------------------0---------------------
(1.4)

You should be careful to properly use dummy variables here. If you now plug such a formula into a computer or a calculator, it will be happy to give you numerical answers.

Exercise 1.4.2: Check that y(x0) = y0   in formula (1.4).

Exercise 1.4.3: Write the solution of the following problem as a definite integral, but try to simplify as far as you can. You will not be able to find the solution in closed form.

          2
y′ + y = ex-x,    y(0) = 10.

Example 1.4.2: The following is a simple application of linear equations and this type of a problem is used often in real life. For example, linear equations are used in figuring out the concentration of chemicals in bodies of water.

A 100 liter tank contains 10 kilograms of salt dissolved in 60 liters of water. Solution of water and salt (brine) with concentration of 0.1 kilograms per liter is flowing in at the rate of 5 liters a minute. The solution in the tank is well stirred and flows out at a rate of 3 liters a minute. How much salt is in the tank when the tank is full?

Let us come up with the equation. Let x denote the kilograms of salt in the tank, let t denote the time in minutes. Then for a small change Δt in time, the change in x (denoted Δx ) is approximately

Δx ≈ (rate in × concentration in)Δt - (rate out × concentration out)Δt.

Dividing through by Δt and taking the limit Δt → 0  we see that

dx-
dt = (rate in × concentration in) - (rate out × concentration out).

In our example, we have

          rate in = 5,

 concentration in = 0.1,
         rate out = 3,
                      x           x
concentration out = -------=  ------------.
                   volume    60 + (5 - 3)t
Our equation is, therefore,
                (        )
dx- = (5 × 0.1) - 3---x--- .
 dt               6 0 + 2t

Or in the form (1.3)

dx-+ ---3---x = 0.5.
dt   60 + 2t

Let us solve. The integrating factor is

          (∫          )       (           )
r(t) = exp    --3---dt  = exp  3ln(60 + 2t) = (60 + 2t)3∕2.
              60 + 2t          2

We multiply both sides of the equation to get

         3∕2dx-          3∕2--3----              3∕2
(60 + 2t)  dt + (60 + 2t)   60 + 2t x = 0.5(60 + 2t) ,
                   d  [       3∕2 ]             3∕2
                   -- (60 + 2t)  x = 0.5(60 + 2t)  ,
                   dt                ∫
                       (60 + 2t)3∕2x =   0.5(60 + 2t)3∕2dt + C,

                                                 ∫  (60 + 2t)3∕2
                                  x = (60 + 2t)-3∕2  -----------dt + C (60 + 2t)-3∕2,
                                                         2
                                             -3∕2-1-        5∕2            -3∕2
                                  x = (60 + 2t)  10(60 + 2t)  + C (60 + 2t)  ,
                                      60 + 2t
                                  x = -------+ C(60 + 2t)-3∕2.
                                        10

We need to find C . We know that at t = 0  , x = 1 0  . So

            60-       -3∕2            -3∕2
10 = x(0) = 10 + C(60)    = 6 + C(60)   ,

or

C  = 4(603∕2) ≈ 1859.03.

We are interested in x when the tank is full. So we note that the tank is full when 60 + 2t = 100  , or when t = 20  . So

x(20) = 60-+-40-+ C(60 + 40)-3∕2 ≈ 10 + 1859.03(100)-3∕2 ≈ 11.86.
           10

The concentration at the end is approximately 0.1186 kg/liter and we started with 1/6 or 0.167 kg/liter.

1.4.1 Exercises

In the exercises, feel free to leave answer as a definite integral if a closed form solution cannot be found. If you can find a closed form solution, you should give that.

Exercise 1.4.4: Solve y′ + xy = x .

Exercise 1.4.5: Solve y′ + 6y = ex  .

Exercise 1.4.6: Solve y′ + 3x2y = sin (x)e-x3    , with y(0) = 1  .

Exercise 1.4.7: Solve y′ + cos(x)y = co s(x)  .

Exercise 1.4.8: Solve -12--y′ + xy = 3
x +1  , with y(0) = 0  .

Exercise 1.4.9: Suppose there are two lakes located on a stream. Clean water flows into the first lake, then the water from the first lake flows into the second lake, and then water from the second lake flows further downstream. The in and out flow from each lake is 500 liters per hour. The first lake contains 100 thousand liters of water and the second lake contains 200 thousand liters of water. A truck with 500 kg of toxic substance crashes into the first lake. Assume that the water is being continually mixed perfectly by the stream. a) Find the concentration of toxic substance as a function of time (in seconds) in both lakes. b) When will the concentration in the first lake be below 0.001 kg per liter. c) When will the concentration in the second lake be maximal.

Exercise 1.4.10: Newton’s law of cooling states that ddxt = - k(x - A )  where x is the temperature, t is time, A is the ambient temperature, and k > 0  is a constant. Suppose that A = A 0cosωt for some constants A
 0   and ω . That is, the ambient temperature oscillates (for example night and day temperatures). a) Find the general solution. b) In the long term, will the initial conditions make much of a difference? Why or why not.

Exercise 1.4.11: Initially 5 grams of salt are dissolved in 20 liters of water. Brine with concentration of salt 2 grams of salt per liter is added at a rate of 3 liters a minute. The tank is mixed well and is drained at 3 liters a minute. How long does the process have to continue until there are 20 grams of salt in the tank?

Exercise 1.4.12: Initially a tank contains 10 liters of pure water. Brine of unknown concentration of salt is flowing in at 1 liter per minute. The water is mixed well and drained at 1 liter per minute. In 20 minutes there are 15 grams of salt in the tank. What is the concentration of salt in the incoming brine?

1.5 Substitution

Note: 1 lecture, §1.6 in EP

Just like when solving integrals, one method to try is to change variables to end up with a simpler equation to solve.

1.5.1 Substitution

The equation

y′ = (x - y + 1)2.

is neither separable nor linear. What can we do? How about trying to change variables, so that in the new variables the equation is simpler. We will use another variable v , which we will treat as a function of x . Let us try

v = x - y + 1.

We need to figure out y′ in terms of v′ , v and x . We differentiate (in x ) to obtain v′ = 1 - y′ . So y′ = 1 - v′ . We plug this into the equation to get

     ′   2
1 - v = v .

In other words, v′ = 1 - v2   . Such an equation we know how to solve.

  1
-----2dv = dx.
1 - v

So

1   ||v + 1||
--ln ||-----||= x + C,
2   ||v - 1||
    ||v-+-1||   2x+2C
    ||v - 1|| = e    ,
or vv+-11 = De 2x  for some constant D . Note that v = 1  and v = -1  are also solutions.

Now we need to “unsubstitute” to obtain

x---y-+-2 = De2x,
  x - y

and also the two solutions x - y + 1 = 1  or y = x , and x - y + 1 = -1  or y = x + 2  . We solve the first equation for y .

    x - y + 2 = (x - y)De 2x,
                   2x      2x
    x - y + 2 = Dxe  - yDe  ,
  -y + yDe 2x = Dxe 2x - x - 2,
          2x       2x
y(- 1 + De ) = Dxe   - x - 2,
               Dxe 2x - x - 2
           y = -----2x------.
                  De  - 1
Note that D  = 0  gives y = x + 2  , but no value of D gives the solution y = x .

Substitution in differential equations is applied in much the same way that it is applied in calculus. You guess. Several different substitutions might work. There are some general things to look for. We summarize a few of these in a table.

When you seeTry substituting


  ′
yy  2
y
 2 ′
y y  3
y
       ′
(co sy)y sin y
(siny)y′ cosy
y′ey  ey

Usually you try to substitute in the “most complicated” part of the equation with the hopes of simplifying it. The above table is just a rule of thumb. You might have to modify your guesses. If a substitution does not work (it does not make the equation any simpler), try a different one.

1.5.2 Bernoulli equations

There are some forms of equations where there is a general rule for substitution which always works. One such example is the so-called Bernoulli equation† .

y′ + p(x)y = q(x)yn.

This equation looks a lot like a linear equation except for the yn  . If n = 0  or n = 1  then the equation is linear and we can solve it. Otherwise, a change of coordinates      1-n
v = y  transforms the Bernoulli equation into a linear equation. Note that n need not be an integer.

Example 1.5.1: Solve

xy′ + y(x + 1 ) + xy5 = 0, y(1) = 1.

First, the equation is Bernoulli (p (x) = (x + 1)∕x and q(x) = -1  ). We substitute

v = y1-5 = y-4,    v′ = - 4y-5y′.

In other words,  -1   5 ′   ′
( ∕4)y v = y . So

     ′              5
   xy  + y(x + 1) + xy = 0,
-xy5 ′              5
-4--v  + y(x + 1) + xy = 0,
  -x
  --v′ + y-4(x + 1) + x = 0,
  4
   --xv′ + v(x + 1) + x = 0,
    4
and finally
 ′  4(x + 1)
v - --------v = 4.
       x

Now the equation is linear. We can use the integrating factor method. In particular, we will use formula (1.4). Let us assume that x > 0  so |x| = x . This assumption is OK, as our initial condition is x = 1  . Let us compute the integrating factor. Here p(s)  from formula (1.4) is -4(s+1)
   s  .

               (∫ x            )                            -4x+4
  ∫x1 p(s)ds          -4(s-+-1)       -4x-4ln(x)+4    -4x+4 -4   e-----
 e       = exp   1      s    ds = e          = e     x  =   x 4 ,
 -∫xp(s)ds   4x+4ln(x)-4   4x-4 4
e  1     = e         = e    x .
We now plug in to (1.4)
               (                  )
        -∫xp(s)ds ∫  x ∫tp(s)ds
v(x) = e  1         e 1    4 dt + 1
             (∫   1            )
        4x-4 4    x e-4t+4
     = e   x      4  t4  dt + 1 .
                1

Note that the integral in this expression is not possible to find in closed form. As we said before, it is perfectly fine to have a definite integral in our solution. Now “unsubstitute”

            (  ∫ x  -4t+4     )
y-4 = e4x-4x4 4     e----dt + 1 ,
                1   t4
             e-x+1
  y = -(--∫-----------)--.
      x 4  x e-44t+4dt + 1 1∕4
          1   t

1.5.3 Homogeneous equations

Another type of equations we can solve by substitution are the so-called homogeneous equations. Suppose that we can write the differential equation as

       ( )
 ′      y-
y  = F  x .

Here we try the substitutions

    y
v = --    and therefore   y′ = v + xv′.
    x

We note that the equation is transformed into

                                                    ′
      ′                   ′                     ---v----   1-
v + xv = F(v)     or    xv =  F(v) - v    or    F (v) - v = x .

Hence an implicit solution is

∫      1
   --------dv = ln|x| + C.
   F (v) - v

Example 1.5.2: Solve

x2y′ = y2 + xy,    y(1 ) = 1.

We put the equation into the form y′ = (y∕x)2 + y∕x  . Now we substitute v = y∕x  to get the separable equation

  ′   2           2
xv = v + v - v = v ,

which has a solution

∫  1
   -2dv = ln|x| + C,
   v
    --1
     v  = ln|x| + C,
             -1
      v = ---------.
          ln|x| + C
We unsubstitute
y-= -----1---,
x   ln|x| + C
    -----x---
y = ln|x| + C .
We want y(1) = 1  , so
             - 1      -1
1 = y(1) =--------- = ---.
          ln |1| + C   C

Thus C  = -1  and the solution we are looking for is

y = ----x---.
    ln|x| - 1

1.5.4 Exercises

Exercise 1.5.1: Solve y′ + y(x2 - 1) + xy6 = 0  , with y(1) = 1  .

Exercise 1.5.2: Solve    ′       2
2yy + 1 = y + x , with y(0) = 1  .

Exercise 1.5.3: Solve y′ + xy = y4   , with y(0) = 1  .

Exercise 1.5.4: Solve          ∘ -------
yy′ + x =  x2 + y2   .

Exercise 1.5.5: Solve y′ = (x + y - 1)2   .

Exercise 1.5.6: Solve  ′   x2-y2
y =   xy  , with y(1) = 2  .

1.6 Autonomous equations

Note: 1 lecture, §2.2 in EP

Let us consider problems of the form

dx-
dt = f(x),

where the derivative of solutions depends only on x (the dependent variable). These types of equations are called autonomous equations. If we think of t as time, the naming comes from the fact that the equation is independent of time.

Let us come back to the cooling coffee problem. Newton’s law of cooling says that

dx-
dt = - k(x - A ),

where x is the temperature, t is time, k is some constant and A is the ambient temperature. See Figure 1.6 for an example with k = 0.3  and A = 5  .

Note the solution x = A (in the example x = 5  ). We call these types of solutions the equilibrium solutions. The points on the x axis where f (x ) = 0  are called critical points. The point x = A is a critical point. In fact, each critical point corresponds to an equilibrium solution. Note also, by looking at the graph, that the solution x = A is “stable” in that small perturbations in x do not lead to substantially different solutions as t grows. If we change the initial condition a little bit, then as t → ∞ we get x →  A . We call such critical points stable. In this simple example it turns out that all solutions in fact go to A as t → ∞ . If a critical point is not stable we would say it is unstable.


PIC

Figure 1.6: Slope field and some solutions of  ′
x  = -0.3(x - 5 )  .
PIC
Figure 1.7: Slope field and some solutions of  ′
x = 0.1x(5 - x)  .


Let us consider the logistic equation

dx-= kx(M  - x),
dt

for some positive k and M . This equation is commonly used to model population if you know the limiting population M , that is the maximum sustainable population. The logistic equation leads to less catastrophic predictions on world population than   ′
x  = kx . In the real world there is no such thing as negative population, but we will still consider negative x for the purposes of the math.

See Figure 1.7 for an example. Note two critical points, x = 0  and x = 5  . The critical point at x = 5  is stable. On the other hand the critical point at x = 0  is unstable.

It is not really necessary to find the exact solutions to talk about the long term behavior of the solutions. For example, from the above we can easily see that

          (
          ||| 5             if x(0) > 0,
          ||||
lim x(t) = |{ 0             if x(0) = 0,
t→∞       ||||| DNE  or - ∞   if x(0) < 0.
          |||(

Where DNE means “does not exist.” From just looking at the slope field we cannot quite decide what happens if x(0) < 0  . It could be that the solution does not exist t all the way to ∞ . Think of the equation x′ = x2   , we have seen that it only exists for some finite period of time. Same can happen here. In our example equation above it will actually turn out that the solution does not exist for all time, but to see that we would have to solve the equation. In any case, the solution does go to - ∞ , but it may get there rather quickly.

Often we are interested only in the long term behavior of the solution and we would be doing unnecessary work if we solved the equation exactly. It is easier to just look at the phase diagram or phase portrait, which is a simple way to visualize the behavior of autonomous equations. In this case there is one dependent variable x . So draw the x axis, mark all the critical points and then draw arrows in between. If f(x) > 0  draw an up arrow and if it is negative draw a down arrow.

PIC

Armed with the phase diagram, it is easy to sketch the solutions approximately.

Exercise 1.6.1: Try sketching a few solutions simply from looking at the phase diagram. Check with the preceeding graphs if you are getting the type of curves.

Once we draw the phase diagram, we can easily classify critical points as stable or unstable.

PIC

Since any mathematical model we cook up will only be an approximation to the real world, unstable points are generally bad news.

Let us think about the logistic equation with harvesting. Logistic equations are commonly used for modeling population. Suppose an alien race really likes to eat humans. They keep a planet with humans on it and harvest the humans at a rate of h million humans per year. Suppose x is the number of humans in millions on the planet and t is time in years. Let M be the limiting population when no harvesting is done. k > 0  is some constant depending on how fast humans multiply. Our equation becomes

dx-
dt =  kx(M  - x) - h.

Multiply out and solve for critical points

dx
--- = -kx2 + kMx - h.
 dt

Critical points A and B are

           ∘-----2------               ∘ ----2-------
A =  kM-+---(kM-)---4hk--    B = kM------(kM--)---4hk.
             2k                          2k

Exercise 1.6.2: Draw the phase diagram for different possibilities. Note that these possibilities are A > B , or A = B , or A and B both complex (i.e. no real solutions). Hint: Fix some simple k and M and then vary h .

For example, let M  = 8  and k = 0.1  . When h = 1  , then A and B are distinct and positive. The graph we will get is given in Figure 1.8. As long as the population stays above B , which is approximately 1.55 million, then the population will not die out. If ever the population drops below B , humans will die out, and the fast food restaurant serving them will go out of business.


PIC

Figure 1.8: Slope field and some solutions of x′ = 0.1x (8 - x) - 1  .
PIC
Figure 1.9: Slope field and some solutions of x′ = 0.1x(8 - x) - 1.6  .


When h = 1.6  , then A = B . There is only one critical point which is unstable. When the population is above 1.6 million it will tend towards this number. If it ever drops below 1.6 million, humans will die out on the planet. This scenario is not one that we (as the human fast food proprietor) want to be in. A small perturbation of the equilibrium state and we are out of business. There is no room for error. See Figure 1.9

Finally if we are harvesting at 2 million humans per year, the population will always plummet towards zero, no matter how well stocked the planet starts. See Figure 1.10.


PIC

Figure 1.10: Slope field and some solutions of  ′
x  = 0.1x (8 - x) - 2  .


1.6.1 Exercises

Exercise 1.6.3: Let  ′    2
x = x   . a) Draw the phase diagram, find the critical points and mark them stable or unstable. b) Sketch typical solutions of the equation. c) Find lim x(t)
t→ ∞  for the solution with the initial condition x (0) = -1  .

Exercise 1.6.4: Let x′ = sinx . a) Draw the phase diagram for - 4π ≤ x ≤ 4π . On this interval mark the critical points stable or unstable. b) Sketch typical solutions of the equation. c) Find limt→ ∞x (t)  for the solution with the initial condition x(0) = 1  .

Exercise 1.6.5: Suppose f(x)  is positive for 0 < x < 1  , it is zero when x = 0  and x = 1  , and it is negative for all other x . a) Draw the phase diagram for x′ = f (x)  , find the critical points and mark them stable or unstable. b) Sketch typical solutions of the equation. c) Find lim  x(t)
t→ ∞  for the solution with the initial condition x(0) = 0.5  .

Exercise 1.6.6: Start with the logistic equation dx
dt = kx(M - x)  . Suppose that we modify our harvesting. That is we will only harvest an amount proportional to current population. In other words we harvest hx per unit of time for some h > 0  (Similar to earlier example with h replaced with hx ). a) Construct the differential equation. b) Show that if kM >  h , then the equation is still logistic. c) What happens when kM  < h ?

1.7 Numerical methods: Euler’s method

Note: 1 lecture, §2.4 in EP

At this point it may be good to first try the Lab II and/or Project II from the IODE website: http://www.math.uiuc.edu/iode/.

As we said before, unless f (x, y)  is of a special form, it is generally very hard if not impossible to get a nice formula for the solution of the problem

y′ = f(x,y),    y(x 0) = y0.

What if we want to find the value of the solution at some particular x ? Or perhaps we want to produce a graph of the solution to inspect the behavior. In this section we will learn about the basics of numerical approximation of solutions.

The simplest method for approximating a solution is Euler’s method‡ . It works as follows: We take x0   and compute the slope k = f(x0,y0)  . The slope is the change in y per unit change in x . We follow the line for an interval of length h on the x axis. Hence if y = y0   at x0   , then we will say that y1   (the approximate value of y at x1 = x0 + h ) will be y1 = y0 + hk . Rinse, repeat! That is, compute x
 2   and y
 2   using x
 1   and y
 1   . For an example of the first two steps of the method see Figure 1.11.


PIC PIC

Figure 1.11: First two steps of Euler’s method with h = 1  for the equation y′ = y2
     3   with initial conditions y(0) = 1  .


More abstractly, for any i = 1, 2,3,... , we compute

xi+1 = xi + h,    yi+1 = yi + hf(xi,yi).

The line segments we get are an approximate graph of the solution. Generally it is not exactly the solution. See Figure 1.12 for the plot of the real solution and the approximation.


PIC

Figure 1.12: Two steps of Euler’s method (step size 1) and the exact solution for the equation  ′   y2
y  = 3   with initial conditions y(0) = 1  .


Let us see what happens with the equation      2
y′ = y∕3   , y(0) = 1  . Let us try to approximate y(2)  using Euler’s method. In Figures 1.11 and 1.12 we have essentially graphically approximated y(2)  with step size 1. With step size 1 we have y(2) ≈ 1.9 26  . The real answer is 3. So we are approximately 1.074 off. Let us halve the step size. Computing y
 4   with h = 0.5  , we find that y(2) ≈ 2.2 09  , so error of about 0.791. Table 1.1 gives the values computed for various parameters.

Exercise 1.7.1: Solve this equation exactly and show that y(2) = 3  .

The difference between the actual solution and the approximate solution we will call the error. We will usually talk about just the size of the error and we do not care much about its sign. The main point is, that we usually do not know the real solution, so we only have a vague understanding of the error. If we knew the error exactly …what is the point of doing the approximation?


h Approximate y(2)  Error PrevEiorruosrerror




1 1.92592592593 1.074074074070
0.5 2.20861152999 0.7913884700130.736809954840
0.25 2.47249414666 0.5275058533350.666557415634
0.125 2.68033658758 0.3196634124230.605990266083
0.0625 2.82040079550 0.1795992044970.561838476090
0.03125 2.90412106479 0.0958789352070.533849442573
0.015625 2.95035498158 0.0496450184220.517788587396
0.0078125 2.97472419486 0.0252758051420.509130743538

Table 1.1: Euler’s method approximation of y(2)  where of y′ = y2∕3   , y(0) = 1  .

We notice that except for the first few times, every time we halved the interval the error approximately halved. This halving of the error is a general feature of Euler’s method as it is a first order method. In the IODE Project II you are asked to implement a second order method. A second order method reduces the error to approximately one quarter every time you halve the interval.

Note that to get the error to be within 0.1 of the answer we had to already do 64 steps. To get it to within 0.01 we would have to halve another three or four times, meaning doing 512 to 1024 steps. That is quite a bit to do by hand. The improved Euler method from IODE Project II should quarter the error every time you halve the interval, so you would have to approximately do half as many “halvings” to get the same error. This reduction can be a big deal. With 10 halvings (starting at h = 1  ) you have 1024 steps, whereas with 5 halvings you only have to do 32 steps, assuming that the error was comparable to start with. A computer may not care about this difference for a problem this simple, but suppose each step would take a second to compute (the function may be substantially more difficult to compute than 2
y∕3   ). Then the difference is 32 seconds versus about 17 minutes. Note: We are not being altogether fair, a second order method would probably double the time to do each step. Even so, it is 1 minute versus 17 minutes. Next, suppose that you have to repeat such a calculation for different parameters a thousand times. You get the idea.

Note that we do not know the error! How do you know what is the right step size? Essentially you keep halving the interval and if you are lucky you can estimate the error from a few of these calculations and the assumption that the error goes down by a factor of one half each time (if you are using standard Euler).

Exercise 1.7.2: In the table above, suppose you do not know the error. Take the approximate values of the function in the last two lines, assume that the error goes down by a factor of 2. Can you estimate the error in the last time from this? Does it (approximately) agree with the table? Now do it for the first two rows. Does this agree with the table?

Let us talk a little bit more about the example  ′  y2
y =  ∕3   , y(0) = 1  . Suppose that instead of the value y(2)  we wish to find y(3)  . The results of this effort are listed in Table 1.2 for successive halvings of h . What is going on here? Well, you should solve the equation exactly and you will notice that the solution does not exist at x = 3  . In fact, the solution goes to infinity when you approach x = 3  .


h Approximate y(3)


1 3.16232281664
0.5 4.54328915766
0.25 6.86078752222
0.125 10.8032064113
0.0625 17.5989264104
0.03125 29.4600446195
0.015625 50.4012144477
0.0078125 87.7576927770

Table 1.2: Attempts to use Euler’s to approximate y(3)  where of      2
y′ = y∕3   , y(0) = 1  .

Another case when things can go bad is if the solution oscillates wildly near some point. Such an example is given in IODE Project II. In this case, the solution may exist at all points, but even a better approximation method than Euler would need an insanely small step size to compute the solution with reasonable precision. And computers might not be able to handle such a small step size anyway.

In real applications you would not use a simple method such as Euler’s. The simplest method that would probably be used in a real application is the standard Runge-Kutta method (we will not describe it here). That is a fourth order method, meaning that if you halve the interval, the error generally goes down by a factor of 16.

Choosing the right method to use and the right step size can be very tricky. There are several competing factors to consider.

You have seen just the beginnings of the challenges that appear in real applications. Numerical approximation of solutions to differential equations is an active research area for engineers and mathematicians. For example, the general purpose method used for the ODE solver in Matlab and Octave (as of this writing) is a method that appeared in the literature only in the 1980s.

1.7.1 Exercises

Exercise 1.7.3: Consider dx
---= (2t - x)2
dt   , x(0) = 2  . Use Euler’s method with step size h = 0.5  to approximate x(1)  .

Exercise 1.7.4: Consider dx
---= t - x
dt , x (0 ) = 1  . a) Use Euler’s method with step sizes h = 1,1∕2,1∕4,1∕8   to approximate x(1 )  . b) Solve the equation exactly. c) Describe what happens to the errors for each h you used. That is, find the factor by which the error changed each time you halved the interval.

Exercise 1.7.5: Approximate the value of e by looking at the initial value problem y′ = y with y(0) = 1  and approximating y(1)  using Euler’s method with a step size of 0.2  .