Chapter 3 Block 3: Dynamical Systems
3.1 USAFACalc for Block 6
In Block 6 we follow the Supplemental Reading available online here. All of the R functionality that we need for this block has now been built into USAFACalc.
Below find the R commands we’ll be using in class. These all have reasonable documentation, you can find documentation and examples by using the ? operator, for example ?plotODEDirectionField.
Students should be able to use plotODEDirectionField, Euler, and plotPhasePlane for the homework, GR, and/or final exam. They’ll use animateODESolution for the project, but that’s it.
3.1.1 plotODEDirectionField
plotODEDirectionField Plots the direction field for a scalar ODE. Direction field, not slope field (direction fields have normalized vectors). There is no command for plotting slope fields, although you could certainly use plotVectorField to achieve it by defining the vector-valued function \((1,f(t,y))\).
In general, the R ODE functions only require the right-hand side of first-order ODEs. So to plot a direction field for the ODE
\[y'=t-y\]
use the following, note that we are only giving it the right-hand side. In general the right-hand side ODEs need to be functions of \(t\) and \(y\), in that order. Don’t reverse \(t\) and \(y\), and don’t omit either of them, even if one of them doesn’t show up in the ODE.
You can vary the arrow thickness and color using lwd and col, as usual. You can plot trajectories that correspond to initial conditions as follows. You can plot one or many trajectories. It is assumed that the initial condition always occurs on the left-hand side of the figure, so you only need to give the \(y\)-value of the initial condition(s).
plotODEDirectionField(t-y~t&y,tlim=c(0,10),ylim=c(-5,5),ics=2)
plotODEDirectionField(t-y~t&y,tlim=c(0,10),ylim=c(-5,5),ics=c(2,-4))
plotODEDirectionField works fine with the add option, or other things can be added to a direction field, so you can overlay plots of whatever function you like on top of the direction field. For instance a solution to the ODE, or a not-the-solution to the ODE.
Consider the ODE
\[y'=y(y-4).\]
Note how there is no \(t\)-dependence on the right-hand side of the ODE. With plotODEDirectionField you can omit the \(t\) variable. You can not omit the \(y\) variable if the ODE is \(y'=t(t-4)\). It is safest to always declare the right-hand side to be a function of \(t\) and \(y\).
Sometimes your direction field will look weird because the aspect ratio is not 1:1 by default. If your vertical arrows look short, but horizontal ones look long, the aspect ratio your problem. You can set asp=1 to fix this.
3.1.2 Euler
Euler performs Euler’s method on a given first-order ODE. The goal is that students should be able to do Euler’s method by hand, but when it comes time to do it for lots of iterations then obviously we’ll just use the computer.
The syntax is very similar to plotODEDirectionField. Specify the right-hand side of the ODE, give the time interval you are interested in and the value of the initial condition, and optionally set the step size. You do not need to stick to \(y\) as the unknown function, you can name it what you like.
## t P
## 1 0.0 0.1000000
## 2 0.1 0.1490000
## 3 0.2 0.2212799
## 4 0.3 0.3270234
## 5 0.4 0.4798406
## 6 0.5 0.6967362
## 7 0.6 0.9965602
## 8 0.7 1.3955271
## 9 0.8 1.8985411
## 10 0.9 2.4873658
## 11 1.0 3.11234983.1.2.1 Fancier Euler
The command Euler has had some extra work put into it. It will let you specify the ODE as an equality, in standard form, so long as you use “_t” to denote derivative. For example, the IVP
\[P'=P(5-P),\quad P(0)=0.1\]
could be approximated with
## t P
## 1 0.0 0.1000000
## 2 0.1 0.1490000
## 3 0.2 0.2212799
## 4 0.3 0.3270234
## 5 0.4 0.4798406
## 6 0.5 0.6967362
## 7 0.6 0.9965602
## 8 0.7 1.3955271
## 9 0.8 1.8985411
## 10 0.9 2.4873658
## 11 1.0 3.1123498Allowing the specification of the entire equation will help with clarity later on.
3.1.3 Euler’s Method for Systems of ODE
Our Euler’s method can approximately solve arbitrarily large systems of ODEs. To solve the system of equations
\[\begin{aligned}
&x'=-y&\quad x(0)=1\\
&y'=x&\quad y(0)=0.5\\
\end{aligned}\]
We can use a syntax very much like findZeros, where we just specify the right-hand side of the ODE, specified as a vector-valued function.
## t x y
## 1 0.0 1.0000000 0.5000000
## 2 0.5 0.7500000 1.0000000
## 3 1.0 0.2500000 1.3750000
## 4 1.5 -0.4375000 1.5000000
## 5 2.0 -1.1875000 1.2812500
## 6 2.5 -1.8281250 0.6875000
## 7 3.0 -2.1718750 -0.2265625
## 8 3.5 -2.0585938 -1.3125000
## 9 4.0 -1.4023438 -2.3417969
## 10 4.5 -0.2314453 -3.0429688
## 11 5.0 1.2900391 -3.1586914This works, however, the order of the variable declaration is very important. It is assumed that because \(x\) is the first variable in t&x&y, that the first function in c(-y,x) is the derivative of \(x\), and that the first value in c(1,0.5) is the initial value of \(x\). This is easy to mess up.
It is better to use the option to specify the entire equation with Euler in this case.
## t x y
## 1 0.0 1.0000000 0.5000000
## 2 0.5 0.7500000 1.0000000
## 3 1.0 0.2500000 1.3750000
## 4 1.5 -0.4375000 1.5000000
## 5 2.0 -1.1875000 1.2812500
## 6 2.5 -1.8281250 0.6875000
## 7 3.0 -2.1718750 -0.2265625
## 8 3.5 -2.0585938 -1.3125000
## 9 4.0 -1.4023438 -2.3417969
## 10 4.5 -0.2314453 -3.0429688
## 11 5.0 1.2900391 -3.1586914Now everything is unambiguous, you can mess up the order of things as much as you like.
3.1.4 plotPhasePlane
For systems of two autonomous ODEs we have plotPhasePlane, which works very much like plotODEDirectionField.
For the system of ODE \[\begin{aligned} &x'=-y\\ &y'=x\\ \end{aligned}\] we can plot the phase plane as below.
plotPhasePlane also allows you to add trajectories to your phase plane. If you wish to plot the trajectory associated with the initial condition \(x(0)=0\), \(y(0)=-2\), use the code
If you want to specify multiple initial conditions, you must put them in a list.
plotPhasePlane works fine with the usual add argument, so you can add a plot of an arbitrary function or whatever you like on this. However, remember that for a phase plane you most likely want a parameteric equation. For example, one solution to the system above is
\[x(t)=\cos(t)\]
\[y(t)=\sin(t).\]
To plot the trajectory on the phase portrait, you’ll need to evaluate \(x\) and \(y\) at a large number of \(t\)-values, and then plot the \(x-y\) pairs. To do that, you’ll want to use plotPoints, likely with the optional argument type="l" to make it connect your points with a line. Example:
#Just make the phase plane, no trajectory.
plotPhasePlane(c(-y,x)~x&y)
#Define my x(t), y(t)
x=makeFun(cos(t)~t)
y=makeFun(sin(t)~t)
ts=seq(0,2*pi,by=0.1)
plotPoints(y(ts)~x(ts),add=TRUE,type="l",col="firebrick",lwd=3)
3.1.5 animateODESolution
Phase planes with trajectories on them are great, but it is hard to see how the solution changes in time. We can animate the solution to our ODE using animateODESolution. To use it, first approximate a solution to the system (or scalar) equation using Euler. Then, tell animateODESolution what relationship you’d like to animate. For instance, to see how the \(x\)/\(y\) coordinates change, i.e., how the particle moves, for this IVP:
\[\begin{aligned} &x'=-y&\quad x(0)=1\\ &y'=x&\quad y(0)=0.5\\ \end{aligned}\]
we can do this:
If you are interested the \(x\) and \(y\) coordinates versus time, you could do
animateODESolution will show up on Project 3, but it is mostly a novelty. It’s hard to get anything quantitative out of these animations. They just look nice and make student’s say “ooh”.
3.2 Lesson 28: Introduction to Ordinary Differential Equations
3.2.1 Objectives
Understand the definition of an ordinary differential equation (ODE), including order.
Understand the difference between solving an equation and verifying a solution.
Understand that a solution of an ODE is a function, rather than a single input number.
Understand the difference between a specific solution and a general solution of an ODE.
Given an ODE and a function (or family of functions), determine whether the given function (or family of functions) is a solution to the ODE.
Given an ODE and a function containing parameter(s), find the parameter(s) that will make the given function a solution of the ODE.
3.2.2 In Class
Intro. This lesson serves as an introduction to ordinary differential equations (ODEs). The concept of a solution being a function may be difficult for some students to grasp.
General Solution. ODEs can have an infinite number of solutions and we describe them all with a parameter, \(C\).
Specific Solution is one particular function that solves the ODE.
Solution Verification. Given an ODE and a potential solution, we can check if the solution solves the ODE by plugging it into the LHS and RHS of the ODE. If LHS=RHS, then the given function is a solution to the ODE. This can be done algebraically, or by using
D, and/or usingplotFun.Order of an ODE. Order of the equation is the order of the highest derivative appearing on the unknown function. On this course we’ll focus on first order ODEs, but have them identify a second and/or third order ODE in this lesson.
3.2.4 Problems & Activities
Emphasize that differential equations have functions for solutions, not numbers. As an example, consider the differential equation \(\frac{dy}{dt}(t) = 2t\) from the reading. This equation has the general solution of \(y(t) = t^{2} + C\), because the derivative of \(y(t)\) is \(2t\), regardless of what value we choose for \(C\). If you want, spend some time on another example, like \(\frac{dy}{dt}(t) = \cos(t)\).
Consider the differential equation given by \(y^{'}(t) = 2ty(t)\). Emphasize that the right hand side of this expression is a function of the input AND output variables. In today’s lesson, we will not use methods for obtaining solutions. Rather, today, we will focus on determining whether a proposed function could be a solution to a given differential equation. In groups/at board/individually, have cadets determine whether the following two candidate functions are solutions to this ODE: \(y_{1}(t) = t^{2}\) and \(y_{2}(t) = e^{t^{2}}\).
Talk about the difference between a general solution and a specific solution. A general solution is a description of all the functions that solve an ODE. If we can integrate the RHS of the ODE, then we can find the general solution by including the constant of integration. We’ll talk about methods for solving ODEs where the RHS includes \(y\) in a later lesson. In the examples above, we found a few specific solutions just by guessing and checking.
There are lots of methods for verifying whether or not a proposed solution is indeed a solution to the ODE. If you can differentiate by hand and show algebraically that LHS=RHS, then great. If you need to use
Dto differentiate, that’s disappointing, but fine. If you need to plot LHS and RHS and see if the graphs are the same or different, fine. Take the ODE \[y'=y+e^t.\] is \(y_1=te^t\) a solution to this ODE? If we calculate LHS, \(y_1'\) we get \[LHS=e^t+te^t.\] We could just as well useDto calculate that fact.
## function (t)
## t * exp(t) + exp(t)If we calculate the RHS, \(y_1+e^t\), we get \[RHS = te^t+e^t.\] It looks to us like LHS=RHS, so \(y_1\) must be a solution to the ODE. If it were harder to see that LHS=RHS, we could plot the two and see if their graphs overlap.
plotFun(t*exp(t)+exp(t)~t,xlim=c(-5,5),lwd=3,col='maroon')
plotFun(y1(t)+exp(t)~t,lwd=3,col='cadetblue2',add=TRUE)
Those graphs seem to be on top of each other, so we reasonably assume that LHS=RHS, and so we think \(y_1\) is a solution to the ODE.
Talk about order of an ODE. The order of an ODE is the order of the highest derivative occurring on the unknown function. The order of \[y'=y+cos(t)\] is one. The order of \[y'=(y'')^4+y\] is two. The order of \[y'''+ty''+\sin(t)y'+4=0\] is three.
In this class we’ll only solve first order ODE, but order is a very common classification of ODEs, so it is good to have at least heard it before.
Have the students work on problems verifying whether or not proposed solutions are actually solutions to the ODE. The end of the chapter has a ton of them.
Also have students work on determining for what parameter value a proposed solution satisfies an ODE. For example, for what value of \(r\) does \(y_1=e^{rt}\) solve the differential equation \[y'=7y?\]
For what value of \(c\) does \(y_2=\cos(ct)\) solve the differential equation \[y''=-9y?\] For what value of \(d\) does \(y_3=6t^2+d\) solve the differential equation \[y'=2\left(\frac{y-4}{t}\right)?\]
3.3 Lesson 29: Analytical Solutions I
3.3.1 Objectives
Understand the definition of an initial value problem (IVP).
Given a pure‐time ODE, find the general solution through integration by hand.
Given an IVP, find the specific solution through integration by hand.
Describe the difference between the solution of an ODE and the solution of an IVP.
3.3.2 In Class
Introduction. Over the next few lessons, we will introduce methods for solving differential equations (analytical, graphical, and numerical). To facilitate their understanding, we will start with ODEs where the right hand side is a function of the input variable (pure-time ODEs). This should be very familiar given the material thus far. After introducing these methods in this context, we will transition to more complicated ODEs later.
Direct Integration. Review basic antiderivatives from the Antidifferentiation block. Remind them that when finding the antiderivative, the result has a \(+C\) component. Remind them why this is. In this lesson, we will practice finding a particular \(C\) given an initial value/condition. This will be known as an Initial Value Problem (IVP). You may see it referred to as either a “specific” or “particular” solution, we are trying to stick with specific in this course, though the two phrases are largely interchangeable.
No Analytical Solutions. We’ll end by introducing an IVP that is not possible to solve analytically. This will segue nicely to the next lessons: graphical and numerical approximation.
Problems & Activities
Start with review problems on ODEs and proposed solutions (the guess-and-check method). This should be done at boards or desks for the first few minutes.
Now transition to finding solutions analytically, rather than just guessing and checking. Have students find the general solution to an ODE that is integrable. (One possible example is \(y^{'} = e^{3t}\)). Find the general solution, and plot several different solutions to show the family of solutions.
If \[y'=e^{3t}\] then \[\int y'(t)\,dt = \int e^{3t}\,dt\] so \[y(t)+C_1=\frac{1}{3}e^{3t}+C_2.\] Combine the constants to get \[y=\frac13 e^{3t}+C\] where \(C=C_2-C_1.\)
Note: In the book, for a couple of examples, we carefully track the constant of integration that appears on the left-hand side of the ODE, as well as the one on the right, and then carefully combine them. However, after about two examples, we drop that practice for the more common practice of only including the constant of integration on one side of the ODE, understanding that we have already combined. Students do not need to write two constants of integration and then combine.
## function (t, C = 0) ## { ## exp(3 * t)/3 + C ## } ## <environment: 0x00000237e957d208>y=makeFun(1/3*exp(3*t)~t) plotFun(y(t)+0.5~t,tlim=c(-3,1),ylim=c(0,5),lwd=3,col="firebrick") plotFun(y(t)+1~t,add=TRUE,lwd=3,col="darkorchid") plotFun(y(t)+2~t,add=TRUE,lwd=3,col="darkgreen") mathaxis.on() grid.on()
Now attack some initial value problems. Unlike ODEs, which have general solutions, IVPs have specific solutions. Take the ODE above, and solve it a few times with different initial conditions. \[\frac{dy}{dt}=e^{3t}\] \[\text{a) } y(0)=2\] \[\text{b) } y(1)=3\] \[\text{c) } y(2)=4\]
After wrapping up that example, cadets should spend time at board/desks on the exercises at the end of the reading.
End class with an example that is not possible to integrate “by hand”, such as the one in the reading: \(y^{'} = e^{- t\sqrt{t}}\). Note that we will have to explore other methods to analyze this equation.
3.4 Lesson 30: Graphical Solutions
3.4.1 Objectives
Given a first‐order ODE, plot a subset of vectors from its slope field by hand.
Describe the difference between a slope field and a direction field.
Given a first‐order ODE, use R to plot a corresponding direction field.
Given a direction field for an ODE and an initial condition, sketch the corresponding specific solution.
Visually match a first‐order ODE to its corresponding direction field, including which are pure‐time and which are autonomous.
3.4.3 In Class
Visualizing ODEs. In this lesson, we’ll introduce methods to visualize first-order ODEs. Visualization is often helpful when dealing with ODEs that don’t have analytical solutions.
Vector Fields. We’ll discuss two types of vector fields – slope fields and direction fields. The only difference between the two is that slope fields show vectors of differing lengths, while direction fields normalize the magnitude of every vector; the only thing you pay attention to in a direction field is the direction the vector heads. Sometimes the vector magnitudes are too large or small which make the graph of a slope field difficult to read; however, the corresponding direction field graph should be legible (this is why we will be mainly focused on creating and using direction fields). We are going to define the slope field vectors to be \((1,f(t,y))^T\), where \(f(t,y)\) is the right-hand side of the ODE. Obviously there are other ways to define a vector with slope \(f(t,y)\), but that definition is common for slope fields.
Direction Fields in R
USAFACalcoffersplotODEDirectionField. Use that in class.
Problems & Activities
Start with a first-order ODE \(y^{'} = t^{2}\). Draw or project a grid on the board (or paper) and ask the students to draw the vectors that occur at various input values. If you want to produce slides with pre-drawn vectors, use the
place.vectorcommand, documented in Block 5 NTIs.Since we are not going to get into how to normalize a vector in this course, this exercise is actually drawing vectors for a slope field. I like to think of it like this: there are infintely many solutions to my ODE. I have no idea of a formula for any of them. However, I know that if the graph of one of the solutions happens to pass through \((t=1,y=0)\), then the ODE tells me that the slope of that solution must be \(y'(1)=1^2=1\). So I’ll draw a vector with slope 1 at the point \((t=1,y=0)\), showing that if a solution passes through this point, it must be headed in this direction.
Some \((t,y)\) input values you could use include (1,0), (1,1), (1,-1), (0,1), (-1,0).
You (or they) may wonder why the axis for the vectors fields are \(t\) and \(y\). We are not actually graphing \(y(t)\), instead we are graphing a snapshot of what what a solution to the ODE must do (i.e. the slope of its graph), if it passes through the vectors base points.
Use
plotODEDirectionFieldto create the direction field for the ODE.
when using plotODEDirectionField, you only specify the right-hand side of the ODE.plotODEDirectionFieldis very general, but as a result you need to make sure your right-hand side is defined as a function of both \(t\) and \(y\), in that order.The syntax of
plotODEDirectionFieldshould be very familiar by now. Adjust the \(t\)-axis withtlim=. Adjust the \(y\)-axis withylim=. Adjust the color withcol=.Now, pick an initial condition, e.g., \(y(0) = 2\), and use the direction field to sketch the solution by hand. Since you know that the vectors tell you the slope of the tangent line of the graph of the solution, you just “follow the arrows”. This is a lot like dropping a volleyball into a river (vector field) and watching the path (solution) taken by the volleyball. The spot where the volleyball first lands in the river is the “initial condition.”
Show them how to get
Rto plot the solutions for you using theic=option. Watch out, it is assumed that the initial condition is specified at the left-hand side of the graph, so if you want a solution corresponding to initial condition \(y(0)=2\), then you’ll need to make yourtlimstart at time 0.
Repeat the exercise: sketch a few vectors by hand, then make a direction field using R, then sketch a specific solution corresponding to an initial condition, then use R to sketch the solution corresponding to the initial condition, for another ODE, this one should have \(y\) on the right-hand side of the ODE. For example, something like: \[y'=ty(1-y)\] \[y(0)=0.1\] \[y(0)=2\] The only difference now is that in order to evaluate \(y'\), you’ll need to involve the \(y\)-coordinate of the base points of your vectors. Otherwise, it is exactly the same.
Talk about the difference between pure-time and autonomous differential equations. In a pure-time ODE, the right hand side only depends on \(t\). We saw how to solve some of those yesterday. In an autonomous ODE, the right-hand side only explicitly depends on \(y\). The direction fields for those will look different. For a pure-time ODE, the direction field will have vectors that do not change in the \(y\)-direction. The autonomous direction field will have vectors that do not change in the \(t\)-direction.
#pure-time
plotODEDirectionField(sin(t)~t&y,main="Pure-Time")
#autonomous
plotODEDirectionField(y*(3-y)~t&y,main="Autonomous")
#Niether
plotODEDirectionField(y*(3-y)+t~t&y,main="Neither pure-time nor autonomous")
- Practice identifying autonomous, pure-time, and neither ODEs from their corresponding direction field. Practice matching a direction field to an ODE. Exercises are in the chapter.
3.5 Lesson 31: Numerical Solutions
3.5.1 Objectives
Understand visually how Euler’s method uses the slope at one point to approximate the next point of the solution a step size away for an IVP.
Given a first‐order IVP and step size, approximate the solution at a specified input (or for a specified number of iterations) using Euler’s method by hand.
Understand that in general the accuracy of Euler’s method is improved by decreasing the step size.
Given a first‐order IVP and step size, use R to implement and/or visualize Euler’s method.
Given a pure‐time first‐order ODE, initial condition, and step size, compare the numerical approximation from Euler’s method to the analytical solution at a specified input in order to compute the error of the numerical approximation.
3.5.3 In Class
Analytical solutions not required. Even though we might not be able to compute the analytical solution to an ODE when the functions involved are not integrable, we can approximate solution values using the information we do have (much like we did with the vector fields).
Pick a starting point and a step size. We do need to specify a starting point (initial condition) and how much we want to increment the independent variable for each iteration of the numerical approximation method.
Euler’s Method. The method we are using in this course, Euler’s method, requires you to know three pieces of information: the current approximated solution value, the step size we want to take to get to the next approximated solution value, and the direction we need to step in given our current position. We can calculate the new “current position” and continue computing the next step using an iterative process.
3.5.5 Problems & Activities
There are no less than one million ways to motivate Euler’s method. I’ll describe two of them here. Do what works for you.
First, the thing that is most successful for me (Joe) personally is to say “Suppose I’m driving on I-25. Right now my car is at mile-marker 152, and I’m heading north (increasing mile markers) at 60 mph. Now, I might change speed as I drive due to traffic conditions, but if you had to guess, given the information you have, what mile marker do you think I’ll be at in 20 minutes (1/3 hour) from now?” They will mostly come up with the correct answer, \[152+1/3\cdot 60 = 172.\] If they can do that, then they have the gist of Euler’s method down. They took our current position, 152, and then computed how much they expected our position to change by taking the time-of-travel, 1/3, and multiplying by our speed, 60.
Let’s suppose that \(y(t)\) is my position (mile marker) on I-25, and that we know that \[y'=60+2\sin(y),\quad y(0)=152\] In this context the differential equation is like a traffic map that Google Maps draws. If you know where you are ,\(y\), the map (ODE) will tell you how fast you are going, \(y'\).
So, at time \(t=0\) we are initially at position \(152\). Assume \(t\) is in hours. Where do we expect to be after 1/6 hour (10 minutes)?
We have our current position, we need our current speed in order to mimic the calculation above. The ODE (traffic map) says that since we are at mile marker 152 our current speed is
\[y'=60+2\sin(152) \approx 61.87.\] Great, if our current speed is 61.87 then in 1/6 hours we expect to be at position
\[152+1/6\cdot 61.87=162.31.\]
So at time \(1/6\), we think we will be at position \(162.31\). So we think that \[y(1/6)\approx 162.31.\]
We can keep going, where do we think we’ll be at time \(1/6+1/6=1/3\). Well, at time \(1/6\) we are at approximately at position 162.31. According to the ODE then we would be travelling at speed
\[y'=60+\sin(162.31) = 59.13\] Therefore our position at time \(2/6\) should be about
\[162.31+1/6\cdot 59.13 \approx 172.16.\] What we are saying is that we think \[y(2/6)\approx 172.16.\]
We could clearly keep doing this over and over again. Formally, we are letting \(\Delta t\) be the step size we are taking, also called step length or time step. In the example above, \(\Delta t=1/6\). We set \(t_n\) be the times at which we are approximating our positions. In the example above, \(t_0=0\), \(t_1=1/6\), \(t_2=2/6\), and so on. We define \(y_n\) to be our approximation to \(y(t_n)\). So we had \[\begin{aligned} &y(0)\approx y_0=152\\ &y(1/6)\approx y_1 = 152+1/6\cdot 61.87=162.31\\ &y(2/6)\approx y_2=162.31+ 1/6\cdot 59.13 = 172.16.\end{aligned}\]
The pattern we are following is that \[y_{n+1} = y_n + \Delta t y'_{n},\] and in order to compute \(y'_n\) we need to use the ODE, using the values of \(y_n\) and \(t_n\).
Students will understand that taking smaller time-steps is better. Long time steps assume that I’m driving at a constant speed for a long time, which might not be true. Smaller time steps give me less time to change my speed.
Another option is to go with the direction-field approach. Start with an ODE that we can solve analytically, such as \[y^{'} = \sin(t),\quad y(0)=1.\] Use
plotODEDirectionFieldto plot the direction field in the neighborhood of \(t=0,y=1\).
We know from the initial condtion that the specific solution we are after passes through the point \((t=0,y=1)\). Place a dot at that point to illustrate.plotODEDirectionField(sin(t)~t&y,tlim=c(-1,1),ylim=c(0,2)) mathaxis.on() grid.on() plotPoints(1~0,pch='.',cex=10,col="red",add=TRUE)
We know that if the solution passes through \((t=0,y=1)\) then it must pass parallel to the vectors in the direction field, whose slopes are computed with the ODE. The slope of the tangent line to our solution must be\[y'=\sin(0)=0.\] Draw a line with slope 0 through the point \((t=0,y=1)\).
plotODEDirectionField(sin(t)~t&y,tlim=c(-1,1),ylim=c(0,2)) mathaxis.on() plotPoints(1~0,pch='.',cex=10,col="red",add=TRUE) plotFun(0*x+1~x,lwd=3,col="dodgerblue4",lty=2,add=TRUE)
We expect our solution to stay close to that line, since it is the tangent line to our solution, at least for a little while. Where is this tangent line at time \(0.25\)? Using the equation for the tangent line (it has slope 0 and passes through (0,1)) we find that we expect it to be at\[y=0\cdot(0.25-0)+1=1.\] So we think that our solution will pass close to the point \((t=0.25,y=1)\).
plotODEDirectionField(sin(t)~t&y,tlim=c(-1,1),ylim=c(0,2)) mathaxis.on() plotPoints(1~0,pch='.',cex=10,col="red",add=TRUE) plotFun(0*x+1~x,lwd=3,col="dodgerblue4",lty=2,add=TRUE) plotPoints(1~0.25,pch='.',cex=10,col="red",add=TRUE)
In order to keep following the slope field, we’ll redraw a new tangent vector to our solution, assuming that our solution really does pass through \((t=0.25,y=1)\). The slope of the tangent line would be\[y'=\sin(0.25)=0.247\] Draw the tangent line to our solution; it passes through \((t=0.25,y=1)\) and has slope 0.247.
plotODEDirectionField(sin(t)~t&y,tlim=c(-1,1),ylim=c(0,2)) mathaxis.on() plotPoints(1~0,pch='.',cex=10,col="red",add=TRUE) plotFun(0*x+1~x,lwd=3,col="dodgerblue4",lty=2,add=TRUE) plotPoints(1~0.25,pch='.',cex=10,col="red",add=TRUE) plotFun(0.247*(x-0.25)+1~x,lwd=3,col="dodgerblue4",lty=2,add=TRUE)
Once again, now we have a tangent line to our solution. Where do we expect our solution to get to at time \(t=0.5\)? If we just use the equation of the tangent line we’ll see that we expect our solution to be at
\[y=0.247\cdot(05-0.25)+1= 1.062.\]So we think our solution will pass close to \((t=0.5,y=1.062)\).
plotODEDirectionField(sin(t)~t&y,tlim=c(-1,1),ylim=c(0,2)) mathaxis.on() plotPoints(1~0,pch='.',cex=10,col="red",add=TRUE) plotFun(0*x+1~x,lwd=3,col="dodgerblue4",lty=2,add=TRUE) plotPoints(1~0.25,pch='.',cex=10,col="red",add=TRUE) plotFun(0.247*(x-0.25)+1~x,lwd=3,col="dodgerblue4",lty=2,add=TRUE) plotPoints(1.062~0.50,pch='.',cex=10,col="red",add=TRUE)
One more time for good measure. If our solution passes through \((t=0.5,y=1.062)\) then it must have slope\[y'=\sin(0.5)=0.479.\]
plotODEDirectionField(sin(t)~t&y,tlim=c(-1,1),ylim=c(0,2)) mathaxis.on() plotPoints(1~0,pch='.',cex=10,col="red",add=TRUE) plotFun(0*x+1~x,lwd=3,col="dodgerblue4",lty=2,add=TRUE) plotPoints(1~0.25,pch='.',cex=10,col="red",add=TRUE) plotFun(0.247*(x-0.25)+1~x,lwd=3,col="dodgerblue4",lty=2,add=TRUE) plotPoints(1.062~0.50,pch='.',cex=10,col="red",add=TRUE) plotFun(0.479*(x-0.5)+1.062~x,lwd=3,col="dodgerblue4",lty=2,add=TRUE)
We expect our solution to stay close to its tangent line for a little bit. At time \(t=0.75\) the tangent line has value
\[y=0.479\cdot(0.75-0.5)+1.062=1.182.\]
So, we think that our solution must pass approximately through \((0.75,1.182)\).plotODEDirectionField(sin(t)~t&y,tlim=c(-1,1),ylim=c(0,2)) mathaxis.on() plotPoints(1~0,pch='.',cex=10,col="red",add=TRUE) plotFun(0*x+1~x,lwd=3,col="dodgerblue4",lty=2,add=TRUE) plotPoints(1~0.25,pch='.',cex=10,col="red",add=TRUE) plotFun(0.247*(x-0.25)+1~x,lwd=3,col="dodgerblue4",lty=2,add=TRUE) plotPoints(1.062~0.50,pch='.',cex=10,col="red",add=TRUE) plotFun(0.479*(x-0.5)+1.062~x,lwd=3,col="dodgerblue4",lty=2,add=TRUE) plotPoints(1.182~0.75,pch='.',cex=10,col="red",add=TRUE)
Now we can hopefully see a pattern emerging. Again, let \(y_n\) be our approximation to \(y(t_n)\). We started with the initial condition, so \(y_0=y(t_0)=y(0)=1\). We chose our step size to be \(\Delta t=0.25\), and estimated the value of \(y\) at times \(t_1=t_0+0.25=0.25\), \(t_2=t_1+0.25=0.5\), \(t_3=t_2+0.75\). We calculated
\[\begin{aligned}
&y(t_1)\approx y_1 = y'_0\cdot \Delta t + y_0 = 1\\
&y(t_2)\approx y_2 = y'_1\cdot \Delta t + y_1 = 1.062\\
&y(t_3)\approx y_3 = y'_2\cdot \Delta t + y_2 = 1.182.\end{aligned}\]So the pattern seems to be \[y_{n+1}=y'_n\Delta t + y_n,\] where \(y'_n\) is calculated from the ODE using the values of \(y_n\) and \(t_n\).
Either of the above introductions is fine, or you can do something else. At this point, have them compute a few iterations of Euler’s method by hand to compute an approximate solution to \[y'=t-y,\quad y(0)=2\] using step sizes of \(\Delta t=0.25\). Have them follow a table like the one in the book to keep things organized. Have them determine how many steps it will take them to approximate the value of the solution at time \(t=4\).
Obviously, we do not want to do all the iterations by hand. Besides, computers that approximate solutions for IVPs can use a much smaller step size in order to increase the approximation’s accuracy. Show them how to use
Eulerto approximate the solution to this problem. You can either only specify the right-hand side:## t y ## 1 0.00 2.0000000 ## 2 0.25 1.5000000 ## 3 0.50 1.1875000 ## 4 0.75 1.0156250 ## 5 1.00 0.9492188 ## 6 1.25 0.9619141 ## 7 1.50 1.0339355 ## 8 1.75 1.1504517 ## 9 2.00 1.3003387 ## 10 2.25 1.4752541 ## 11 2.50 1.6689405 ## 12 2.75 1.8767054 ## 13 3.00 2.0950291 ## 14 3.25 2.3212718 ## 15 3.50 2.5534538 ## 16 3.75 2.7900904 ## 17 4.00 3.0300678or you can specify the full equation, using
_tto denote the derivative.## t y ## 1 0.00 2.0000000 ## 2 0.25 1.5000000 ## 3 0.50 1.1875000 ## 4 0.75 1.0156250 ## 5 1.00 0.9492188 ## 6 1.25 0.9619141 ## 7 1.50 1.0339355 ## 8 1.75 1.1504517 ## 9 2.00 1.3003387 ## 10 2.25 1.4752541 ## 11 2.50 1.6689405 ## 12 2.75 1.8767054 ## 13 3.00 2.0950291 ## 14 3.25 2.3212718 ## 15 3.50 2.5534538 ## 16 3.75 2.7900904 ## 17 4.00 3.0300678Euler’s method produces a series of points that approximate the values of the solution to the IVP. Show students how the approximation improves as \(\Delta t\) gets smaller. The analytical solution to \[y'=\cos(t),\quad y(0)=0\] is \(y=\sin(t)\). Plot the analytical solution as well as the output of Euler’s method for various step sizes.
plotFun(sin(t)~t,tlim=c(0,10),lwd=4,col='darkmagenta') mathaxis.on() solution=Euler(cos(t)~t&y,tlim=c(0,10),ic=0,stepSize = 0.5) plotPoints(y~t,data=solution,col='darkseagreen4',add=TRUE) solution=Euler(cos(t)~t&y,tlim=c(0,10),ic=0,stepSize = 0.1) plotPoints(y~t,data=solution,col='darkred',add=TRUE) solution=Euler(cos(t)~t&y,tlim=c(0,10),ic=0,stepSize = 0.05) plotPoints(y~t,data=solution,col='deepskyblue4',add=TRUE)
Use the example above to calculate the error in Euler’s method. The error at time step \(n\) is defined to be \(|y(t_n)-y_n|\). Consider the error in Euler’s method at time 9. Using step size 0.25, time 9 is \(t_{36}\).
The value of the analytical solution at time 9 is 0.4121185. The Euler approximation to \(y(9)\) is \(y_{36}=0.6488611\). The error in Euler’s method at time 9 using step size 0.25 is \[|0.4121185-0.6488611|=0.2367426.\]
## t y ## 32 7.75 1.1014389 ## 33 8.00 1.1273875 ## 34 8.25 1.0910125 ## 35 8.50 0.9945755 ## 36 8.75 0.8440725 ## 37 9.00 0.6488611## [1] 0.2367426If we reduce the step size to 0.1, then the Euler approximation we are after is \(y_{90}\). The value of \(y_{90}\) is \(y_{90}=0.5073315\). The error in Euler’s method at time 9 using step size 0.1 is \[|\sin(9)-y_{90}|=|0.4121185 - 0.5073315|=0.095213\]
## t y ## 896 8.95 0.4665841 ## 897 8.96 0.4576902 ## 898 8.97 0.4487510 ## 899 8.98 0.4397674 ## 900 8.99 0.4307403 ## 901 9.00 0.4216707## [1] 0.09521301Remember, if we actually had access to the analytical solution of an IVP, then we wouldn’t need Euler’s method. This calculation is for demonstration purposes only. In real life, you never have an exact formula for how much error there is in your Euler approximation.
Repeat this process with an IVP that cannot be solved analytically. Perhaps \[y^{'} = y\sin(t + \sqrt{t}) \text{ with } y(0) = 1\].
3.6 Lesson 32: Analytical Solutions II
3.6.1 Objectives
Recognize and distinguish between autonomous ODEs of the form \(y' = f(y)\) and pure‐time ODES of the form \(y' = f(t)\).
Understand the definition of a separable ODE. Given a first‐order ODE, determine if it is separable.
Given a separable ODE, find the general solution using separation of variables by hand.
Given an IVP, use separation of variables to find the general and specific solutions by hand.
3.6.3 In Class
Separable differential equations. Up until now, we’ve only worked with ODEs whose rate of change was only dependent upon the variable \(t\), which we can solve by integration techniques from Block 4 or using
antiD. The next two days will be spent working with separable ODEs. We define a first-order ODE as separable if it can be written in the form \[y'=f(t)g(y).\]Analytical solutions. In order to solve an ODE of the form \(y^{'} = f(t)g(y)\) analytically, we must use algebra to rewrite the equation. We need to get all the terms involving \(y\) multiplied by \(\frac{dy}{dt}\), and all of the \(t\) terms on the other side of the equation, not multiplied by \(\frac{dy}{dt}\). A common error cadets will make is that they will use algebra to create a sum of functions of \(y\) on the LHS which does not work. For example, consider the ODE \(y^{'} = y + 1\). In order to set up the problem correctly for integration using separation of variables so that you have a product of functions of \(y\) and \(\frac{dy}{dt}\) on the LHS, you must divide both sides by \((y + 1)\). Simply subtracting \(y\) from both sides does not result in a product of functions of \(y\) on the LHS. It is also important to remember that y is an unknown function. The operation \(\int y\,dt\) makes perfectly good sense. It means “whatever the antiderivative of the function \(y\) is”. However, you aren’t going to a more simple answer than what you have – you don’t know what \(y\) is. The answer is certainly not \(\frac{1}{2}y^2\), as you can check in the case that \(y=\cos(t)\).
Separation of Variables. Assuming we can algebraically separate the ODE into the product of functions of \(y\) on the LHS and only constants on the RHS, we can attempt to find a solution to the ODE by integrating both sides of the equation. We use a trick involving the chain rule to get a usable expression on the LHS of the ODE.
3.6.5 Problems & Activities
Start with a reminder of where we’ve been in this block. We’ve studied ODEs of the form \(y^{'} = f(t)\). Explain that today we’ll widen our toolbox by exploring differential equations in the form \(y'=f(t)g(y)\).
Practice identifying whether or not an equation is separable. We won’t intentionally try to trick them by asking them if something difficult is separable, but they should be able to recognize that \[\frac{dy}{dt}=yt-t\] is separable by factoring out at \(t\).
Following the motivating examples in the book, solve a couple of problems using integration by parts. The big thing is to separate the variables, and then let the chain rule save the day when you are faced with integrating something like \(y^3\frac{dy}{dt}\). Even though you don’t know what function \(y\) is, the chain rule tells you that the anti-derivative of \(y^3\frac{dy}{dt}\) is \(\frac14 y^4\).
After a couple of examples, it is completely fine to only write the constant of integration on one side or the other of the ODE.
After a couple of examples, if you want to think of “multiplying both sides by \(dt\)” rather than “canceling \(dt\)” as in the book, that is completely fine.
After a couple of examples, it is completely fine to “cancel” the \(dt\)’s and write \[\int y \frac{dy}{dt}\,dt=\int y\,dy.\]
In this section it is completely fine to write \[\int \frac{1}{y}\,dy = \ln(y).\] When this comes up, we’ll always assume that \(y>0\). Once you do that, then you can omit any sign analysis on terms like \(e^{C_1}\). Just turn \(e^{C_1}=C\), as in the book. Yes, I know perfectly well that it should be \(\ln|y|\) and that then we should deal with cases to turn \(e^{C_1}\) into a constant \(C\) with unrestricted sign. For this class that is more trouble than it is worth. They’ll deal carefully with this in an ODE course.
It is totally fine for students to use
antiDto do integration. Disappointing, but fine.Though we are doing separation of variables, we should focus on autonomous differential equations. After a couple non-autonomous examples, switch over to autonomous.
Have students work on problems from exercises in the section. There are plenty of them. Find both general solutions and specific solutions to IVPs.
3.7 Lesson 33: Analytical Solutions II
3.7.1 Objectives
Given an IVP, use separation of variables to find the general and specific solutions by hand.
Given an IVP in the context of a real‐world scenario, find the solution and interpret its results within the context of the scenario.
3.7.3 In Class
Review. This is a continuation of last lesson. Review separability, and how we can not find a more simple expression for \(\int y(t)\,dt\), at least without knowing what \(y(t)\) is, but we *can find a simpler expression for \(\int y(t)\frac{dy}{dt}\,dt\).
Practice. Much of this section should be devoted to practicing finding solutions to ODEs of the form \(y^{'} = f(y)\). There is time for you to introduce your own examples, or you can work with the ones from the book.
Applications Spend a significant amount of time on the application problems in the back of the book. Students should be able to
- Identify the correct initial conditions from a word problem.
- Use the solution to the IVP to answer a context-specific question.
- Students are welcome to solve these questions using either algebra, or by using tools like
findZeros. Both are demonstrated in the book.
3.7.5 Problems & Activities
Spend about half the class having the students solve exercises from the back of the chapter. They can use
antiDif they need to.Cover the two “application” problems from the book, 7.4.3 and 7.4.4, paying special attention to the “Additional Insights”. Students will need to be able to identify the correct initial conditions from a word problem, solve the IVP, and then use the solution to the IVP to solve an additional question. Try to cover the full examples, as they will be asked similar questions in the homework. More examples can be found in HW 35, or swing by my office and borrow an ODE book.
3.8 Lesson 34: Modeling with ODEs I
3.8.1 Objectives
Given a description of a rate and an initial condition, identify an IVP that matches that description.
Set up and/or solve an IVP that models the changing temperature of an object placed in a fixed environment using Newton’s law of heating/cooling.
Set up and/or solve an IVP that models a growing population without a carrying capacity.
Set up but do not solve an IVP that models a growing population with a carrying capacity.
Set up an IVP that models the changing amount of an investment or debt account, given an interest rate and repeated investments or repeated payments
3.8.3 In Class
Review of ODEs. Start with a reminder of where we’ve been the last few lessons. We’ve worked with separable ODEs including pure-time and autonomous ODEs. We’ve found solutions to IVPs involving these using analytical methods, and approximated solutions using graphical and numerical methods. Today we’ll start talking about why we care about ODEs. It turns out that ODEs provide us with a language that clearly describes many natural phenomena. Many formulas which students are familiar with (projectile motion with no wind resistance, population growth being two) aren’t obvious at all. Just by looking at the formula, it isn’t clear that the formula describes anything at all. Those formulas are just solutions of ODEs, and if you know the language of ODEs it is clear that the ODE is describing something.
Putting it Together. Now we will focus on using ODEs for the purpose of modeling. Some of this should be review from block 1, but we will discuss basic growth/decay problems, and then move to growth problems with carrying capacity using the logistic equation. Demonstrate using the squirrel example from the textbook.
The book contains lots of material. This is partly due to Rob’s original material, and partly due to my own tendencies. Obviously you won’t be able to cover all of those models in one day. Focus on exponential growth (population and interest), logistic growth, and Newton’s Law of Cooling.
Today is not about solving. It is about a) interpreting what the ODE says, so that students can set up similar IVPs when given information, and b) using the solution to an IVP to answer a question. You shouldn’t spend a lot of time solving these equations. They should, at this point, remember the solution to \(y'=ky\). The solution to \(T'=k(T_a-T)\) you can do quickly, or assign as an exercise. Don’t spend time solving logistic models. We didn’t do partial fractions decomposition in this class. Plus, that’s not the point. Writing down and interpreting the models is and using the solutions of the models to answer questions is the point.
Examples Work through as many examples as you can. Blow over the solving of the ODEs. Talk about interpreting what the ODEs say and using the solutions of ODEs to answer questions.
3.8.5 Problems & Activities
Start with exponential growth. Interpret the differential equation \(P'=kP\) as “the rate of growth of the population is always proportional to the current population”. That means that \(k\) must be of dimension time\(^{-1}\). Write down the general solution, \(P=P_0e^{kt}\), and remind them, from Block 1, of the relationship between \(k\) and doubling time. With the general solution in hand, work through a few of the questions in Example 8.1.2.1.
Another nice example of exponential growth is continuously compounded interest. Interpret the differential equation \(P'=rP+I\) as “the rate at which the account value grows is the rate at which interest is accumulated, plus the rate at which money is deposited”. The rate at which interest is accumulated is modeled as \(rP\), where \(r\) is the interest rate. The solving of that ODE should be a forgone conclusion now. Solve if if you feel you have time, but feel free to write down the general solution, \(P=\left(P_0+\frac{I}{r}\right)e^{rt}-\frac{I}{r}\). Then, do Example 8.1.3.1. Note that using Euler’s method to solve the ODE numerically, and interpret the output of that is also totally fine.
Now discuss a similar growth process, but note that in many cases, a quantity will realistically not be allowed to grow infinitely. At some point, a “carrying capacity” is reached, and there is simply no more room to grow. Consider plants or animals in a finite space. These problems are often modeled using a logistic equation: \[P^{'} = kP\left(1 - \frac{P}{M}\right)\], where \(k\) is the growth rate when the population is small and growing exponentially, and \(M\) is the carrying capacity. Don’t try to solve this ODE. Instead follow Example 8.2.0.1, and use
Eulerto do the solving work for you.Now, discuss Newton’s Law of Cooling. Again, solving this analytically in class isn’t the point, it should be do-able on homework. Feel free to either assert that the solution is \(T=(T0 - T_A)e^{-kt} +T_A\), or just skip completely! It is interesting how much information we can extract without solving the ODE analytically. Cover Example 8.4.1.
As time remains, cover more examples. There are plenty in the book.
3.9 Lesson 35: Modeling with ODEs II
3.9.1 Note
Today is labeled as a modeling day, but in reality it is about equilibrium solutions, phase lines, and stability. That’s just how things shook out.
3.9.2 Objectives
Understand the definition of an equilibrium solution and its stability (stable or unstable).
Given a direction field for an ODE, visually estimate the equilibrium solutions of the ODE and determine the stability of each equilibrium.
Given an autonomous first‐order ODE, algebraically solve for the equilibrium solutions of the ODE by hand.
Given an autonomous first‐order ODE, create a phase line and use it to determine the stability of each equilibrium solution.
Given an IVP, describe the long‐term behavior of the solution (as time goes to infinity).
3.9.4 In Class
Autonomous ODEs. Start with a reminder that we have been focusing on ODEs of the form \(y' = f(y)\). These are known as autonomous ODEs. In order to find the value of \(y'\), you only need to know the current value of \(y\), not what time it currently is.
Equilibrium Solutions. Define an equilibrium solution as functions that a) solve the ODE, and b) are constant. In order to find equilibrium solutions, we need to solve \[0=y'=f(y)\] for \(y\). Since we are looking for constant solutions, these should turn out as constant numbers.
Direction Fields and Phase Lines. Remind them that we can visualize these ODEs using
plotODEDirectionField. Recall that autonomous ODE have direction fields that do not change along the \(t\) axis, so we can compress all of the information in the direction field into a phase line.Once we have phase lines, define an equilibrium solution as stable (note, not asymptotically stable) if initial conditions near the equilibrium solution correspond to solutions that stay near the equilibrium solution (don’t tend toward another equilibrium solution, nor \(\pm \infty\)). In the scalar case, that will mean that the solution starting near the equilibrium solution will tend toward the equilibrium solution. Otherwise, the equilibrium solution is unstable.
3.9.6 Problems & Activities
Start with a fairly simple autonomous ODE, such as \(y^{'} = y^{2} - 1\). Using
plotODEDirectionField, note that the behavior of this ODE does not vary based on \(t\). Using \(N=21\) as an option puts arrows right on the equilibrium solution. It will look ok if you omit that option.
Add a few initial conditions to the direction field. Note that they all have one of a few basic behaviors: converge to -1, go to \(\infty\), go to \(-\infty\), and the very special one that stays +1.
Since the behavior of the solutions to autonomous ODEs is so simple, we summarize the information in the direction field with a phase line. Draw a phase line for this ODE.
Note that the initial conditions +1 and -1 were very special and lead to constant solutions. Constant solutions have zero derivative. Constant solutions to the ODE are called equilibrium solutions. To find the equilibrium solutions we could
Guess, based off the phase line. Using basic algebra or viewing the direction field or phase line, have the cadets tell you what the equilibrium solutions are for this ODE. In this case, the equilibrium solutions are \(y = - 1\) and \(y = 1\).
Use algebra. We would need to solve \[0=y'=y^2-1\] for \(y\). That would give us \[y^2=1\] so \[y=\pm 1.\]
Use
findZeros. While this is a little disappointing, it is fine.
## y ## 1 -1 ## 2 1There are two types of equilibrium solutions we’ll discuss for autonomous ODEs: stable and unstable. An equilibrium solution is stable if solutions starting nearby stay nearby (doesn’t go to another equilibrium solution, nor \(\pm \infty\). Similarly, an equilibrium solution is unstable if nearby solutions veer away from that equilibrium. Using our example, have the cadets classify \(y = - 1\) and \(y = 1\).
If time permits, talk about stability in the context of modeling, say for population models. You could certainly return to the logistic equation example, Example 8.2.0.1 and draw the phase line for that model. Talk about the difference between the stable equilibrium solution at \(P=225\) and the unstable one at \(P=0\). Students love considering zombies, so don’t shy away from considering \(P<0\).
The rest of the class should be dedicated to board/group work. There are plenty of exercises in the book.
3.10 Lesson 36: Systems of ODEs I
3.10.1 Objectives
Understand that a system of ODEs is a collection of multiple ODEs with a single input and multiple outputs.
Understand that the solution of a system of ODEs is a vector of functions, rather than a single function.
Given a system of two first‐order ODEs and a proposed solution, determine if the proposed solution satisfies the system.
Given a system of two first‐order ODEs, find all equilibrium solutions of the system algebraically by hand or using R.
3.10.2 Note
The official syllabus says that we should find equilibrium solutions by hand. That is a typo – using R is fine also.
3.10.4 In Class
- Systems of ODEs. Introduce systems of ODEs. Just like with systems of equations, you want functions that satisfy both equations.
For example, if a system of equations is \[x+y=10\] \[x-y=0\] you want to find numbers, \(x\) and \(y\) with the property that both \(x+y=10\) and \(x-y=0\).
Similarly, if the system of ODEs is \[x'=y\] and \[y'=x+4\] you are looking for two functions \(x\) and \(y\), that satisfy both equations.
General Solutions. Once we move to systems of ODEs, our general solutions will have more constants in them. This is a fact, you can try to motivate that fact or not – I would not. In general, there are as many constants in the general solution as there are equations in the system. We will stick to two equations in the system.
Verifying Proposed Solutions. Spend some time verifying that proposed solutions (which are now pairs of functions) are indeed solutions to the system of ODEs.
Initial conditions. Now that we have multiple constants in our general solution, we’ll need to provide multiple initial conditions in order to pin down a specific solution. It is fine if they use
findZerosto solve for the constants needed to solve the initial conditions.Equilibrium Solutions. Very analogous so equilibrium solutions for scalar ODEs, equilibrium solutions are pairs of functions that solve the system of ODEs and are both constant. Spend some time finding equilibrium solutions. It is completely fine to rely on
findZeroshere.
3.10.6 Problems and Activities
This is a pretty full lesson, you’ll need to spend your time wisely.
Start with a system of ODEs, and verifying a proposed solution. Examples 1-4 provide this aplenty. Only spend enough time here to make sure that they get the idea that solutions to systems of ODEs are pairs of functions. Perhaps 2 examples.
Initial Conditions. Again, do just barely enough here. This doesn’t even align with an objective, but it does help remind the students that we are now dealing with pairs of functions. Example 5 is a good one, or add some initial conditions to any of Example 1-4. I would do 1-2 examples here.
Now move on to equilibrium solutions. Analogously to finding equilibrium solutions for scalar ODE, you are now looking for solutions to the system (so, pairs of functions) that are both constant. If both functions are constant, then both of their derivatives are zero.
For instance, consider the system of ODEs \[\begin{aligned} x'&=y-4\\ y'&=x^2-9\end{aligned}\]
If we have two constant functions (equilibrium solutions) that solve this system, then \[\begin{aligned} 0&=y-4\\ 0&=x^2-9\end{aligned}\]
This one is pretty easy to solve, there are two solutions: \(x=3,y=4\) and \(x=-3 ,y=4\). However, it is easy to make up much harder problems. It is fine to use
findZerosto solve these problems.## x y ## 1 -3 4 ## 2 3 4Now move on to more examples. Example #6 is a predator-prey example. Focus on finding the equilibrium solutions, and then interpreting the equilibrium solutions as populations. If the population ever manages to hit these states, then the population will never change. Example 7 is a competing species example. Interpret it similar to above.
Example 8 is a classic pendulum problem. I particularly like this one, though I expect it is more difficult for students to interpret. Don’t shy away! For this problem I would start by asking students, based on physical intuition, what they expect the equilibrium solutions to be. That is, what configuration could the pendulum get into such that the pendulum will never move from the configuration? Someone will surely say hanging straight down and at rest. That gets you \(a=0\), \(v=0\). Next, ask about \(a=2\pi\), \(v=0\), is that just as good? As far as the equations are concerned, yes, that’s just as good (but not a physically different configuration). Finally, after some prompting, hopefully someone will say that the pendulum might be balanced upside down. That’s correct, and corresponds to \(a=\pi\), \(v=0\). Finally, use
findZeros(or trigonometry) to solve for the equilibrium solutions.Example 9 or any of the exercises at the end of the chapter, as time allows.
3.11 Lesson 37: Systems of ODEs II
3.11.1 Objectives
Given a system of two first‐order ODEs, find all equilibrium solutions of the system algebraically by hand or using R.
Given a phase plane for a system of two first‐order ODEs, sketch a specific solution given an initial condition, visually identify any equilibrium solutions and determine the stability of each equilibrium (stable or unstable).
Visually match a system of two first‐order ODEs to its corresponding phase plane.
Use R to produce a phase plane for a given system of two first‐order ODEs.
3.11.2 Note
The official course objectives say that students should be able to find equilibrium solutions by hand. That is a typo. It is fine if they use
findZeros.The objectives say they should be able to visually match a system of ODEs to the phase plane. This is very similar to this example from direction fields, as well as exercises 5-12 from Exercises 5.3. If you have time, do a few of these, there are problems at the end of the chapter. However, I’m pretty comfortable asking the cadets to do it on their own.
3.11.4 In Class
Do one warm up problem of finding an equilibrium solution to a system of ODEs. Pick any example from the text that you haven’t used yet.
Now start sketching a phase plane. Follow the development in the book. Each point in the plane is used to represent a possible state (in the book, populations). From each potential state that we could reach, we ask the ODE how the population would change from here. Then we draw a vector heading in the direction the ODE prescribed. This would be a good exercise to do on the board.
Once you’ve done a few vectors in the phase plane, use
plotPhasePlaneto get the full picture.Talk about the difference between phase planes and direction fields. They look a lot the same, they are used to do the same thing. However, there is no time axis on a phase plane. You can indeed “follow the flow” on a phase plane to get an approximate solution to an IVP, but you won’t be able to read off the phase plane the time at which the solutions achieve each state.
Finally, stability. Follow the books lead and use examples with physical interpretations. Interpret the stability of an equilibrium solution as “If the initial condition changes from the equilibrium initial condition by a little, the solution to the IVP will/won’t change a lot”.
3.11.6 Problems and Activities
- After a warmup of finding an equilibrium solution, develop the notion of the phase plane. It is much like the phase line, but since there are two quantities changing we need two axes, thus the phase plane. Follow the development in the book.
- Once you’ve drawn a few vectors on the board, use
plotPhasePlaneto help you along.
Once you have a phase plane in mind, ask how the populations might evolve if the initially start at \(H=3\), \(S=100\)? Just like a direction field, the solution will “flow along” the arrows. This is a place where I would project the phase plane on the board, draw a big dot at position \(H=3\), \(S=100\), and then have them tell me how to flow along the directions.
Now, interpret this curve you have drawn. The curve passes through several points, for example \(H=4.26\), \(S=55.32\). That means that at some time we expect there to be 4.26 hawks and 55.32 squirrels living in the park. Interpret a few other points. Have the students describe in words how this population behaves. Something like “Initially, the number of hawks in the park increases, while the number of squirrels decreases. Once the population reaches roughly 4.25 hawks and 50 squirrels, though, the number of hawks starts decreasing (they are running out of squirrels to eat), and while the number of squirrels continues to decline. Once the population reaches about 2 hawks and 30 squirrels, the number of squirrels start to rebound (there aren’t as many hawks to eat them), while the hawk population continues declining”, and so on.
You can actually draw the curve you sketched on the board by giving
plotPhasePlanean initial condition to start from.
plotPhasePlane(c(-0.5*H+0.009*S*H,
0.3*S*(1-S/90) - 0.09*S*H)~H&S,
xlim=c(0,6),ylim=c(0,130),ics=c(3,100))
Finally, interpret that point in the center of the swirl. It is an equilibrium solution!
Now, ask what would happen if the initial condition were not exactly on the equilibrium solution. This one is a little hard to tell, but it looks like no matter where you start, the populations will eventually swirl in to the equilibrium state. That makes it a stable equilibrium solution.
Consider the other equilibrium solutions by plotting a the phase plane over a bigger window, as done in the book. Ask them to visually identify the equilibrium solutions, and ask them to visually classify the stability of the equilibrium solutions.
- As time allows, pursue some more examples. They might especially like Lanchester’s Laws, or the Pendulum Example
3.12 Lesson 38: In-Class Project
3.12.1 Objectives
Given a real-world scenario and a corresponding system of first-order ODEs, interpret the meaning of the parameter(s) and/or term(s) in the system.
Given a real-world scenario and a corresponding system of first-order ODEs, determine how the parameter(s) and/or term(s) of the system should change if the scenario changes.
Given a real-world scenario and a corresponding system of first-order ODEs, find all equilibrium solutions algebraically by hand.
Given a real-world scenario and a corresponding system of first-order ODEs, use R to produce a phase plane and visually determine the stability of each equilibrium.
Given a system of first-order ODEs, an initial condition, and a step size, use Euler’s method to approximate the solution of the system at a specified input (or for a specified number of iterations), by hand or in R.