Systems of First-Order Equations

systems.mws

Douglas B. Meade (meade@math.sc.edu)

1 Overview

First-order linear systems of ordinary differential equations (and initial value problems) can be solved in a number of different ways. The most efficient technique depends on your background in linear algebra, and the specific problem being studied. The general idea, as introduced in Section 9.3, is to write the system in normal form, then find a fundamental matrix for the homogeneous system and, if the system is non-homogeneous, a particular solution to the nonhomogeneous system. This worksheet contains one example, solved using three different approaches to determine the fundamental matrix.

Maple’s linear algebra package, linalg, will be used extensively to perform a variety of matrix operations, including matrix multiplication, inverse of a matrix, eigenvalue decomposition, and matrix exponentials.

2 Objectives

• to demonstrate the critical role that the fundamental matrix plays in the solution of linear systems of ODEs

• to present multiple methods for the (theoretical) computation of the fundamental matrix

• to illustrate techniques that can be used to realize the different methods for computing the fundamental matrix

3 New Maple Commands

• diag returns the diagonal matrix with specified diagonal entries

• det computes the determinant of a matrix • diag creates a diagonal matrix

• eigenvects determines the eigenvalues (with multiplicity) and, eigenvectors for a square matrix

• evalm evaluates a given matrix expression

• exponential the matrix exponential

• inverse computes the inverse of a matrix

• linalg Maple’s linear algebra package

• map applies an or function to all elements in a set, list, equation, or matrix (see also, map2)

• matrix creates a matrix data structure of the given size

• multiply matrix-matrix or matrix-vector multiplication

• nullspace determines a basis for the nullspace of a given matrix

• transpose returns the transpose of a matrix

• &* non-commutative matrix multiplication (replaces * for matrices)

4 Example 1 ( Nagle/Saff, p. 550 )

> restart;
> with( plots ):
> with( DEtools ):

Note that, in addition to the standard initialization, the linalg package is also loaded for future use.

> with( linalg );
Warning, new definition for norm

Warning, new definition for trace

[BlockDiagonal, GramSchmidt, JordanBlock, LUdecomp, QRdecomp, Wronskian, addcol, addrow, adj , adjoint, angle, augment, backsub, band, basis, bezout, blockmatrix, charmat, charpoly, cholesky, col, coldim, colspace, colspan, companion, concat, cond, copyinto, crossprod, curl, definite, delcols, delrows, det, diag, diverge, dotprod, eigenvals, eigenvalues, eigenvectors, eigenvects, entermatrix, equal, exponential, extend, ffgausselim, fibonacci, forwardsub, frobenius, gausselim, gaussjord, geneqns, genmatrix, grad, hadamard, hermite, hessian, hilbert, htranspose, ihermite, indexfunc, innerprod, intbasis, inverse, ismith, issimilar, iszero, jacobian, jordan, kernel, laplacian, leastsqrs, linsolve, matadd, matrix, minor, minpoly, mulcol, mulrow, multiply, norm, normalize, nullspace, orthog, permanent, pivot, potential, randmatrix, randvector, rank, ratform, row, rowdim, rowspace, rowspan, rref , scalarmul, singularvals, smith, stack, submatrix, subvector, sumbasis, swapcol, swaprow, sylvester, toeplitz, trace, transpose, vandermonde, vecpotent, vectdim, vector, wronskian]

4.1 Explanation of Warning Messages

The two warning messages indicate that the norm and trace commands now have new definitions. (To see the on-line help for the old and new definitions, see ?norm and ?linalg,norm and ?trace and ?linalg,trace.)

Find a general solution to the first-order linear system in normal form ( x’ = Ax + f ) for the special case where

4.2 x is a three-dimensional (column) vector

> x := matrix( 3, 1 );
> evalm( x );

x := array(1..3, 1..1, [])

[ \begin{matrix}x_{1, 1} \\ x_{2, 1} \\ x_{3, 1} \end{matrix} ]

4.3 A is the 3 x 3 matrix

