Last time, we learned how to solve systems of differential equations using the elimination method, and we saw that the algebra can be quite involved for even fairly simple systems. Now we will learn a completely different approach using the techniques of linear algebra and matrices, which will also help give us better insight on the questions of long term behaviour of dynamical systems. Let’s start with the general form of a coupled system of linear first order differential equations with constant coefficients, which can be written as x prime equals a x plus b y, and y prime equals c x plus d y. Notice that this sort of resembles the setup for simultaneous algebraic equations, and we can summarize these equations using a column vector for the unknown variables x and y. We can also pack the coefficients a, b, c and d into a two by two matrix. If we call the vector x, this time using a bold or underlined symbol to indicate that this is a vector quantity, we can write down a compact form of this system: x prime equals capital A times x, where capital A is the matrix a b c d. If we were faced with the scalar-valued equivalent of this system, which would be x prime equals a times x, where x and a are just single variable quantities, we know that the solution to this would be of the form x of t equals e to the a times t, times a constant which will be equal to the initial value of x at t equals zero. We’re working with a vector-valued differential equation, and we can still use the same idea, but we will find ourselves trying to evaluate the exponential of a matrix if we try the same solution formula. We actually can do this, because the Taylor series of an exponential function is just a sum of powers. If we sub in the matrix A t, we get a series made up of powers of the matrix A. If A is a diagonalizable matrix, we can write the nth power of A in terms of the eigenvectors U and eigenvalues D as A to the n equals U times D to the n times the inverse of U. Replacing all the powers of A in our series with this leads to the useful fact that the exponential of a matrix A t is equal to U times the exponential of D t times the inverse of U. By analogy, when we solved homogeneous second-order DEs, we can say that the ansatz solutions to our system of DEs will be of the form e to the lambda t times the vector u, where lambda is an eigenvalue of A, found in the matrix D, and vector u is the corresponding eigenvector of A, found in the matrix U. This leads to three different cases depending on the values of the eigenvalues lambda, which are worth knowing for reference. Notice that when matrix A has real eigenvalues, we can use our ansatz solution directly, and when the eigenvalues are complex, the solution is a little more cumbersome, but it can be done using Euler’s formula all the same. Lastly, there is the third rare case of real repeated defective eigenvalues. Defective means that we are unable to find a full set of linearly independent eigenvectors for an eigenvalue, so a diagonalization of A does not exist, and this will turn out to be the equivalent of the repeated roots case: although here, we need to find an additional constant vector v such that A minus lambda i times v equals our one eigenvector u, and then we use the similar trick of multiplying the ansatz by t to get the general solution. Note that if we find that we have repeated eigenvalues but that they still have a full set of eigenvectors, we don’t need to use this method, and we can instead still use the simpler first case. So that’s how to solve the general form of a homogeneous linear system of differential equations: x prime equals A x, where x is a vector and A is a matrix. Now let’s see how to solve the nonhomogeneous case, which will occur if either variable has an explicit dependence on the independent variable t. We can write this in the form x prime equals A times x plus f of t, where f is a vector-valued function. The method that we’ve just gone through will give us the complementary solution, and we need to add on a particular integral, of which we can use the two methods that we learned about: the method of undetermined coefficients and variation of parameters. The method of undetermined coefficients approach is pretty much identical to what we already know, we choose a trial function based on the form of f, and then substitute it in with some multiplier coefficients to be solved for. For example, if we have the following system, where the f of t term contains linear polynomials in t, then our trial function will also be a linear polynomial in t. The only difference is that this time, those unknown coefficients will be vectors instead of scalars. To be thorough and briefly illustrate this example to completion, we can differentiate our particular integral, substitute it into our system of x prime equals Ax plus f of t, multiply out and then equate the like terms in our vector equation, and then solve the linear system of equations to find the resulting constants. This will give us our particular integral. Now remember that this particular integral must be combined with the complementary function in order to get the general solution, so in order to find the complementary function we would find the eigenvalues and eigenvectors of A, which will be as follows. Since both eigenvalues are real we can use the first complementary solution case, which gives us this, and finally we combine that with the particular integral we found just a moment ago, which gives us this for our general solution. If we want to use variation of parameters, we need to collect the independent solutions from our complementary solution into a matrix, capital X, which will act kind of like the Wronskian for our problem, defined like this, where each column of the matrix X is a linearly independent part of the complementary solution. The formula for variation of parameters then looks like this, giving us the particular integral as the matrix X times the integral of the inverse matrix of X times f of t dt, which is a lot easier to remember than the formula from before, but that doesn’t mean it will be any less work, since it requires inverting matrices and taking multiple integrals. Either way, both of these methods allow us to solve nonhomogeneous linear systems, where the substitution method learned in the last tutorial might not have been possible at all, and they also extend to higher order systems too. Let’s now study the concept of system stability. Similar to the concept of a tangent field for representing the trajectories of solutions to first-order differential equations that we learned in the numerical methods tutorial, systems of differential equations that are autonomous, meaning that they do not depend explicitly on t, also have an equivalent notion, that of the phase plane. For systems with two dependent variables x and y, we can have a family of general solutions as parametric equations in t. If we draw the tangents to these curves at every point in the plane, we get the phase plane, which shows us graphically how a system will evolve in time from any given initial state, where the state of the system follows trajectories in the phase plane. For example, here is the phase plane vector field for this nonlinear system of DEs. The regions where dx/dt and dy/dt are zero, meaning the trajectory is stationary in that direction, are shown in red and yellow, and these are called the nullclines of the system. The equilibrium points, also known as fixed points, are the points where all nullclines intersect each other, since it is at those points that the rate of change of the state vector x is zero, so a trajectory that reaches there will never move. It is useful to know that for homogeneous linear systems, the only equilibrium point will be the origin. We can then consider how trajectories around equilibrium points behave, which is the concept of stability. A stable equilibrium point tends to attract nearby trajectories, while an unstable point repels them away, and these can be identified by the eigenvalues of the system, as shown in this table. Negative eigenvalues give exponential decay terms in the solution, while positive eigenvalues give exponentials that grow to infinity. For 2 by 2 linear systems, we can also plot the different types of system based on the trace of the matrix A and the determinant of A, since these uniquely define the eigenvalues of A. This alternative way of representing the possible behaviors of our system is called a Poincaré diagram. When sketching a phase plane, it is often useful to draw the eigenvectors of A, if they are real, as straight lines which are shown in blue on these diagrams, as these will act as asymptotic trajectories that divide the plane up into regions that behave similarly. For nonlinear systems, the nullclines can be curved, which can produce more complicated dynamics that is often hard to investigate directly. However, given any autonomous nonlinear system of the form x prime equals f of x y and y prime equals g of x y, we can perform an operation known as linearization, which is essentially a linear approximation to our system. Given such a system, we can linearize it about a given point x one y one to get an ‘effective’ A matrix given by the Jacobian matrix of our system, which is essentially the gradient operator from vector calculus, but acting on each component separately. The Jacobian is defined as a matrix of all the partial derivatives of our nonlinear functions f and g, evaluated at the point x one y one. This matrix can then be used in place of our nonlinear system for small deviations away from the point x one y one, and we can proceed to assess stability using the eigenvalue theory that we mentioned earlier. The specialized study of nonlinear dynamical systems is quite deep, involving topics such as chaos theory, and this goes beyond the scope of this series, so we won’t go any further in this tutorial. Before we move on, let’s check comprehension.