3.1 Introduction to Systems of ODEs

Often we do not have just one dependent variable and just one differential equation, we may end up with systems of several equations and several dependent variables even if we start with a single equation.

If we have several dependent variables, suppose $y sub 1$ , $y sub 2$ ,..., $y_{n}$ , then we can have a differential equation involving all of them and their derivatives. For example, $y 1 ″ equals f open parenthesis y 1 ′, y 2 ′, y sub 1, y sub 2, x close parenthesis$ . Usually, when we have two dependent variables we have two equations such as

$\begin{matrix} y_{1}^{″} = f_{1} (y_{1}^{'}, y_{2}^{'}, y_{1}, y_{2}, x) \\ y_{2}^{″} = f_{2} (y_{1}^{'}, y_{2}^{'}, y_{1}, y_{2}, x) \end{matrix}$

for some functions $f sub 1$ and $f sub 2$ . We call the above a system of differential equations. More precisely, the above is a second order system of ODEs as second order derivatives appear. The system

$\begin{matrix} x_{1}^{'} = g_{1} (x_{1}, x_{2}, x_{3}, t), \\ x_{2}^{'} = g_{2} (x_{1}, x_{2}, x_{3}, t), \\ x_{3}^{'} = g_{3} (x_{1}, x_{2}, x_{3}, t), \end{matrix}$

is a first order system, where

x sub 1, x sub 2, x sub 3

are the dependent variables, and

t

is the independent variable.

The terminology for systems is essentially the same as for single equations. For the system above, a solution is a set of three functions $x sub 1 open parenthesis t close parenthesis$ , $x sub 2 open parenthesis t close parenthesis$ , $x sub 3 open parenthesis t close parenthesis$ , such that

$\begin{matrix} x_{1}^{'} (t) = g_{1} (x_{1} (t), x_{2} (t), x_{3} (t), t), \\ x_{2}^{'} (t) = g_{2} (x_{1} (t), x_{2} (t), x_{3} (t), t), \\ x_{3}^{'} (t) = g_{3} (x_{1} (t), x_{2} (t), x_{3} (t), t) . \end{matrix}$

We usually also have an initial condition. Just like for single equations we specify $x sub 1$ , $x sub 2$ , and $x sub 3$ for some fixed $t$ . For example, $x sub 1 open parenthesis 0 close parenthesis equals a sub 1$ , $x sub 2 open parenthesis 0 close parenthesis equals a sub 2$ , $x sub 3 open parenthesis 0 close parenthesis equals a sub 3$ . For some constants $a sub 1$ , $a sub 2$ , and $a sub 3$ . For the second order system we would also specify the first derivatives at a point. And if we find a solution with constants in it, where by solving for the constants we find a solution for any initial condition, we call this solution the general solution. Best to look at a simple example.

Example 1

Sometimes a system is easy to solve by solving for one variable and then for the second variable. Take the first order system

$\begin{matrix} y_{1}^{'} = y_{1}, \\ y_{2}^{'} = y_{1} - y_{2}, \end{matrix}$

with initial conditions of the form $y sub 1 open parenthesis 0 close parenthesis equals 1$ and $y sub 2 open parenthesis 0 close parenthesis equals 2$ .

Show solution

We note that $y_{1} = C_{1} e^{x}$ is the general solution of the first equation. We then plug this $y sub 1$ into the second equation and get the equation $y_{2}^{'} = C_{1} e^{x} - y_{2}$ , which is a linear first order equation that is easily solved for $y sub 2$ . By the method of integrating factor we obtain

$e^{x} y_{2} = \frac{C_{1}}{2} e^{2 x} + C_{2}$

$y_{2} = \frac{C_{1}}{2} e^{x} + C_{2} e^{- x} .$

The general solution to the system is, therefore,

$y_{1} = C_{1} e^{x}, and y_{2} = \frac{C_{1}}{2} e^{x} + C_{2} e^{- x} .$

We now solve for $C sub 1$ and $C sub 2$ given the initial conditions. We substitute $x equals 0$ and find that $C sub 1 equals 1$ and $C_{2} = \frac{3}{2}$ . Thus the solution is: $y_{1} = e^{x},$ and $y_{2} = \frac{1}{2} e^{x} + \frac{3}{2} e^{- x}$ .

Generally, we will not be so lucky to be able to solve for each variable separately as in the example above, and we will have to solve for all variables at once. While we won’t generally be able to solve for one variable and then the next, we will try to salvage as much as possible from this technique. It will turn out that in a certain sense we will still (try to) solve a bunch of single equations and put their solutions together. Let’s not worry right now about how to solve systems yet.

We will mostly consider the linear systems. The example above is a so-called linear first order system. It is linear as none of the dependent variables or their derivatives appear in nonlinear functions or with powers higher than one ( $x$ , $y$ , $x^{'}$ and $y^{'}$ , constants, and functions of $t$ can appear, but not $x y$ or ${(y^{'})}^{2}$ or $x^{3}$ ). Another, more complicated, example of a linear system is

$\begin{matrix} {y_{1}^{'}}^{'} = e^{t} y_{1}^{'} + t^{2} y_{1} + 5 y_{2} + \sin (t), \\ {y_{2}^{'}}^{'} = t y_{1}^{'} - y_{2}^{'} + 2 y_{1} + \cos (t) . \end{matrix}$

Applications

Let us consider some simple applications of systems and how to set up the equations.

Example 2

First, we consider salt and brine tanks, but this time water flows from one to the other and back. We again consider that the tanks are evenly mixed.

Figure

1

: A closed system of two brine tanks.

Suppose we have two tanks, each containing volume $V$ liters of salt brine. The amount of salt in the first tank is $x sub 1$ grams, and the amount of salt in the second tank is $x sub 2$ grams. The liquid is perfectly mixed and flows at the rate $r$ liters per second out of each tank into the other. See Figure $1$ .

The rate of change of $x sub 1$ , that is $x 1 ′$ , is the rate of salt coming in minus the rate going out. The rate coming in is the density of the salt in tank 2, that is $\frac{x_{2}}{V}$ , times the rate $r$ . The rate coming out is the density of the salt in tank 1, that is $\frac{x_{1}}{V}$ , times the rate $r$ . In other words it is

$x_{1}^{'} = \frac{x_{2}}{V} r - \frac{x_{1}}{V} r = \frac{r}{V} x_{2} - \frac{r}{V} x_{1} = \frac{r}{V} (x_{2} - x_{1}) .$

Similarly we find the rate

x 2 ′

, where the roles of

x sub 1

and

x sub 2

are reversed. All in all, the system of ODEs for this problem is

$\begin{matrix} x_{1}^{'} = \frac{r}{V} (x_{2} - x_{1}), \\ x_{2}^{'} = \frac{r}{V} (x_{1} - x_{2}) . \end{matrix}$

In this system we cannot solve for

x sub 1

x sub 2

separately. We must solve for both

x sub 1

and

x sub 2

at once, which is intuitively clear since the amount of salt in one tank affects the amount in the other. We can’t know

x sub 1

before we know

x sub 2

, and vice versa.

We don’t yet know how to find all the solutions, but intuitively we can at least find some solutions. Suppose we know that initially the tanks have the same amount of salt. That is, we have an initial condition such as $x_{1} (0) = x_{2} (0) = C$ . Then clearly the amount of salt coming and out of each tank is the same, so the amounts are not changing. In other words, $x sub 1 equals C$ and $x sub 2 equals C$ (the constant functions) is a solution: $x_{1}^{'} = x_{2}^{'} = 0$ , and $x_{2} - x_{1} = x_{1} - x_{2} = 0$ , so the equations are satisfied.

Let us think about the setup a little bit more without solving it. Suppose the initial conditions are $x sub 1 open parenthesis 0 close parenthesis equals A$ and $x sub 2 open parenthesis 0 close parenthesis equals B$ , for two different constants $A$ and $B$ . Since no salt is coming in or out of this closed system, the total amount of salt is constant. That is, $x sub 1 plus x sub 2$ is constant, and so it equals $A + B$ . Intuitively if $A$ is bigger than $B$ , then more salt will flow out of tank one than into it. Eventually, after a long time we would then expect the amount of salt in each tank to equalize. In other words, the solutions of both $x sub 1$ and $x sub 2$ should tend towards $\frac{A + B}{2}$ . Once you know how to solve systems you will find out that this really is so.

Example 3

As an example application, let us think of mass and spring systems again.

Figure

2

As an example application, let us think of mass and spring systems again. Suppose we have one spring with constant $k$ , but two masses $m sub 1$ and $m sub 2$ . We can think of the masses as carts, and we will suppose that they ride along a straight track with no friction. Let $x sub 1$ be the displacement of the first cart and $x sub 2$ be the displacement of the second cart. That is, we put the two carts somewhere with no tension on the spring, and we mark the position of the first and second cart and call those the zero positions. Then $x sub 1$ measures how far the first cart is from its zero position, and $x sub 2$ measures how far the second cart is from its zero position. The force exerted by the spring on the first cart is $k open parenthesis x sub 2 minus x sub 1 close parenthesis$ , since $x sub 2 minus x sub 1$ is how far the string is stretched (or compressed) from the rest position. The force exerted on the second cart is the opposite, thus the same thing with a negative sign.

Newton’s second law states that force equals mass times acceleration. So the system of equations governing the setup is

$\begin{matrix} m_{1} x_{1}^{″} = k (x_{2} - x_{1}) \\ m_{2} x_{2}^{″} = - k (x_{2} - x_{1}) \end{matrix}$

In this system we cannot solve for the $x sub 1$ or $x sub 2$ variable separately. That we must solve for both $x sub 1$ and $x sub 2$ at once is intuitively clear, since where the first cart goes depends exactly on where the second cart goes and vice-versa.

Changing to First Order

Before we talk about how to handle systems, let us note that in some sense we need only consider first order systems. Let us take an $n^{t h}$ order differential equation

$y^{(n)} = F (y^{(n - 1)}, . . ., y^{'}, y, x) .$

We define new variables $u sub 1,..., u sub n$ and write the system

$\begin{matrix} u_{1}^{'} = u_{2} \\ u_{2}^{'} = u_{3} \\ ⋮ \\ u_{n - 1}^{'} = u_{n} \\ u_{n}^{'} = F (u_{n}, u_{n - 1}, \dots, u_{2}, u_{1}, x) \end{matrix}$

We solve this system for $u sub 1, u sub 2, and so on, u sub n$ . Once we have solved for the $u$ ’s, we can discard $u sub 2$ through $u_{n}$ and let $y equals u sub 1$ . We note that this $y$ solves the original equation.

A similar process can be followed for a system of higher order differential equations. For example, a system of $k$ differential equations in $k$ unknowns, all of order $n$ , can be transformed into a first order system of $n \times k$ equations and $n \times k$ unknowns.

Example 6

Sometimes we can use this idea in reverse as well. Let us take the system

$x^{'} = 2 y - x, y^{'} = x,$

where the independent variable is $t$ . We wish to solve for the initial conditions $x open parenthesis 0 close parenthesis equals 1$ and $y open parenthesis 0 close parenthesis equals 0$ .

If we differentiate the second equation we get $y^{″} = x^{'}$ . We know what $x^{'}$ is in terms of $x$ and $y$ , and we know that $x = y^{'}$ .

$y^{″} = x^{'} = 2 y - x = 2 y - y^{'} .$

We now have the equation $y^{″} + y^{'} - 2 y = 0$ . We know how to solve this equation and we find that $y = C_{1} e^{- 2 t} + C_{2} e^{t}$ . Once we have $y$ we use the equation $y^{'} = x$ to get $x$ .

$x = y^{'} = - 2 C_{1} e^{- 2 t} + C_{2} e^{t}$

We solve for the initial conditions $1 = x (0) = - 2 C_{1} + C_{2}$ and $0 = y (0) = C_{1} + C_{2}$ . Hence, $C_{1} = - C_{2}$ and $1 equals 3 C sub 2$ . So $C_{1} = - \frac{1}{3}$ and $C_{2} = \frac{1}{3}$ . Our solution is

$x = \frac{2 e^{- 2 t} + e^{t}}{3}, y = \frac{- e^{- 2 t} + e^{t}}{3}$

It is useful to go back and forth between systems and higher order equations for other reasons. For example, software for solving ODE numerically (approximation) is generally for first order systems. To use it, you take whatever ODE you want to solve and convert it to a first order system. It is not very hard to adapt computer code for the Euler or Runge–Kutta method for first order equations to handle first order systems. We simply treat the dependent variable not as a number but as a vector. In many mathematical computer languages there is almost no distinction in syntax.

Autonomous Systems and Vector Fields

A system where the equations do not depend on the independent variable is called an autonomous system. For example the system $x^{'} = 2 y - x$ , $y^{'} = x$ is autonomous as $t$ is the independent variable but does not appear in the equations.

For autonomous systems we can draw the so-called direction field or vector field, a plot similar to a slope field, but instead of giving a slope at each point, we give a direction (and a magnitude). The previous example, $x^{'} = 2 y - x$ , $y^{'} = x$ , says that at the point $(x, y)$ the direction in which we should travel to satisfy the equations should be the direction of the vector $open parenthesis 2 y minus x, x close parenthesis$ with the speed equal to the magnitude of this vector. So we draw the vector $open parenthesis 2 y minus x, x close parenthesis$ at the point $(x, y)$ and we do this for many points on the $x y$ -plane. For example, at the point $open parenthesis 1, 2 close parenthesis$ we draw the vector $(2 (2) - 1, 1) = (3, 1)$ , a vector pointing to the right and a little bit up, while at the point $open parenthesis 2, 1 close parenthesis$ we draw the vector $(2 (1) - 2, 2) = (0, 2)$ a vector that points straight up. When drawing the vectors, we will scale down their size to fit many of them on the same direction field. If we drew the arrows at the actual size, the diagram would be a jumbled mess once you would draw more than a couple of arrows. So we scale them all so that not even the longest one interferes with the others. We are mostly interested in their direction and relative size. See Figure $3$ .

Function graph showing the direction field of the system x′ = k*y - x, y′ = x, drawn from 308 sample arrows, and 9 trajectories through (1, 0), (1, 1), (-1, -1), (2, -1), (-2, 1), (-1, 1), (1, -1), (0, 2) and (0, -2). Adjustable parameter: Coefficient of y in x' (Lebl's system is k = 2) (k) = 2. Viewing window: x from -4.93 to 5.93, y from -2.86 to 3.86. — This section's own direction field, x′ = 2y − x, y′ = x, carrying nine trajectories instead of the single one Lebl's Figure 4 shows. The curve seeded at (1, 0) is his: follow it two units of t and it arrives at (2.475, 2.457), exactly the value the text computes. Four seeds sit on the two directions the section's own y″ + y′ − 2y = 0 hands you — the e^t part of the solution runs out along (1, 1), the e^(-2t) part runs in along (2, −1) — so those four draw as perfectly straight rays, the skeleton one trajectory cannot show. Now drag k, which stands in for Lebl's 2: the characteristic equation becomes λ² + λ − k = 0, so crossing k = 0 destroys the saddle outright (the whole y-axis turns into equilibria) and leaves a stable node, and below k = −0.25 the discriminant 1 + 4k goes negative and that node becomes a spiral sink — three regimes in one drag.

We can draw a path of the solution in the plane. Suppose the solution is given by $x = f (t)$ , $y = g (t)$ . We pick an interval of $t$ (say $0 less than or equal to t less than or equal to 2$ for our example) and plot all the points $(f (t), g (t))$ for $t$ in the selected range. The resulting picture is called the phase portrait (or phase plane portrait). The particular curve obtained is called the trajectory or solution curve. See an example plot in Figure $4$ . In the figure the solution starts at $open parenthesis 1, 0 close parenthesis$ and travels along the vector field for a distance of 2 units of $t$ . We solved this system precisely, so we compute $x open parenthesis 2 close parenthesis$ and $y open parenthesis 2 close parenthesis$ to find $x open parenthesis 2 close parenthesis approximately equals 2.475$ and $y open parenthesis 2 close parenthesis approximately equals 2.457$ . This point corresponds to the top right end of the plotted solution curve in the figure.

Notice the similarity to the diagrams we drew for autonomous systems in one dimension. But note how much more complicated things become when we allow just one extra dimension.

We can draw phase portraits and trajectories in the $x y$ -plane even if the system is not autonomous. In this case, however, we cannot draw the direction field, since the field changes as $t$ changes. For each $t$ we would get a different direction field.

Direction field for x' = 2y-x and y' = x — Figure $3$ : The direction field for $x^{'} = 2 y - x, y^{'} = x$ .

The direction field for x' = 2y - x, y' = x with the trajectory of the solution starting at (1, 0) for 0 <= t <= 2. — Figure $4$ : The direction field for $x^{'} = 2 y - x, y^{'} = x$ with the trajectory of the solution starting at $open parenthesis 1, 0 close parenthesis$ for $0 less than or equal to t less than or equal to 2$ .

Picard’s theorem

Perhaps before going further, let us mention that Picard’s theorem on existence and uniqueness still holds for systems of ODE. Let us restate this theorem in the setting of systems. A general first order system is of the form

$\begin{matrix} x_{1}^{'} = F_{1} (x_{1}, x_{2}, \dots, x_{n}, t), \\ x_{2}^{'} = F_{2} (x_{1}, x_{2}, \dots, x_{n}, t), \\ ⋮ \\ x_{n}^{'} = F_{n} (x_{1}, x_{2}, \dots, x_{n}, t) . \end{matrix}$

(3.1.1)

That is, a unique solution exists for any initial condition given that the system is reasonable ( $F_{j}$ and its partial derivatives in the $x$ variables are continuous). As for single equations we may not have a solution for all time $t$ , but at least for some short period of time.

As we can change any $n$ th order ODE into a first order system, then we notice that this theorem provides also the existence and uniqueness of solutions for higher order equations that we have until now not stated explicitly.

Adapted from Differential Equations for Engineers by Jiří Lebl (Oklahoma State University), hosted on LibreTexts (math.libretexts.org) and licensed under CC BY-SA 4.0. Changes were made. License: CC-BY-SA-4.0.