> A := matrix( 3, 3, [[1, -2, 2], [-2, 1, 2], [2, 2, 1]] );

A := [ \begin{matrix}1 & −2 & 2 \\ −2 & 1 & 2 \\ 2 & 2 & 1 \end{matrix} ]

4.4 and, f = t g is the three-dimensiona (column) vector

> f := t -> evalm( t * g );

f := t → evalm(t g)

4.4.1 Note concerning Maple’s evaluation of matrices ( evalm )

Note that Maple does not automatically “evaluate” matrices and matrix operations. This is done for a number of reasons, including the promotion of high-level operations and avoidance of huge amounts of useless output.

To demonstrate the need for the evalm in the definition of f, consider the following three commands:

> t*x;

t x

> t*evalm(x);

t [ \begin{matrix}x_{1, 1} \\ x_{2, 1} \\ x_{3, 1} \end{matrix} ]

> evalm( t*x );

[ \begin{matrix}t x_{1, 1} \\ t x_{2, 1} \\ t x_{3, 1} \end{matrix} ]

Evaluation of matrices within expressions is controlled, by the user, through the use of the evalm command (see the on-line help for evalm for further discussion of this issue).

where g is the (column) vector

> g := matrix( 3, 1, [ -9, 0, -18 ] );

g := [ \begin{matrix}−9 \\ 0 \\ −18 \end{matrix} ]

Just to check that all is as it should be,

   > f(t);

[ \begin{matrix}−9 t \\ 0 \\ −18 t \end{matrix} ]

5 Solution by Eigenvalue Decomposition and Variation of Parameters

The first approach to this problem is most similar to the one presented in Nagle and Saff. The principle difference is the use of variation of parameters in the place of the method of undetermined coefficients. The key ingredient needed for this approach is a fundamental matrix, X(t), for the homogeneous system. One formula for a fundamental matrix is: X(t) = T e^{(Λ t)} where Λ is the diagonal matrix containing the eigenvalues of A and T is the corresponding matrix containing the eigenvectors of A. (That is, A = T Λ T^{(−1)}.)

5.1 Comment on computation of eigenvalue decompositions in Maple

There are a number of ways to find the eigenvalue decomposition of a matrix in Maple. The approach used here most closely follows the hand calculations. The eigenvects command in the linalg package will usually be more efficient (see ?linalg,eigenvects for full details).

Similar high-level solutions to other parts of this problem also exist, e.g., using &*() for the identity matrix, charpoly for the characteristic polynomial, exponential for the matrix exponential. Descriptions of these, and other, commands can be found in the on-line help for the linalg package.

The characteristic matrix, A − λ I, and characteristic polynomial, det(A − λ I), are

   > CHARm := evalm( A - lambda*diag(1,1,1) );

CHARm := [ \begin{matrix}1 − λ & −2 & 2 \\ −2 & 1 − λ & 2 \\ 2 & 2 & 1 − λ \end{matrix} ]

> CHARp := det(CHARm);

CHARp := −27 + 9 λ + 3 λ^{2} − λ^{3}

The eigenvalues of A are the roots of the characteristic polynomial

> EIGvals := solve( CHARp=0, lambda );

EIGvals := −3, 3, 3

Thus, the diagonal matrix of eigenvalues is

> Lambda := diag( EIGvals );

Λ := [ \begin{matrix}−3 & 0 & 0 \\ 0 & 3 & 0 \\ 0 & 0 & 3 \end{matrix} ]

The eigenspace for λ = −3 is the null space of the characteristic matrix with λ = −3, that is,

> EIGvects[-3] := nullspace( subs( lambda=-3, evalm(CHARm) ) );

EIGvects_{−3} := {[1, 1, −1]}

For λ = 3 the eigenspace could be two-dimensional:

> EIGvects[3] := nullspace( subs( lambda=3, evalm(CHARm) ) );

EIGvects_{3} := {[1, 0, 1], [−1, 1, 0]}

The eigenvector matrix is, therefore,

> T := transpose( matrix( 3, 3, [ op(EIGvects[-3]), op(EIGvects[3])
] ) );

