Chapter 6 Introduction to Systems of ODEs
So far we have discussed ordinary differential equations where the function we have been looking for was a scaler function (y(x):R→R). Where the unknown function is a vector (→y(x):R→Rn), we are dealing with systems of ODEs.
Definition Systems of Ordinary Differential Equations have the following general form:
G1(x,y1,y2,⋯,yn,dy1dx,⋯,dyndx,⋯,dk1y1dxk1,⋯, dknyndxkn)=0,
⋮
Gn(x,y1,y2,⋯,yn,dy1dx,⋯,dyndx,⋯,dk1y1dxk1,⋯, dknyndxkn)=0.
The system is ordinary as we still have one independent variable x, but now in contrast to single ODEs, we have n functions of independent variables y1(x),y2(x),⋯,yn(x) to solve for. This is the implicit form but the systems of ODEs can be written in explicit from as well. Many problems in physics and biology give rise to systems of ODEs. Here are few examples:
These models can be used to predict the dynamics of predator and prey systems such as rabbits (x(t)) and foxes (y(t)). A classic model is Lotka-Volterra model (1925/26) that can exhibit a periodic solution of the population of predator and preys as one goes up and the other goes down.
dxdt=ax−bxy,dydt=dxy−cy.
The chemical rate equations for a set of chemical reactions. For example, consider the reversible binary reaction A+B⇌C with forward rate of k1 and backward rate of k2. We have the following rate equations for the concentrations [A], [B] and [C].
d[A]dt=d[B]dt=k2[C]−k1[A][B],d[C]dt=−d[A]dt=−k2[C]+k1[A][B]. These kind of equations are used in mathematical modelling of biochemical reaction networks in systems and synthetic biology.
This is an example from mechanics. Consider the system of 3 masses and 4 springs that are fixed between two walls. We can write equations of motions that describe the position of the 3 masses as a function of time (x1(t), x2(t) and x3(t)) as seen in Figure 6.1.

