MATH 302-010
Ordinary Differential Equations
Prof. D. A. Edwards
Copyright 2024

Introduction to ODE Tools (10/1 Version)

This command makes sure that Mathematica is starting with all variables cleared:

Quit[]

If you look after the title of the section headings below, you will see an arrow.  This indicates that the section has been collapsed.  Click on it to expand the section.

Direction Fields

To plot the direction field for an ODE, we use VectorPlot[{1,y’},trange,yrange].  This will plot a series of vectors whose slope is y', and whose color changes depending on the size of |y'|.  For instance, for the equation y'=2sin(3y), we have the following direction field:

rhs1 = 2*Sin[3*y]
fieldplot = VectorPlot[{1,%},{t,0,1},{y,0,2}]

Recall that we can enter multiple definitions on separate lines by hitting <RETURN>; only when we hit <ENTER> (or <SHIFT>-<RETURN>) are the commands evaluated.  Also recall that we can use the % operator to refer to the line immediately before.

The direction field suggests convergence to a few equilibrium solutions, which can be found by setting the right-hand side of the equation equal to zero.

Solve[rhs1==0,y]

When setting up an equation to solve, we use two equals signs (==), since one equals sign means a definition.  In this case, we get an infinite list of solutions.  To find the ones on the graph, we restrict the domain of our search and use the command FindRoot[equation,{variable,guess}] instead.

equil1 = FindRoot[rhs1==0,{y,0}]
equil2 = FindRoot[rhs1==0,{y,1}]

Note that FindRoot returns a “rule”.  To do our plots, we only need the right-hand side, which we can find by substituting the rule into y.  We use /. to do the substitution:

equil1a = y/.equil1
equil2a =  y/.equil2

Then we just plot the two values together (in a list, so in braces) to see the integral curves:

eqplot = Plot[{equil1a,equil2a},{t,0,1},PlotStyle->{Red,Red}]

Note that we use the PlotStyle command to select the line colors. Then we use the command Show command to make everything display together:

Show[fieldplot,eqplot]

Exercise

1.  Graph the direction field and the equilibrium solution for the differential equation y'(t)+3y(t)=1.

Solution

rhs2=1-3*y
Solve[rhs2==0,y]
equil2 = y/.%

field2 = VectorPlot[{1,rhs2},{t,0,1},{y,0,2}]
eq2plot = Plot[equil2,{t,0,1},PlotStyle->Red]
Show[field2,eq2plot]

Basic ODE Tools

We begin by defining a differential equation.  The easiest way is just to use the prime notation:

correcteq=y'[t]+7*y[t]==3

When setting up an equation to solve, we use two equals signs (==), since one equals sign means a definition.

Once we have a differential equation, we use the command DSolve to solve it, if possible.  The basic syntax (no initial condition) is DSolve[equation,function,independent variable]:

DSolve[correcteq,y[t],t]

Here the (Subscript[c, 1])(or C[1], depending on the computer) represents the arbitrary constant, which we know occurs if no initial condition is specified.

If we have an initial condition, we define that as a separate equation, and then use DSolve.  The syntax in this case is DSolve[{equation,IC},function,independent variable]:

initcon=y[0]==5;
f1=DSolve[{correcteq,initcon},y[t],t]

Note that we make a set of equations by enclosing them in braces.  

WARNING: It is critical that you use two equals signs when defining the initial conditions.  If instead you accidentally type something like y[0]=5, you’ve now defined a variable, instead of set up an equation.  This is exceedingly hard to undo; the easiest way I’ve found is to restart the computational engine with the following command:

Quit[]

Then you’ll need to re-evaluate any of the previous commands.

Note that in each of the examples above, DSolve returns a rule y[t]->....  For plotting, it is sometimes easier just to have the right-hand side, which we may obtain by using DSolveValue instead of DSolve:

plotf1=DSolveValue[{correcteq,initcon},y[t],t]

One thing that we have to be careful about when solving nonlinear differential equations is to make sure the exponent is placed in the correct place.  For instance consider the following example:

correcteq2=y'[t]+t*(y[t])^2==0

Here we have placed the entire function y(t) in parentheses, so the whole function is squared.  Mathematica handles the correct equation just fine:

DSolve[correcteq2,y[t],t]

However, there are other cases where it won’t give you all the solutions (see the homework), so you have to be careful.

If we put the exponent just on the y, Mathematica doesn’t handle it properly:

wrongeq=y'[t]+t*y^2[t]==0
DSolve[wrongeq,y[t],t]

Exercises

1.  Solve the differential equation y'(t)+3y(t)=1, both in general and with the initial condition y(0)=1.

Solution

eq1=y'[t]+3*y[t]==1
ic1=y[0]==1

DSolve[eq1,y[t],t]
DSolve[{eq1,ic1},y[t],t]

2.  Solve the differential equation y'(t)=(t/(sin(y))), both in general and with the initial condition y(0)=0.