T := [ \begin{matrix}1 & −1 & 1 \\ 1 & 1 & 0 \\ −1 & 0 & 1 \end{matrix} ]

Just to check, let’s see that A = T Λ T^{{−1}}.

> evalm( A = T &* Lambda &* inverse(T) );

[ \begin{matrix}1 & −2 & 2 \\ −2 & 1 & 2 \\ 2 & 2 & 1 \end{matrix} ] = [ \begin{matrix}1 & −2 & 2 \\ −2 & 1 & 2 \\ 2 & 2 & 1 \end{matrix} ]

Looks good!

5.2 A high-level approach to the eigenvalue decomposition – eigenvects

Note that the eigenvalues and eigenvectors can be obtained in one step with the use of the eigenvects command from the linalg package.

> EIGSYS := eigenvects( A );

EIGSYS := [3, 2, {[1, −1, 0], [0, 1, 1]}], [−3, 1, {[−1, −1, 1]}]

It can be a challenge to extract the eigenvalues and eigenvectors from the list returned by eigenvects. (Try it!) See ?linalg,eigenvects for a full description of the meaning of each element of this list.

A fundamental matrix for the homogeneous differential equation can now be constructed

> X := t -> evalm( T &* diag( exp(-3*t), exp(3*t), exp(3*t) ) ):
> evalm( X(t) );

[ \begin{matrix}e^{(−3 t)} & −e^{(3 t)} & e^{(3 t)} \\ e^{(−3 t)} & e^{(3 t)} & 0 \\ −e^{(−3 t)} & 0 & e^{(3 t)} \end{matrix} ]

5.3 Verification of fundamental matrix

Prior to finding a particular equation, by variation of parameters, let’s check that X actually solves the homogeneous equation

> CHECK1 := evalm( map( diff, X(t), t ) - A &* X(t) );

CHECK1 := [ \begin{matrix}0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{matrix} ]

The general solution to the homogeneous equation is X(t) c for any threedimensional vector c. That is,

> C := matrix( 3, 1 ):
> evalm( C );

[ \begin{matrix}C_{1, 1} \\ C_{2, 1} \\ C_{3, 1} \end{matrix} ]

> xHOMO := evalm( X(t) &* C );

xHOMO := [ \begin{matrix}e^{(−3 t)} C_{1, 1} − e^{(3 t)} C_{2, 1} + e^{(3 t)} C_{3, 1} \\ e^{(−3 t)} C_{1, 1} + e^{(3 t)} C_{2, 1} \\ −e^{(−3 t)} C_{1, 1} + e^{(3 t)} C_{3, 1} \end{matrix} ]

A particular solution is, thus, given by x_{p}(t) = X(t) ∫_{0}^{t} X(s)^{(−1)} f(s) ds:

> xPART := evalm( X(t) &* map( int, evalm( inverse( X(s) ) &* f(s)
), s=0..t ) );

xPART := [ \begin{matrix}%1 − e^{(3 t)} ((t + \frac{1}{3}) e^{(−3 t)} − \frac{1}{3}) + e^{(3 t)} ((5 t + \frac{5}{3}) e^{(−3 t)} − \frac{5}{3}) \\ %1 + e^{(3 t)} ((t + \frac{1}{3}) e^{(−3 t)} − \frac{1}{3}) \\ −%1 + e^{(3 t)} ((5 t + \frac{5}{3}) e^{(−3 t)} − \frac{5}{3}) \end{matrix} ]

%1 := e^{(−3 t)} (e^{(3 t)} t − \frac{1}{3} e^{(3 t)} + \frac{1}{3})

Further simplification occurs by normalizing each term, then combining the powers:

> xPART := map( normal, evalm(xPART) );

