Calculus II:  Differential Equations Introduction

An Introduction to Differential Equations.

Modeling with Differential Equations

First of all-I do not expect you to understand all of the Mathematica that is used here. I don't expect to recreate all of this, although you should be able to understand it as we go through it. Some of the things you will need to learn (notably how to create vector fields) and we will continue to discuss those techniques in class in the coming weeks.

Clear[P,k,K,x,m]

Now, let's get started with our introduction to differential equations. We will base our introduction on understanding different models of population dynamics.

Exponential Growth

Let's start by looking at the exponential model of population growth. In this model, we assume that the rate of the population growth is proportional to the population:

exponential = P'[t] ==  k* P[t]

We may not know how to solve this differential equation (an equation that involves a derivative of a function), but we can get Mathematica to do it for us. Do you think it is wise to use Mathematica without understanding completely what you are doing mathematically?

solexponential = DSolve[exponential,P[t],t]
P[t_] = P[t]/.First[solexponential]

The solution is a family of curves. C[1] is the constant of integration, and as we vary C[1] we will get different curves. Remember that k is simply the proportionality constant, so we will set it equal to 1 here. The way to do a substitution in Mathematica is with the /. command.

functionlist=P[t]/. {C[1]-> Range[-10,10,1] ,k->1}

familyexp = Plot[Evaluate[functionlist],{t,-2,2},
PlotRange -> {-10,10},AxesLabel-> {"t","Solutions to DE"},PlotStyle->{Red}]

What have we done? We have solved our mathematical model, but we need to interpret the result. The result has to make sense relative to the physical system it is modeling. We are modeling population growth, so what restrictions should we impose? I think that the time should be greater than or equal to zero, and the population should be greater than or equal to zero also. Therefore, in terms of the underlying model, the only part of the graph we are interested in is:

lines = Plot[Evaluate[functionlist],{t,0,2},
PlotRange -> {0,10},AxesLabel-> {"t","Solutions to DE"},PlotStyle->{Red}]

If we are told what an initial condition is for the system, we can pick out the particular curve which models the system. If we are told that initially the population is 3, then we need to force our model to satisy P[0]=3. This determines the constant of integration, and we get a single solution:

IC = Solve[P[0]== 3,C[1]]
P[t_] = P[t]/.First[IC]

single = Plot[P[t] /.k->1,{t,0,2},
PlotRange -> {0,10},AxesLabel-> {"t","Solution to IVP"},PlotStyle->{Blue,Thick}]

Mathematica can solve an initial value problem (a problem with a differential equation and an initial condition) in one command:

