Earlier we saw how to use systems of first order differential equations as a way to represent second order differential equations. This was an essential step in order to use numerical methods such as Euler's Method and RK4 to estimate specific solutions. [br][br]Now we'll study systems of first order differential equations in their own right, independent of any relationship they have to second order differential equations. [br][br]The theory of systems of first order differential equations is beautiful and parallels the theory of eigenvalues and eigenvectors from linear algebra. We'll take a look at this in later lessons. For today we will start by looking at an application to get ourselves acquainted with the topic: [br][br][list][*]A linear constant coefficient model of two species -- predator and prey. [/*][/list][br]There are more sophisticated non-linear, non constant coefficient models of this process, but we won't worry about it while we're just getting acquainted with systems of differential equations.

The [url=https://www.stowelandtrust.org/conserved/properties/bouchardfarmlandmarkmeadow]Landmark Meadow[/url] in Stowe, VT is a great habitat for owls and mice. The two species have a predator and prey relationship; the mice are the prey of the owl, who are predators. Naturally, an abundance of mice is good for the owls, while a dearth of owls is good for the mice. Similarly, too many owls is bad for the mice. In this way, the populations of mice and owls are inter-related. This makes these populations a perfect example of something that can be modeled with a system of equations.[br][br]To do so, we need make a few assumptions about the way the presence of each species affects the rate of change of each species. We can and [i]will[/i] change them later to be more robust and reflective of reality, but for now this is a fine start:[br][br][list][*]On average, an owl eats 2 mice per day. So each owl represents a reduction in the mouse population of 2 mice per day.[/*][*]On average, each owl reproduces 0.003 times per day (in other words, each owl is the parent of 1 new baby owl per year). So each owl represents an increase in the owl population of 0.003 owls per day.[/*][*]On average, each mouse reproduces 0.055 times per day (in other words, each mouse is the parent of about 20 new baby mice per year). So each mouse represents an increase in the mouse population of 0.055 mice per day.[/*][/list][br]We could add additional assumptions about how these owl and mouse populations impact the rate of change of each other's populations, but this is enough to get started. [br][br]Now let's attempt to build a model of the mouse population with a function [math]M(t)[/math], and a model of the owl population by a function [math]O(t)[/math]. The independent variable, [math]t[/math], is in days in both models. Remember: these are only [i]models. [/i]They are not reality. So don't overthink them.[br][br]The rates of change of these models are governed by the above assumptions. In particular:[br][br][math]\frac{dM}{dt}=0.055M-2O[/math][br][br][math]\frac{dO}{dt}=0.003O[/math][br][i][br]This[/i] is a system of first order differential equations. Its full name is a [b]homogeneous linear system of differential equations with constant coefficients[/b].

There is a beautiful theory that involves the eigenvalues and eigenvectors of matrices that is used to construct solutions to these types of systems of differential equations. We'll learn about it soon. For now, we'll just feast on the results of this theory, and check that it does in fact produce solutions.[br][br]I contend that these two functions are general solutions of the above system of equations:[br][br][math]M\left(t\right)=c_1\left(1\right)e^{\frac{11}{200}t}+c_2\left(-500\right)e^{\frac{3}{1000}t}[/math][br][br][math]O\left(t\right)=c_1\left(0\right)e^{\frac{11}{200}t}+c_2\left(-13\right)e^{\frac{3}{1000}t}[/math][br][br]To check, let's enter these two functions, being sure to agree to create sliders for [code]c_1[/code] and [code]c_2[/code].

Now we need to calculate the first derivatives of [code]M(t)[/code] and [code]O(t)[/code] with [code]Derivative()[/code] and check that they solve the differential equation by verifying that both[br][br][code]Simplify(M'-0.055M+2O)[/code][br][br]and[br][br][code]Simplify(O'-0.003O) [br][/code][br][br]are equal to 0.

[i]Voila[/i]! [br][br]We have verified that these two functions are solutions of the system of two first order differential equations!

Just as we used slope fields to visualize first order differential equations, and systemifications of second order differential equations, we can use slope fields to visualize first order systems of differential equations.[br][br]To visualize the system from above[br][br][math]\frac{dM}{dt}=0.055M-2O[/math][br][br][math]\frac{dO}{dt}=0.003O[/math][br][br]we'll make the decision to think of [math]M[/math] as [code]y[/code] and [math]O[/math] as [code]x[/code], and plot the slope field of [math]\frac{\left(\frac{dM}{dt}\right)}{\left(\frac{dO}{dt}\right)}[/math]. [i]Note[/i]: we could make the opposite decision about how to think of [math]M[/math] and [math]O[/math], but if we do so we will need to remember to respect this decision later when we get to plotting solution functions as a parametric curve.[br][br]With this decision about [math]M[/math] and [math]O[/math], the code to visualize the system is[br][code][br]Slopefield((0.055y-2x)/(0.003x))[/code][br]

Just as in the case of systemifications of second order differential equations, the general solution functions[br][br][math]M\left(t\right)=c_1\left(1\right)e^{\frac{11}{200}t}+c_2\left(-500\right)e^{\frac{3}{1000}t}[/math][br][br][math]O\left(t\right)=c_1\left(0\right)e^{\frac{11}{200}t}+c_2\left(-13\right)e^{\frac{3}{1000}t}[/math][br][br]can be visualized in the slope field as parametric curves. To do so in this case, first declare the functions (agreeing to create sliders for [code]c_1[/code] and [code]c_2[/code]), and then plot the parametric curve of the solution functions in the slope field with[br][br][code]curve(O(t),M(t),t,-100,1000)[br][/code][br]Note that it's essential that [math]O(t)[/math] is the [code]x[/code] position, and [math]M(t)[/math] is in the [code]y[/code] position of [code]curve()[/code] since we made the decision earlier that [math]O[/math] would play the role of [code]x[/code], and [math]M[/math] the role of [code]y[/code]. We could have made a different decision, but we need to be careful to respect that decision here.[br][br]Also, note that I chose [code]t[/code] to go from -100 to 1000 to give us a wider view of the solution functions.

You may be wondering: where's the curve? Well we just need to zoom out a bit. [br][br]Here it is:

Be sure to explore the variety of solutions by adjusting the sliders for [code]c_1[/code] and [code]c_2[/code]. [br][br]Remember: [math]M(t)[/math], mice, is vertical, or [code]y[/code], and [math]O(t)[/math], owls, is horizontal, or [code]x[/code]. Consequently, at the outset, this solution makes no sense since it represents negative mouse [i]and [/i]owl populations. [br][br]Can you find a set of coefficients that generates a model of mice and owls that does make sense? Playing with these coefficients leads us to thinking for a moment about how to work with initial conditions on our functions, and how initial conditions lead to specific solutions of the equation.

One last item to consider in this example are initial conditions on [math]M(t)[/math] and [math]O(t)[/math]. For instance, perhaps we know that at some point in time when we want to start a model there were 13 owls and 499 mice living in Landmark Meadow. In other words [math]M(0)=499[/math] and [math]O(0)=13[/math].[br][br]We can solve a system of [i]algebraic[/i] equations to discover the values of [code]c_1[/code] and [code]c_2[/code] that lead to specific functions [math]M(t)[/math] and [math]O(t)[/math] that both solve the system of [i]differential[/i] equations, and match these initial conditions. [br][br]To this, we need to solve the system of algebraic equations:[br][br][math]M\left(0\right)=c_1\left(1\right)e^{\frac{11}{200}0}+c_2\left(-500\right)e^{\frac{3}{1000}0}=c_1-500c_2=499[/math][br][br][math]O\left(t\right)=c_1\left(0\right)e^{\frac{11}{200}0}+c_2\left(-13\right)e^{\frac{3}{1000}0}=-13c_2=13[/math][br][br]Take a moment to think of a few ways to solve this system. In the end, note that [math]c_2=-1[/math] and [math]c_1=-1[/math].[br][br]Be sure to use the above applet to visualize these constants, the solution in the slope field. What can you say about the future of the mouse population in these initial conditions?

Note that for all selections of coefficients, the curve [code]a[/code] always moves [i]away [/i]from the asymptotic linear axes of the slope field as [code]t[/code] increases. When solutions do this, they are described as [b]unstable[/b] solutions since the populations never stabilize towards a proportional relationship with each other, and instead diverge from a proportional relationship.[br][br]Describing all the different types (stable and unstable) of a system based on the variety of different combinations of the constants is called a [b]phase space analysis [/b](or [b]phase portrait analysis[/b]), and we'll discuss it in detail after we get through a bit of theory of solving systems of differential equations with constant coefficients.