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):RR). Where the unknown function is a vector (y(x):RRn), 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:

Example 6.1 (Predator-prey systems)

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=axbxy,dydt=dxycy.

Example 6.2 (Chemistry and biochemistry)

The chemical rate equations for a set of chemical reactions. For example, consider the reversible binary reaction A+BC 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.

Example 6.3 (Coupled spring-mass systems)

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.

Diagram of the mass and spring system

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(x1l1)+k2(x2x1l2),d2x2dt2=k2(x2x1l2)+k3(x2x2l3),d2x3dt2=k3(x3x2l3)+k4(l1+l2+l3x3).

Example 6.4 (SIR model of an epidemic)

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).

Example 6.5 (Turning a higher order ODE into systems of 1st order ODEs)

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)ηukx].

We can use a general vector notation to write systems of 1st order ODEs as dyn×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=ni=1α1iyi+g1(t),dyndt=ni=1αniyi+gn(t), where αijR are the constant coefficients forming the matrix An×n. We can write the system of linear ODEs in matrix form: [dy1dtdy2dtdyndt]=[α11α12α1nα21α22α2nαn1αn2αnn][y1y2yn]+[g1(t)g2(t)gn(t)]. We can write this briefly as dydt=Ay+g(t).

Now if we define the linear operator L[y]=[ddtA]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[λ1y1+λ2y2]=λ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=ni=1ciyi, 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

  1. Obtain complimentary function yCF by solving the corresponding homogenous systems of ODEs (L[yCF]=0)

  2. 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]=0dyHdt=AyH 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: V1AV=Λ=[λ1000000λn].

where, ith column of the matrix V is the eigenvector of matrix A corresponding to the eigenvalue λi: Avi=λivi.

We use the eigenvectors matrix V to obtain the solution of the homogenous system. dydt=AyV1dydt=V1AVV1y. Letting z=V1y, we can write the system of ODEs as dzdt=Λz. The ith row of this equation gives us dzidt=λizi, which can be solved, so we obtain: Z=[c1eλ1tc2eλ2tcneλnt]yH=VZ=[v1v2vn][c1eλ1tc2eλ2tcneλnt]. Therefore, we obtain the general solution of the homogenous system of first order linear ODEs to be: yCF=yHGS=c1eλ1tv1+c2eλ2tv2++cneλntvn. 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.

Example 6.6 Solve the system of ODEs for {x(t),y(t)}.

dxdt=4x3y5,dydt=2x+3y2.

We can write this ODE in vector form as: dydt=Ay+g(t)=[4323][xy]+[52].

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=[12],v2=[31].

yCF=yHGS(t;c1,c2)=c1e2t[12]+c2e3t[31].

2nd step is to find any particular integral yPI(t)} that satisfies:

L[yPI]=g(t)=[52]. We use the ansatz: yPI=[ab]. By plugging the ansatz into the ODE we obtain the undetermined coefficients a and b: yPI=[ab]=A1[52]=16[3324][52]=[7/23]. So we have: yGS(t;c1,c2)=yCF+yPI=c1e2t[12]+c2e3t[31]+[7/23].

Example 6.7 Solve the problem of damped Harmonic spring with zero forcing md2xdt2+ηdxdt+kx=0

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=η±η24km2m. 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=ηmukmx. 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:Av1=λ1v1[01kmηm][v1xv1u]=λ1[v1xv1u],

v1u=λ1v1xv1=[1λ1]. Similarly, we obtain for the second eigenvector: v2:Av2=λ2v2v2=[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.
Example 6.8 dydt=Ay;A=[1113].

Matrix A in this case has repeated roots: (1λ)(3λ)+1=0λ1=λ2=2.

We next find the Eigenvector(s): Av1=λ1v1[1113][v1xv1y]=2[v1xv1y].
Which gives us: v1x=v1yv1=[11]. 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: W1AW=[2102]=J, where, J is the Jordan normal form for a 2 by 2 matrix. We have: Aw=w[2102][1113][1α1β]=[1α1β][2102]. So, we obtain the following equations for α and β. αβ=1+2αα+3β=1+2βα+β=1w2=[12], giving us W=[1112], and we can check that W1AW=J. The Jordan normal form allows us to solve the non-diagonalizable systems of ODEs: dydt=AyW1dydt=[W1AW]W1y. Letting z=W1y, we obtain dzdt=Jz, 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=WzGS=[v1w2][c1e2t+c2te2tc2e2t]=(c1e2t+c2te2t)v1+c2e2tw2.

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: W1AW=J=[λ10000λ1000010000λ10000λ]. Letting z=W1y, we obtain dzdt=Jz, so we have: dzndt=λznzn=cneλtdzn1dt=λzn1+znzn1=cn1eλt+cnteλtdzn2dt=λzn2+zn1zn2=cn2eλt+cn1teλt+cnt22eλtdz1dt=λz1+z2z1=c1eλt+c2teλt++cntn1(n1)!eλt So, we can obtain yGS=WzGS.

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].