Clear[P]
DSolve[{P'[t]==k P[t],P[0]==3},P[t],t]

Solving an Initial Value Problem by Hand (Separable Equations)

We can do everything that Mathematica just did for us ourselves. Obviously, that is important, since we never want to use Mathematica to do our thinking for us! So here is an example of how to solve a particular initial value problem by hand. This problem falls into a class of problems which are known as separable, and this is how all problems in that class are solved.

P'(t) = k P(t), P(0) = (Subscript[P, 0])

We begin by solving the differential equation part of the initial value problem.

((d P)/(d t))= k P

Separate, which means put all the P on one side of the equation, all the t on the other:

((d P)/P) = k d t

Now, integrate both sides of the equation:

∫((d P)/P) = ∫k d t 
ln |P| = k t + (Subscript[C, 1])
|P| = (e^(k t + (Subscript[C, 1])))= Exp( k t +(Subscript[C, 1]))
|P| = (e^(k t ))(e^(Subscript[C, 1]))
|P|=(Subscript[C, 2]) (e^(k t )), where (Subscript[C, 2]) = (e^(Subscript[C, 1]))>0.

Now, there some subtle manipulations here, but as always if we let the math lead us we arrive at the correct solution. First, our constant (Subscript[C, 2]) is greater than zero, and that makes sense since we have the absolute value of P. We want P, so we write

P(t)="±"(Subscript[C, 2]) (e^(k t )), where (Subscript[C, 2]) = (e^(Subscript[C, 1]))>0.
P(t) = C (e^(k t )), where C="±"(Subscript[C, 2]) ∈(-∞,∞).

Now, we use the initial condition to determine the constant of integration. Note that you don't want to just say that C=(Subscript[P, 0])! There is a reason why this is true in this case, and it doesn't happen in general.

P(0) = C (e^(k (0) ))= C 
P(0) = (Subscript[P, 0])
C = (Subscript[P, 0])

And the solution to the initial value problem is

P(t)=(Subscript[P, 0]) (e^(k t )).

We can check that this is a solution by taking a derivative and verifying the differential equation is satisfied.

P'(t)=k (Subscript[P, 0]) (e^(k t ))= k P(t), and P(0)=(Subscript[P, 0]).

Problems with the Exponential Model of Population Growth

Plot P[t] for a large time period, where (Subscript[P, 0])=3, and the growth rate k=1.

Plot[3Exp[1*t],{t,0,20},AxesLabel-> {"t","Population"},PlotStyle->{Blue}]

Does P[t] seem to be a realistic view of population growth? I don't think it is, since the population grows forever, and that isn't possible, due to a lack of resources, space, and other factors. This leads us to look at models which more accurately reflect the true nature of population growth.

The Logistic Equation: An Improved Population Model

A better model of population growth which does not have the problem of runaway growth is the logistic equation. The logistic equation will not have population grow to infinity, yet it will have exponential growth for small populations. We will be talking more about it later in Sec 9.5.

Logistic = P'[t] == k (1 - P[t]/K) P[t]

Let's solve this equation and see what we get:

solLogistic = DSolve[Logistic,P[t],t]
P[t_] = P[t]/.First[solLogistic]

Let's plot the family of solutions for K=7 and k=0.5:

functionlist=P[t]/. {C[1]-> Range[-10/7,10/7,1/7] ,k->0.5,K->7}

familylog = Plot[Evaluate[functionlist],{t,0,2},
PlotRange -> {0,10},AxesLabel-> {"t","Solutions to DE"},PlotStyle->{Red}]

HEY! Where did the solutions for y less than 7 go? The answer is somewhat tricky, and one of the reasons you need to understand what you are doing before you rush off into Mathematica land.

The equation in the given form is in terms of the constant C[1]. We would much rather work with the initial value P(0) = (Subscript[P, 0])=P0, so we will make the following substitution:

P[t] = P[t] /. K C[1] -> Log[P0/(P0-K)]

Why did we make that particular substitution? Here's why:

P(0)=(((e^(k (0)+K C[1])) K)/(-1+(e^(k (0)+K C[1]))))=(Subscript[P, 0]) (solve for K C[1])
(((e^(K C[1])) K)/(-1+(e^(K C[1]))))=(Subscript[P, 0])
(e^(K C[1])) K=(Subscript[P, 0])(-1+(e^(K C[1])))
(e^(K C[1]))( K-(Subscript[P, 0]))=-(Subscript[P, 0])
(e^(K C[1]))=((-(Subscript[P, 0]))/( K-(Subscript[P, 0])))
(e^(K C[1]))=((Subscript[P, 0])/( (Subscript[P, 0])-K))
K C[1]=ln(((Subscript[P, 0])/( (Subscript[P, 0])-K)))

You should compare the previous two representations of the solution. In the first,  (e^(C[1])) can only be positive, and we cannot generate all solutions of the differential equation for all initial conditions. Specifically, 0 <P(0) < K would be missing if we just assigned values to C[1]. This is a way that our knowledge of differential equations allows us to interpret the results that Mathematica gives us. And this is a very important interpretation indeed!

functionlist=P[t]/. {P0-> Range[0.1,15,1] ,k->0.5,K->7}

familylog = Plot[Evaluate[functionlist],{t,0,10},
PlotRange -> {0,10},AxesLabel-> {"t","Solutions to DE"},PlotStyle->{Red}]

Can you identify the equilibrium solutions? The equilibrium solutions are where the solution is a constant for all time. Can you prove what the equilibrium solutions are from the differential equation?

The equilibrium solutions will not change (hence the name), so their rate of change with time will be zero. That means that

((d P)/(d t))= k (1 - (P/K))P = 0

So the equilibrium solutions are P=0, P=K. These are, necessarily, constant value solutions.

Direction Fields: Getting a Picture of the Familyof Curves Without Solving the Differential Equation

Clear[P,k,K,x,m,y]

The importance of direction fields in the solution of differential equations can not be understated. The direction field is found from plotting short arrows which represent the slope of the solution at that point in space. We know the slope, since for a first order differential equation like we are looking at, the slope is given by the derivative which we know:

((d y)/(d t))=f(t,y)

For example, in exponential growth, the differential equation is:

exponential

and we can plot P as a function of t.  The vector field command is given below (remember we loaded this package right at the start!), although you can see that it doesn't look like much.

fieldexp=VectorPlot[{1,Last[exponential] /. k-> 1},{t,-2,2},{P[t],-10,10}, ImageSize->{300,300},AspectRatio->1,VectorScale->{Tiny,Tiny,None}]

The option I highlighted in red makes the vector field much easier to interpret.

If we put our family of solutions against the vector field, the meaning of a direction field becomes clear:

Show[familyexp,fieldexp]

Note the command

Last[exponential]

just allows me to quickly get the right hand side of the differential equation. For the differential equation dy/dx=sin(x)+(y^2)we can get the direction field using

fieldexp=VectorPlot[{1,Sin[x]+y^2},{x,-2,2},{y,-2,2}, ImageSize->{300,300},AspectRatio->1,VectorScale->{Tiny,Small,None}]

Direction fields often look best when the region in the x direction is the same length as in the y direction.

Another Direction Field  Example

Let's look at the direction fields and family of solutions for the differential equation given by

eqn = y'[x] == x+y[x]

The direction field looks like:

fieldeqn=VectorPlot[{1,Last[eqn]},{x,-2,2},{y[x],-5,5}, ImageSize->{300,300},AspectRatio->1,VectorScale->{Tiny,Tiny,None}]

The solutions are given by

soleqn = DSolve[eqn,y[x],x]

y[x_] = y[x]/.First[soleqn]

functionlist=y[x]/. C[1]-> Range[-10,10,1]

family = Plot[Evaluate[functionlist],{x,-2,2},
PlotRange -> {-5,5}, ImageSize -> {300,200},AxesLabel-> {"x","Solutions to DE"},PlotStyle->{Red}]

Show[family,fieldeqn]

What is the equilibrium solution?

((d y)/(d x))=0=x+y
y(x) = - x ≠ constant.

There is no equilibrium solution. Remember, the equilibrium solutions must be constant!

Extending to Higher Order Derivatives: Motion of a Spring

A second order (two derivatives, not one) differential equation will need two integrations to solve, and so will have two constants of integrations. An initial value problem for a second order differential equation will require two initial conditions to determine the two constant of integration.

Spring = m x''[t] == -k x[t]

Let's solve this equation and try see what we get:

solSpring = DSolve[{Spring,x[0]==10,x'[0]== 3},x[t],t]
x[t_] = x[t]/.First[solSpring]

springline = Plot[x[t] /. k ->2 /. m->3,{t,0,16},
PlotRange -> {-20,20},AxesLabel-> {"t","Solution to IVP"}]

Problem

Consider the differential equation:

Clear[y,x,eqn]
eqn= y'[x] == Cos[x]^2

Can you use Mathematica to create the direction field for the differential equation?
Can you use Mathematica to find the family of curves that satisfies this differential equation, and plot them? Can you overlay it on the direction field?
Can you solve the differential equation by hand using the separation of variables technique?

Try it yourself before looking at my solution!

Direction Field

fieldeqn=VectorPlot[{1,Cos[x]^2},{x,-2,2},{y,-2,2}, ImageSize->{300,300},AspectRatio->1,VectorScale->{Tiny,Small,None}]

DSolve Solution and Plots

soleqn = DSolve[y'[x] == Cos[x]^2,y[x],x]

y[x_] = y[x]/.First[soleqn]

functionlist=y[x]/. C[1]-> Range[-10,10,1]

family = Plot[Evaluate[functionlist],{x,-2,2},
PlotRange -> {-2,2}, ImageSize -> {300,200},AxesLabel-> {"x","Solutions to DE"},PlotStyle->{Red}]

Show[family,fieldeqn]

Hand solution

Here is the solution by hand (notice that integration techniques will be important here!)

y'[x]==((Cos[x])^2)
dy==((Cos[x])^2) dx
∫dy==∫((Cos[x])^2) dx
y==∫(1/2)(1+Cos[2x]) dx
y==(1/2)x+(1/2)∫Cos[2x] dx        let u=2x, du=2 dx
y==(1/2)x+(1/4)∫Cos[u] du       
y==(1/2)x+(1/4)Sin[u]+ C
y==(1/2)x+(1/4)Sin[2x]+ C