xPART := [ \begin{matrix}5 e^{(−3 t)} e^{(3 t)} t + e^{(−3 t)} e^{(3 t)} + \frac{1}{3} e^{(−3 t)} − \frac{4}{3} e^{(3 t)} \\ 2 e^{(−3 t)} e^{(3 t)} t + \frac{1}{3} e^{(−3 t)} − \frac{1}{3} e^{(3 t)} \\ 4 e^{(−3 t)} e^{(3 t)} t + 2 e^{(−3 t)} e^{(3 t)} − \frac{1}{3} e^{(−3 t)} − \frac{5}{3} e^{(3 t)} \end{matrix} ]

> xPART := map( combine, xPART, power );

xPART := [ \begin{matrix}5 t + 1 + \frac{1}{3} e^{(−3 t)} − \frac{4}{3} e^{(3 t)} \\ 2 t + \frac{1}{3} e^{(−3 t)} − \frac{1}{3} e^{(3 t)} \\ 4 t + 2 − \frac{1}{3} e^{(−3 t)} − \frac{5}{3} e^{(3 t)} \end{matrix} ]

Therefore, the general solution to the nonhomogeneous equation is

> x := evalm( xHOMO + xPART );

x := [ \begin{matrix}e^{(−3 t)} C_{1, 1} − e^{(3 t)} C_{2, 1} + e^{(3 t)} C_{3, 1} + 5 t + 1 + \frac{1}{3} e^{(−3 t)} − \frac{4}{3} e^{(3 t)} \\ e^{(−3 t)} C_{1, 1} + e^{(3 t)} C_{2, 1} + 2 t + \frac{1}{3} e^{(−3 t)} − \frac{1}{3} e^{(3 t)} \\ −e^{(−3 t)} C_{1, 1} + e^{(3 t)} C_{3, 1} + 4 t + 2 − \frac{1}{3} e^{(−3 t)} − \frac{5}{3} e^{(3 t)} \end{matrix} ]

The terms in this solution can be re-written in a variety of ways. Two of the most effective are to group the terms according to the type of exponential term that appears in the solution

> x := map( collect, x, {exp(-3*t), exp(3*t)} );

x := [ \begin{matrix}(\frac{1}{3} + C_{1, 1}) e^{(−3 t)} + (−C_{2, 1} + C_{3, 1} − \frac{4}{3}) e^{(3 t)} + 1 + 5 t \\ (\frac{1}{3} + C_{1, 1}) e^{(−3 t)} + (−\frac{1}{3} + C_{2, 1}) e^{(3 t)} + 2 t \\ (−C_{1, 1} − \frac{1}{3}) e^{(−3 t)} + (C_{3, 1} − \frac{5}{3}) e^{(3 t)} + 4 t + 2 \end{matrix} ]

or by the unknown coefficient for each eigenvector

> x := map( collect, x, {C[1,1], C[2,1], C[3,1]} );

x := [ \begin{matrix}e^{(−3 t)} C_{1, 1} − e^{(3 t)} C_{2, 1} + e^{(3 t)} C_{3, 1} + 5 t + 1 + \frac{1}{3} e^{(−3 t)} − \frac{4}{3} e^{(3 t)} \\ e^{(−3 t)} C_{1, 1} + e^{(3 t)} C_{2, 1} + 2 t + \frac{1}{3} e^{(−3 t)} − \frac{1}{3} e^{(3 t)} \\ −e^{(−3 t)} C_{1, 1} + e^{(3 t)} C_{3, 1} + 4 t + 2 − \frac{1}{3} e^{(−3 t)} − \frac{5}{3} e^{(3 t)} \end{matrix} ]

5.4 Point to Ponder

The solution found above is not the same as the one found in Nagle and Saff. Find the relationship between the constants C[1,1], C[2,1], and C[3,1] and the constants c_{1}, c_{2}, and c_{3} in Nagle and Saff.

Note that X(0) ≠ I. A fundamental matrix with X(0) = I is

> XX := evalm( X(t) &* inverse( X(0) ) );

XX := [ \begin{matrix}\frac{1}{3} e^{(−3 t)} + \frac{2}{3} e^{(3 t)} & \frac{1}{3} e^{(−3 t)} − \frac{1}{3} e^{(3 t)} & %1 \\ \frac{1}{3} e^{(−3 t)} − \frac{1}{3} e^{(3 t)} & \frac{1}{3} e^{(−3 t)} + \frac{2}{3} e^{(3 t)} & %1 \\ %1 & %1 & \frac{1}{3} e^{(−3 t)} + \frac{2}{3} e^{(3 t)} \end{matrix} ]