Solution

eq2=y'[t]==t/Sin[y[t]]
ic2=y[0]==0

DSolve[eq2,y[t],t]
DSolve[{eq2,ic2},y[t],t]

Integral Curves

To plot one solution, we use the Plot command.  If we have a rule we obtained using DSolve, we evaluate y(t) (with the /. command) using the solution to the ODE:

Plot[y[t]/.f1,{t,0,1}]

If we have the right-hand side from DSolveValue, we just plot it directly:

Plot[plotf1,{t,0,1}]

Now suppose we want to plot a series of integral curves of the solution to correcteq for various values of y(0).  First, we set up the solution:

genic=y[0]==y0
f2=DSolve[{correcteq,genic},y[t],t]

If we want to plot a series of integral curves together, we have to create a table of the values.  (This is equivalent to creating a looping structure.)  The syntax is Table[expression,{variable to substitute,lower limit,upper limit}].  Therefore, this generates a table of solutions with (Subscript[y, 0]) running from 1 to 5, incrementing by 1:

f2tab = Table[y[t]/.f2,{y0,1,5}]

Then we just plot the table to see the integral curves:

Plot[f2tab,{t,0,1},PlotRange->All]

Here we use the PlotRange->All option because otherwise the graph cuts off.  Note the initial conditions are {1,2,3,4,5}, as expected.

We can also use the same technique to plot the solutions to different equations.

parameq = {y'[t]+a*y[t]==Exp[-3*t],y[0]==-1}

So now we have an equation that depends on a parameter.  But again we can evaluate the solution for various vales of the parameter (in this case, for a={-2,-1,0,1,2}).

solset=Table[DSolve[parameq,y[t],t],{a,-2,2}]
Plot[y[t]/.solset,{t,0,1}]

Exercise

1.  Graph five integral curves for the differential equation y'(t)+3y(t)=1 with the initial conditions y(0)=-1, 0, 1, 2, and 3.

Solution

eq1=y'[t]+3*y[t]==1
genic=y[0]==y0

DSolve[{eq1,genic},y[t],t]
Table[y[t]/.%,{y0,-1,3}]
Plot[%,{t,0,1},PlotRange->All]

Advanced ODE Tools

Sometimes there are first-order equations for which a closed-form solution can’t be written in a nice form:

numeq=y'[t]+Sin[t^3]*Exp[-t]*y[t]==0
DSolve[numeq,y[t],t]

For these types of equations, we use NDSolve instead.  This will force Mathematica to solve the problem numerically.  The syntax in this case is NDSolve[{equation,IC},function,{independent variable,lower limit,upper limit}].  IMPORTANT: Note that you do NOT put the independent variable with the function.

numic = y[0]==1
numsolve = NDSolve[{numeq,numic},y,{t,0,3}]

Note that the output isn’t helpful, so we usually suppress it with a semicolon.  However, with the expression now defined (as a rule), we can compute y at various points using evaluation.

y[3]/.numsolve

Plot[y[t]/.numsolve,{t,0,3}]

Note that y(3) in the graph is around 0.83, as expected.

The reason why we have to put a list of t-values is so that Mathematica won’t try to compute the solution numerically forever.  We have to make sure we use the right values, because going outside the range won’t work:

y[4]/.numsolve
Plot[y[t]/.numsolve,{t,0,4}]

To do integral curves, we again create a table.  Here the syntax is slightly different.  We do the table first:

genic=y[0]==y0
num2=Table[NDSolve[{numeq,genic},y,{t,0,3}],{y0,-2,2}]

Then we do the evaluation at the plot step:

Plot[Evaluate[y[t]/.num2],{t,0,3}]

Here the extra Evaluate forces the curves to be different colors.

Exercise

1.  Consider the differential equation y'(t)+cos((t^2))(e^(-t))y(t)=sin(t). Solve the problem for y(0)=1 and find y(3).  Then plot integral curves for   y(0)=0, 1, 2, 3, and 4.

Solution

numeq1=y'[t]+Cos[t^2]*Exp[-t]*y[t]==Sin[t]
numic1=y[0]==1
numsol = NDSolve[{numeq1,numic1},y,{t,0,3}];
y[3]/.numsol

genic1=y[0]==y0
icur1=Table[NDSolve[{numeq1,genic1},y,{t,0,3}],{y0,0,4}];
Plot[y[t]/.icur1,{t,0,3}]

Inserting Comments

Usually when creating a Mathematica notebook that you want to show to someone else (e.g., the grader), it’s wise to add comments to explain what you are doing.  That way others can follow your work and track down any possible errors.  To insert text, from a blank prompt hit CMD-7 (on a Mac), ALT-7 (on a PC), or Select Format..Style..Text.  Try this below:

I expect you to insert comments in any of your Mathematica printouts so that the grader may follow your reasoning.