Quick Review of Differential Equations from Calculus II

Barry McQuarrie
Aug, 2014

Let's get started with our introduction to differential equations. This basically reviews the concepts seen in calculus, and introduces you to some of the Mathematica syntax that is commonly used. We will not be using a great deal of Mathematica in this class, but it will show up on assignments. The topics we touch on today will be discussed again in the coming weeks.

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:

L[P_]:=  D[P,t]==  k* P

This differential equation (an equation that involves a derivative of a function) is separable, and therefore easy to solve by hand. The differential equation from above is

L[P[t]]

We can also get Mathematica to solve it for us. Do you think it is wise to use Mathematica without understanding completely what you are doing mathematically?

sol = DSolve[L[P[t]],P[t],t]
P[t_] = P[t]/.First[sol]

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 (not a constant of integration), 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},PlotStyle->Red,AxesLabel-> {"t","General Solutions"},ImageSize->{500,Automatic}]

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:

general = Plot[Evaluate[functionlist],{t,0,2},
PlotRange -> {0,10},PlotStyle->Red,AxesLabel-> {"t","General Solutions"},ImageSize->{500,Automatic}]

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},PlotStyle->{Thick,Blue},AxesLabel-> {"t","Particular Solution"},ImageSize->{500,Automatic}]

Show[general,single]

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

Clear[P]
DSolve[{L[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 by hand. 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]).

Motion of a Spring -- Second Derivative in Equation

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.

Clear[x]

L[x_]:=  m D[x,{t,2}]==-k x

L[x[t]]

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

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

Here is a plot for specific values of k and m.

springline = Plot[x[t] /. k ->2 /. m->3,{t,0,16},
PlotRange -> {-20,20},AxesLabel-> {"t","Particular Solution"},PlotStyle->{Red,Thick},ImageSize->{400,Automatic}]

We will spend considerable time in this course understanding solutions to differential equations like this one (second order linear).

We can solve a differential equation, and then examine the solution graphically subject to varying some parameters using Manipulate:

Clear[x]

sol = DSolve[{L[x[t]],x[0]==10,x'[0]== 3},x[t],t]

Manipulate[Plot[((10 (Sqrt[k]) Cos[(((Sqrt[k]) t)/(Sqrt[m]))]+3 (Sqrt[m]) Sin[(((Sqrt[k]) t)/(Sqrt[m]))])/(Sqrt[k])),{t,0,16},PlotRange->{{0,16},{-20,20}},PlotStyle->Thick,AxesLabel->{"t","Particular Solution"}],{m,3,10},{k,2,10}]

Direction Fields

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

The importance of direction fields in the solution of first order 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:

L[P_]:=  D[P,t]==  k* P

and the function f(t,y) can be obtained using the Mathematica command Last.

Last[L[P]]

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[Evaluate[{1,Last[L[P]] /. k-> 1}],{t,-2,2},{P,-10,10},VectorScale->{Tiny,Tiny,None}]

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

Show[familyexp,fieldexp]

Another Example

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

L[y_]:=  D[y,x] == x+y

The direction field looks like:

fieldeqn=VectorPlot[Evaluate[{1,Last[L[y]]}],{x,-2,2},{y,-5,5}, ImageSize->{400,300},VectorScale->{Tiny,Small,None}]

The solutions are given by

soleqn = DSolve[L[y[x]],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},AxesLabel-> {"x","General Solutions"}]

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!

Euler's Method: Recursion Relations

Solve IVP:  (dy/dt)=f(t,y), y((Subscript[t, 0]))=(Subscript[y, 0]) in the interval t∈((Subscript[t, 0]),b).

Euler's Method is a numerical technique to solve initial value problems. You probably saw it in Calculus II, and it is discussed in our text in Section 2.7. We will improve on Euler's Method by deriving Taylor's Method of Order n, which is Euler's method when n=1. Using n=2 will give us much more accurrate results.

y[12]

Clear[f,y,t,"(Subscript[t, 0])","(Subscript[y, 0])",h,b,data]

f[t_,y_]=y-(t^2)*Cos[3y];

(Subscript[t, 0])=0.0;
(Subscript[y, 0])=0.72;
b=4.0;

num=100;
h=((b-(Subscript[t, 0]))/num);

t[0]=(Subscript[t, 0]);
y[0]=(Subscript[y, 0]);
t[n_]:=t[n]=(Subscript[t, 0])+n*h
y[n_]:=y[n]=y[n-1]+h*f[t[n-1],y[n-1]]

data=Table[{t[n],y[n]},{n,0,num}];

euler=ListPlot[data,PlotRange->{{0,5},{0,5}},Joined->True,Frame->True,FrameLabel->{"t","y"},PlotStyle->{Red,Thick}]

field=VectorPlot[{1,y-(t^2)Cos[3y]},{t,0,4},{y,0,5},VectorScale->{Tiny,Small,None}];
Show[field,euler,ImageSize->{500,500}]

With numerical techniques, be sure to investigate changing the step size (h), since you can easily get crappy results if you choose too large a step size.

Clear[f,y,t,"(Subscript[t, 0])","(Subscript[y, 0])",h,b,data]

f[t_,y_]=y-(t^2)*Cos[3y];

(Subscript[t, 0])=0.0;
(Subscript[y, 0])=0.72;
b=4.0;

num=50;
h=((b-(Subscript[t, 0]))/num);

t[0]=(Subscript[t, 0]);
y[0]=(Subscript[y, 0]);
t[n_]:=t[n]=(Subscript[t, 0])+n*h
y[n_]:=y[n]=y[n-1]+h*f[t[n-1],y[n-1]]

data=Table[{t[n],y[n]},{n,0,num}];

euler2=ListPlot[data,PlotRange->{{0,5},{0,5}},Joined->True,Frame->True,FrameLabel->{"t","y"},PlotStyle->{Blue,Thick}];

Show[field,euler,euler2,ImageSize->{500,500}]