Figure 6.1: Diagram of the mass and spring system
Using the second Newton law for each mass (F=ma) and the Hook’s law (F=−kΔx) for the springs and given the relaxed lengths of each spring (l1 to l4), and assuming the distance between the walls is the sum of the the relaxed lengths of the springs, we have the following system of ODEs:
d2x1dt2=−k1(x1−l1)+k2(x2−x1−l2),d2x2dt2=−k2(x2−x1−l2)+k3(x2−x2−l3),d2x3dt2=−k3(x3−x2−l3)+k4(l1+l2+l3−x3).
One of the simplest but influential models of an epidemic is the compartmental SIR model. In this model, disease propagates between the population compartments through susceptible individuals (S) becoming infected (I) after encountering other infected individuals with rate β and finally moving to a recovered population (R) with rate γ. The following simple system of ODEs describes the dynamics of this compartmental model.
dSdt=−βSI,dIdt=βSI−γI,dRdt=γI.
This model and its variations are the basis of a lot of mathematical modelling that has been performed for predicting the effect of different interventions in the Covid-19 world-wide pandemic. If you like to learn a bit more about the SIR ODE system above check out this video. Also, the SIR model can be modeled using a so-called agent-based stochastic approach, where the individuals and their random interactions are specifically followed. You can check out this interesting video from the 3blue1brown series to see some cool exploration of this approach to SIR models.
6.1 Systems of first order ODEs
Systems of ODEs of general order can be rewritten in terms of systems of first order ODEs, so the following system written in explicit form is more general than it seems.
dy1dx=F1(x,y1,y2,⋯,yn), ⋮ dyndx=Fn(x,y1,y2,⋯,yn).
The following second order ODE is the equation of motion for damped harmonic oscillator for a mass m that is attached to an ideal spring which follows the Hook’s law (FS=−kx, where x is position of the mass measured from the spring relaxed position) and has a damping friction force that is propotional and opposite in direction to its velocity FD=−ηdxdt and is acted on by a deriving force F(t).
md2xdt2+ηdxdt+kx=F(t).
By defining new variable u=dx/dt as the velocity of the mass m, we can turn this second order ODE to a system of two first order ODEs:
dxdt=u,dudt=1m[F(t)−ηu−kx].
We can use a general vector notation to write systems of 1st order ODEs as d→yn×1dt=→Fn×1(t,→yn×1).
Here n is the number of equations, t is the independent variable and →y is the function we are looking for. In the next section, we discuss an important subclass of these systems.
6.2 Systems of linear 1st order ODEs with constant coefficients
Systems of linear 1st order ODEs with constant coefficients is an important class that we will discuss their solutions in detail. They have the following general form. dy1dt=n∑i=1α1iyi+g1(t),⋮dyndt=n∑i=1αniyi+gn(t), where αij∈R are the constant coefficients forming the matrix An×n. We can write the system of linear ODEs in matrix form: [dy1dtdy2dt⋮dyndt]=[α11α12…α1nα21α22…α2n⋮⋮⋱⋮αn1αn2…αnn][y1y2⋮yn]+[g1(t)g2(t)⋮gn(t)]. We can write this briefly as d→ydt=A→y+→g(t).
Now if we define the linear operator L[→y]=[ddt−A]→y then we can write the system of linear first order ODEs as L[→y]=→g(t).
Since both matrix A and derivative ddt are linear operators, the operator associated to systems of linear ODEs L is also a linear operator. By linearity we have: L[λ1→y1+λ2→y2]=λ1L[→y1]+λ2L[→y2]. Therefore, the solutions of the homogenous systems of n linear ODEs L[→yH]=0 forms a vector space of dimension n. So a set of linearly independent solutions B={→yi}ni=1 form a basis for this space. Therefore, similar to linear ODEs, the general solution can be written as →yHGS=n∑i=1ci→yi, where cis are n arbitrary constants of integration.
The general solution of non-homogenous systems of 1st order linear ODEs
Similar to the case of linear ODEs, here also we find the general solution in two steps
Obtain complimentary function →yCF by solving the corresponding homogenous systems of ODEs (L[→yCF]=0)
Find a particular integral →yPI that satisfies the full non-homogenous systems of ODEs (L[→yPI]=→g(t)).
Then, for the general solution →yGS we have: →yGS(t;c1,c2,⋯,cn)=→yCF+→yPI.
Solving the homogenous problem
L[→yH]=0⟹d→yHdt=A→yH First, we consider the case where matrix A has n distinct roots and therefore is diagonalizable.
That means that there exists a matrix V where, we have: V−1AV=Λ=[λ1000⋱000λn].
where, ith column of the matrix V is the eigenvector of matrix A corresponding to the eigenvalue λi: A→vi=λi→vi.
We use the eigenvectors matrix V to obtain the solution of the homogenous system. d→ydt=A→y⟹V−1d→ydt=V−1AVV−1→y. Letting →z=V−1→y, we can write the system of ODEs as d→zdt=Λ→z. The ith row of this equation gives us dzidt=λizi, which can be solved, so we obtain: →Z=[c1eλ1tc2eλ2t⋮cneλnt]⟹→yH=V→Z=[⋮⋮⋮→v1→v2⋯→vn⋮⋮⋮][c1eλ1tc2eλ2t⋮cneλnt]. Therefore, we obtain the general solution of the homogenous system of first order linear ODEs to be: →yCF=→yHGS=c1eλ1t→v1+c2eλ2t→v2+⋯+cneλnt→vn. Finding particular integrals
L[→yPI]=→g(t). We will use Ansatz and use the methods of undetermined coefficients and variation of parameters as done for the linear ODEs.
dxdt=−4x−3y−5,dydt=2x+3y−2.
We can write this ODE in vector form as: d→ydt=A→y+→g(t)=[−4−323][xy]+[−5−2].
First step is to obtain →yCF. We obtain the eigenvalues λ1 and λ2 and eigenvectors →v1 and →v2.
λ2+λ−6=0⟹λ1=2,λ2=−3. We have the corresponding eigenvectors: →v1=[1−2],→v2=[3−1].
→yCF=→yHGS(t;c1,c2)=c1e2t[1−2]+c2e−3t[3−1].
2nd step is to find any particular integral →yPI(t)} that satisfies:
L[→yPI]=→g(t)=[−5−2]. We use the ansatz: →yPI=[ab]. By plugging the ansatz into the ODE we obtain the undetermined coefficients a and b: →yPI=[ab]=A−1[52]=−16[33−2−4][52]=[−7/23]. So we have: →yGS(t;c1,c2)=→yCF+→yPI=c1e2t[1−2]+c2e−3t[3−1]+[−7/23].
This is a second order linear ODE and we can solve it using the methods discussed in the last chapter. Using the ansatz eλt, we obtain the following characteristic equation: mλ2+ηλ+k=0⇒λ1,2=−η±√η2−4km2m. If the roots are distinct we obtain: xGS=c1eλ1t+c2eλ2t.
Alternatively, we can transform the second order ODE to systems of first order ODEs as seen before: dxdt=u,dudt=−ηmu−kmx. We obtain the same λ1,2 as above for the eigenvalues of corresponding matrix A for this system of ODEs.
For the eigenvectors (→v1 and →v2) of A, we have:
→v1:A→v1=λ1→v1⟹[01−km−ηm][v1xv1u]=λ1[v1xv1u],
v1u=λ1v1x⟹→v1=[1λ1]. Similarly, we obtain for the second eigenvector: →v2:A→v2=λ2→v2⟹→v2=[1λ2].
This gives us the following general solution: →yGS=[xGSuGS]=c1eλ1t[1λ1]+c2eλ2t[1λ2], which gives us the same general solution for xGS as obtained using the previous method.
When A has repeated eigenvalues
Case 1: A is still diagonalizable (it has n linearly independent eigenvectors). Then we can still use the method described. For example for n=2 we have: A=[λ00λ]⇒→v1=[10],→v2=[01]. So, we have →yCF=c1eλt[10]+c2eλt[01].
Case 2: A is not diagonalizable (it has less than n linearly independent eigenvectors). Then we will use the Jordan normal form. We first discuss an example and then see the general case.Matrix A in this case has repeated roots: (1−λ)(3−λ)+1=0⇒λ1=λ2=2.
We next find the Eigenvector(s):
A→v1=λ1→v1⟹[1−113][v1xv1y]=2[v1xv1y].
Which gives us:
v1x=−v1y⇒→v1=[1−1].
So, this matrix only has one eigenvector and the matrix is not diagonalizable.
We look for similarity transformation to a Jordan normal form (J; an almost diagonal form as defined below). We look for a matrix of the form: W=[1α−1β]. So that: W−1AW=[2102]=J, where, J is the Jordan normal form for a 2 by 2 matrix. We have: Aw=w[2102]⇒[1−113][1α−1β]=[1α−1β][2102]. So, we obtain the following equations for α and β. α−β=1+2αα+3β=−1+2β⟹α+β=−1⟹→w2=[1−2], giving us W=[11−1−2], and we can check that W−1AW=J. The Jordan normal form allows us to solve the non-diagonalizable systems of ODEs: d→ydt=A→y⟹W−1d→ydt=[W−1AW]W−1→y. Letting →z=W−1→y, we obtain d→zdt=J→z, so we have: dz1dt=2z1+z2,dz2dt=2z2. We can solve the second ODE to obtain z2=c2e2t and we can then use this solution to solve for z1 from the first equation above: z1=c1e2t+c2te2t. So in vector form: →zGS=[z1z2]=[c1e2t+c2te2tc2e2t]. Now we can obtain the general solution for →yGS. →yGS=W→zGS=[→v1→w2][c1e2t+c2te2tc2e2t]=(c1e2t+c2te2t)→v1+c2e2t→w2.
The case of non-diagonalizable An×n with one repeated eigenvalue λ
Assume λ is associated with only a single eigenvector. We can use the Jordan normal form (J) to obtain a solution to the associated systems of linear ODEs. We look for a similarity transformation W to transform A to J: W−1AW=J=[λ10000λ10000⋱10000λ10000λ]. Letting →z=W−1→y, we obtain d→zdt=J→z, so we have: dzndt=λzn⇒zn=cneλtdzn−1dt=λzn−1+zn⇒zn−1=cn−1eλt+cnteλtdzn−2dt=λzn−2+zn−1⇒zn−2=cn−2eλt+cn−1teλt+cnt22eλt⋮dz1dt=λz1+z2⇒z1=c1eλt+c2teλt+⋯+cntn−1(n−1)!eλt So, we can obtain →yGS=W→zGS.
There could be situations where the matrix has some distinct eigenvalues and some repeated eigenvalues, which will result in different Jordan normal forms. For example, consider a matrix A3×3 with two distinct eigenvalues one repeated. The suitable Jordan normal form would have the following form: J=[λ1100λ1000λ2].