%1 := −\frac{1}{3} e^{(−3 t)} + \frac{1}{3} e^{(3 t)}

The fact that this matrix is a fundamental matrix can be checked as before:

   > CHECK2 := evalm( map( diff, XX, t ) - A &* XX );

CHECK2 := [ \begin{matrix}0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{matrix} ]

Further, checking that this matrix reduces to the identity when t=0 we use, for the first time, map2 (see ?map):

   > CHECK3 := evalm( map2( subs, t=0, XX ) );

CHECK3 := [ \begin{matrix}e^{0} & 0 & 0 \\ 0 & e^{0} & 0 \\ 0 & 0 & e^{0} \end{matrix} ]

5.5 Point to Ponder

The matrix XX was formed by matrix multiplication, on the right, by the inverse of the matrix at the initial time (t=0) . Is the order of these terms important? (Explain.)

6 Fundamental Solution by Matrix Exponentials

An alternate means of obtaining a fundamental solution for a system is to use the matrix exponential: X(t) = e^{(A t)}. The matrix exponential, introduced in Section 9.7 of Nagle and Saff, is found using the exponential command from the linalg package:

   > X2 := exponential( A*t );

X2 := [ \begin{matrix}\frac{1}{3} e^{(−3 t)} + \frac{2}{3} e^{(3 t)} & \frac{1}{3} e^{(−3 t)} − \frac{1}{3} e^{(3 t)} & %1 \\ \frac{1}{3} e^{(−3 t)} − \frac{1}{3} e^{(3 t)} & \frac{1}{3} e^{(−3 t)} + \frac{2}{3} e^{(3 t)} & %1 \\ %1 & %1 & \frac{1}{3} e^{(−3 t)} + \frac{2}{3} e^{(3 t)} \end{matrix} ]

%1 := −\frac{1}{3} e^{(−3 t)} + \frac{1}{3} e^{(3 t)}

6.1 Point to Ponder

The matrix exponential is the same as the fundamental matrix, XX, found at the end of the previous section.

• Is this a coincidence, or is this a general property of fundamental matrices and matrix exponentials? (Explain.)

The remainder of the compuations (homogeneous and particular solutions) are identical to the ones presented earlier, so are omitted. Note that the exact form of the general solution will be different from, but equivalent to, the one found above.

The only work that needs to be done to make this result easier to use is to convert it into a function. The unapply command must be applied to each element of the fundamental matrix, so map is needed:

   > XX2 := map( unapply, evalm(X2), t ):
   > XX2(t);

[ \begin{matrix}\frac{1}{3} e^{(−3 t)} + \frac{2}{3} e^{(3 t)} & \frac{1}{3} e^{(−3 t)} − \frac{1}{3} e^{(3 t)} & %1 \\ \frac{1}{3} e^{(−3 t)} − \frac{1}{3} e^{(3 t)} & \frac{1}{3} e^{(−3 t)} + \frac{2}{3} e^{(3 t)} & %1 \\ %1 & %1 & \frac{1}{3} e^{(−3 t)} + \frac{2}{3} e^{(3 t)} \end{matrix} ]

%1 := −\frac{1}{3} e^{(−3 t)} + \frac{1}{3} e^{(3 t)}

7 Fundamental Solution by Laplace Transforms

The Laplace transform is an altogether different approach to the problem of finding a fundamental matrix. If X’ = A X, then s XX(s)−X(0) = A XX(s) and so XX(s) = −(A − sI )^{(−1)} X(0), where the Laplace transform of X(t) is denoted by XX(s). The fundamental matrix is obtained by applying the inverse Laplace transform componentwise to the negative of the inverse of the characteristic matrix, - (A − s I)^{(−1)}.

Here is how this could be implemented:

Recall that the characteristic matrix is

> evalm( CHARm );

[ \begin{matrix}1 − λ & −2 & 2 \\ −2 & 1 − λ & 2 \\ 2 & 2 & 1 − λ \end{matrix} ]

The inverse matrix is, therefore:

> invCHARm := inverse( CHARm );

invCHARm := [ \begin{matrix}−\frac{λ + 1}{λ^{2} − 9} & \frac{2}{λ^{2} − 9} & −\frac{2}{λ^{2} − 9} \\ \frac{2}{λ^{2} − 9} & −\frac{λ + 1}{λ^{2} − 9} & −\frac{2}{λ^{2} − 9} \\ −\frac{2}{λ^{2} − 9} & −\frac{2}{λ^{2} − 9} & −\frac{λ + 1}{λ^{2} − 9} \end{matrix} ]

The inverse Laplace transform can be determined by inspection to the entries of this matrix once the negative of the inverse of the characteristic matrix is written with the partial fraction decomposition

> map( convert, evalm(-invCHARm), parfrac, lambda);

[ \begin{matrix}\frac{2}{3} \frac{1}{λ − 3} + \frac{1}{3} \frac{1}{λ + 3} & −\frac{1}{3} \frac{1}{λ − 3} + \frac{1}{3} \frac{1}{λ + 3} & %1 \\ −\frac{1}{3} \frac{1}{λ − 3} + \frac{1}{3} \frac{1}{λ + 3} & \frac{2}{3} \frac{1}{λ − 3} + \frac{1}{3} \frac{1}{λ + 3} & %1 \\ %1 & %1 & \frac{2}{3} \frac{1}{λ − 3} + \frac{1}{3} \frac{1}{λ + 3} \end{matrix} ]

%1 := \frac{1}{3} \frac{1}{λ − 3} − \frac{1}{3} \frac{1}{λ + 3}

or, letting Maple do the work for us,

> with( inttrans ):
Warning, new definition for hilbert

> X3 := map( invlaplace, evalm(-invCHARm), lambda, t );

X3 := [ \begin{matrix}\frac{1}{3} e^{(−3 t)} + \frac{2}{3} e^{(3 t)} & \frac{1}{3} e^{(−3 t)} − \frac{1}{3} e^{(3 t)} & %1 \\ \frac{1}{3} e^{(−3 t)} − \frac{1}{3} e^{(3 t)} & \frac{1}{3} e^{(−3 t)} + \frac{2}{3} e^{(3 t)} & %1 \\ %1 & %1 & \frac{1}{3} e^{(−3 t)} + \frac{2}{3} e^{(3 t)} \end{matrix} ]

%1 := −\frac{1}{3} e^{(−3 t)} + \frac{1}{3} e^{(3 t)}

Once again, this result agrees with the previous results. The remainder of the steps in the determination of the general solution are the same – don’t forget to convert X3 to a function before attempting to do too many computations with this matrix.

8 Concluding Remarks

It has already been noted that XX = X2 = X3 . In general, any of the three methods will lead to equivalent general solutions for the nonhomogeneous problem.

If A has an eigenvalue which does not have a full set of eigenvectors, then the explicit computation of a fundamental matrix requires special handling. Likewise, when A has complex eigenvalues (and eigenvectors), real-valued solutions are obtained only after careful handling of the real and imaginary parts of the eigenvalues and eigenvectors. In either case, the matrix exponential, as computed by exponential, and the Laplace transform avoid this complication.

9 Copyright

systems.mws – a Maple worksheet for use with systems of ODEs and IVPs (for use with section 9.4–9.7, pp. 533–564, of Nagle and Saff) Prepared by: Douglas B. Meade (7 March 1996) as part of the Maple Supplement for Nagle/Saff (Addison-Wesley) Updated to Maple V, Release 4, by Douglas B. Meade (27 December, 1996) For more information, please contact the author at one of following addresses: Department of Mathematics, Univ. of South Carolina, Columbia, SC 29208 (803) 777-6183 meade@math.sc.edu (copyright, Douglas B. Meade, 1996)