ORDINARY DIFFERENTIAL EQUATIONS POWERTOOL

Unit 11 -- First-Order Linear Systems

Outline of Unit 11

Initialization

> restart;

> with( DEtools ):

> with( plots ):

> with( linalg ):

> with( student ):

Warning, the name changecoords has been redefined

Warning, the name adjoint has been redefined

Warning, the protected names norm and trace have been redefined and unprotected

> pvac := true:

>

A nullcline for a first-order system of differential equations is a curve along which at least one of the variables does not change, e.g., where or .The key is to assume the first derivative of each unknown function is zero and to find all curves on which the resulting equation is satisfied. While it can be difficult to give a general outline for the determination of the nullclines for a general nonlinear system, for a linear system, the nullclines are straight lines.

Consider the general first-order linear system with constant coefficients,

> ode1 := diff( x(t), t ) = a*x(t) + b*y(t) + f:

> ode2 := diff( y(t), t ) = c*x(t) + d*y(t) + g:

> sys := { ode1, ode2 };

>

The nullclines of the system are

> nullcline := eval( sys, { diff(x(t),t)=0, diff(y(t),t)=0, x(t)=x, y(t)=y } );

>

The change from to and from to is made to simplify the graphical display of the nullclines.

For a concrete example, let

> value_of_constants := {a=1,b=-2,c=3,d=1,f=2,g=1};

>

then the system is

> sys1 := eval( sys, value_of_constants );

>

and the nullclines are

> nullcline1 := eval( sys1, { diff(x(t),t)=0, diff(y(t),t)=0, x(t)=x, y(t)=y } );

>

or, expressing each line in slope-intercept form.

> map( isolate, nullcline1, y );

>

A graphical view of the nullclines is obtained by creating an implicitplot of each nullcline. (The implicitplot command is needed to handle cases where the nullcline is a vertical line.)

> WINDOW := x=-5..5, y=-5..5:

> plot_nullcline := display( map( implicitplot, nullcline1, WINDOW, color=GREEN ) ):

> plot_nullcline;

>

Equilibrium solutions (if they exist) occur at the intersection of nullclines corresponding to each of the differential equations in the system. For a general nonlinear system, special care must be taken to ensure that an intersection point of nullclines is actually an equilibrium solution. Fortunately, for a two-dimensional linear system, any intersection of the nullclines is automatically an equilibrium solution.

In terms of the example introduced above,

> sys1;

>

with nullclines

> nullcline1;

>

The intersection of the nullclines can be approximated from the plot of the nullclines. The exact solution is found to be

> equil_soln := solve( nullcline1, {x,y} );

>

or, expressed as the coordinates of a point,

> equil_point := eval( [x,y], equil_soln );

>

Graphically, the equilibrium solution is the intersection of the two nullclines.

> plot_equil := pointplot( equil_point, symbol=CIRCLE, symbolsize=18, color=BLUE ):

> display([plot_nullcline,plot_equil]);

>

The direction field can be a source of much information about a system of ODEs without knowing an explicit formula for the solution. For starters, the direction field for a system shows the direction in which the solution moves at any point in space.

A phase portrait for an autonomous system of ODEs displays solutions to the system in "phase space". That is, as curves in the space of the unknown functions ( and ) and not the independent variable ( ).

To create a phase portrait in Maple, the DEplot command is recommended. (The phaseportrait command, also from the DEtools package, is essentially the same, but there is no reason to learn and remember a new command when simple modifications to DEplot suffice.)

One solution curve is generated for each initial condition. For a two-dimensional system in x(t) and y(t), each initial condition can be specified either as triples of numbers or as pairs of equations . For example, initial conditions at the points with integer coordinates along the axes can be specified as

> ICy := [0,0,i] $ i=-5..5;

> ICx := [x(0)=i,y(0)=0] $ i=-5..5;

> IC := ICx, ICy:

>

As in the previous uses of DEplot, an interval of values for the independent variable is required. Unlike the use of DEplot for the creation of direction fields, this argument does have a meaning for phase portraits. The time interval provides limits for the numerical methods used to obtain approximate solutions for each initial condition.

> DOMAIN;

>

The phase portrait for this system and these initial conditions is created with

> DEplot( sys1, [x(t),y(t)], DOMAIN, [ IC ],

> arrows=NONE, scene=[ x, y ], linecolor=BLUE );

>

To restrict the solutions to a pre-determined "window", include the viewing window for the unknown functions:

> WINDOW;

>

> plot_phase_portrait := DEplot( sys1, [x(t),y(t)], DOMAIN, [ IC ], WINDOW,

> arrows=NONE, scene=[ x, y ], linecolor=BLUE ):

> plot_phase_portrait;

>

A simple change to the arrows= option includes the direction field in the phase portrait.

> DEplot( sys1, [x(t),y(t)], DOMAIN, [ IC ], WINDOW,

> arrows=SMALL, scene=[ x, y ], linecolor=BLUE );

>

This plot can be superimposed on the direction field, nullclines, and equilibrium solution with

> display( [plot_direction_field, plot_nullcline,

> plot_equil, plot_phase_portrait] );

>

Even though the the independent variable ( ) is not explicitly displayed in a phase portrait, it is implicitly present in that the solution curves are traveled at different speeds. To show the coordination between solutions from different initial conditions, specify the linecolor= option with a function or expression that depends on the independent variable. For example,

> DEplot( sys1, [x(t),y(t)], DOMAIN, [ IC ],

> arrows=NONE, scene=[ x, y ], linecolor=t );

Note that it is not necessary to be too concerned about the specific range of values of the linecolor= option; DEplot automatically normalizes these values to the interval [0,1].

>

> WINDOW2 := x=-20..20, y=-20..20:

> phase_portrait := T -> DEplot( sys1, [x(t),y(t)], t=0..T, [ IC ],

> WINDOW2, arrows=SLIM, scene=[ x, y ],

> linecolor=BLUE, stepsize=0.1 ):

> display( [ phase_portrait(0.1), seq( phase_portrait(i/4), i=1..10) ], insequence=true );

>

Note that to provide a consistent resolution for all solution curves in the animation, the stepsize= option has been used.

>

A solution curve is a plot of one component of the solution as a function of the independent variable and is produced using the DEplot command. The only changes to the arguments are removing the specification of the display window size and the arrows= argument and changing the scene= argument to specify the component to be plotted.

For example, when the initial condition is

> IC := [x(0)=0,y(0)=1];

>

graphs of the - and -components of the particular solution satisfying this initial condition are created and displayed with

> plot_soln_x := DEplot( sys1, [x(t),y(t)], DOMAIN, [ IC ],

> scene=[ t, x ], linecolor=BLUE ):

> plot_soln_y := DEplot( sys1, [x(t),y(t)], DOMAIN, [ IC ],

> scene=[ t, y ], linecolor=GREEN ):

> display( [ plot_soln_x, plot_soln_y ], labels=[`t`,``] );

If multiple initial conditions are specified, one solution curve is produced for each initial condition.

> IC2 := [x(0)=0,y(0)=1], [x(0)=0,y(0)=-1];

>

> plot_soln_x2 := DEplot( sys1, [x(t),y(t)], DOMAIN, [ IC2 ],

> scene=[ t, x ], linecolor=[BLUE,GREEN] ):

> plot_soln_y2 := DEplot( sys1, [x(t),y(t)], DOMAIN, [ IC2 ],

> scene=[ t, y ], linecolor=[BLUE,GREEN] ):

>

These plots can be displayed as before

> display( [ plot_soln_x2, plot_soln_y2 ] , labels=[`t`,``] );

>

To eliminate some of the clutter from the multiple plots, it is sometimes advantageous to display each component in side-by-side plots.

> display( array([plot_soln_x2,plot_soln_y2]) );

Error, (in plot/plotarray/niceticks) cannot evaluate boolean: -trunc(1/5*(undefined-5*10^floor((undefined+undefined*I)/ln(10))*ceil(undefined/(10^floor((undefined+undefined*I)/ln(10)))))/(10^floor((undefined+undefined*I)/ln(10)))) < -3

> plot_soln_x2; plot_soln_y2;

>

Note how the colors are used to indicate pairs of solutions corresponding to the same initial condition.

>

For a first-order linear system of ODEs, Maple's dsolve command should be able to find the general solution to the system and the particular solution for any initial condition.

For example, for the system of ODEs

> sys1;

>

the general solution is

> gen_soln := dsolve( sys1, {x(t),y(t)} );

>

The solution satisfying an initial condition, say

> IC;

>

can be found by substituting the initial condition into the general solution

> q1 := eval( eval( gen_soln, t=0 ), IC );

>

and solving for the constants of integration

> q2 := solve( q1, {_C1,_C2} );

>

The corresponding particular solution to this system of differential equations that satisfies this initial condition is

> part_soln := eval( gen_soln, q2 );

>

Alternatively, the solution to the initial value problem

> IVP := sys1 union convert(IC,set);

>

could be found with the single command

> ivp_soln := dsolve( IVP, {x(t),y(t)} );

>

While the odetest command is unable to check solutions for a system of ODEs, it is not difficult to verify that the two solutions presented above are identical

> q3 := eval( x(t), part_soln ) = eval( x(t), ivp_soln );

> evalb(q3);

> q4 := eval( y(t), part_soln ) = eval( y(t), ivp_soln );

> evalb(q4);

>

and satisfy both the system of differential equations

> q5 := simplify( eval( convert(sys1,list), ivp_soln ) );

> map( evalb, q5 );

>

and the initial condition

> q6 := eval( ivp_soln, t=0 );

> map( evalb, eval( IC, q6 ) );

>

The solution of a first-order linear system of ODEs is possible because the problem can be reduced to finding the eigenvalues and eigenvectors of an appropriate matrix and then knowing how to use this information to construct the general or particular solution. To put the problem in the context of linear algebra and eigenvalues, the system of ODEs should be written in vector form.

The general form for a linear system of ODEs is

> Diff( X, t ) = A*X + b;

>

where is the vector of unknown functions, is the coefficient matrix and is the non-homogeneous ("forcing") vector.

>

Example 1: Real and Distinct Eigenvalues

Consider the system of ODEs

> sys2 := [ diff( x(t), t ) = x(t) + y(t), diff( y(t), t ) = x(t) + y(t) ];

>

Let the vector of unknown functions be

> X := vector( 2, [ x(t), y(t) ] );

>

The coefficient matrix, , and forcing vector, , can be extracted via the genmatrix command from the linalg package:

> A := genmatrix( map(rhs,sys2), [x(t),y(t)], `-b` );

> b := evalm( -`-b` );

>

Thus, the system can be expressed in vector form as

> map( diff, evalm(X), t ) = evalm(A) * evalm(X) + evalm(b);

>

The general solution to this system of ODEs is

> gen_soln := dsolve( sys2, {x(t),y(t)} );

>

which, when written in vector form

> X_soln := eval( evalm(X), gen_soln );

>

and separated into terms involving at most one of the two constants of integration, can be written as

> X1 := map( coeff, X_soln, _C1 ):

> X2 := map( coeff, X_soln, _C2 ):

> Xp := eval( evalm(X_soln), [_C1=0,_C2=0] ):

> evalm(X) = _C1 * evalm(X1) + _C2 * evalm(X2) + evalm(Xp);

>

Note that, because this is a homogeneous system, the zero vector is a particular solution and the general solution is written as a linear combination of two vectors. The significance of this form of the solution appears when the terms in the solution are examined side-by-side with the eigenvalue decomposition of the coefficient matrix.

At first glance, the output from Maple's eigenvectors command is hopelessly complicated.

> e_decomp := [ eigenvectors( A ) ];

>

However, it is not difficult to extract the eigenvalues ( ) and corresponding basis of eigenvectors ( ) to obtain

> for i from 1 to nops(e_decomp) do

> e_val[i] := op(1,e_decomp[i]);

> e_vec[i] := op(3,e_decomp[i]);

> e_sol[i] := exp( e_val[i]*t) * op(e_vec[i]);

> print( lambda[i]=e_val[i], E[i]=evalm(e_vec[i]), XX[i] = evalm(e_sol[i]) );

> od:

>

Note that the vectors and are the basis vectors for the linear combination that is the general solution of this homogeneous system.

>

Example 2: Complex Eigenvalues

Returning to the system introduced at the beginning of this unit,

> sys1;

>

recall that this system is non-homogeneous. In fact, the coefficient matrix and forcing term are

> A := genmatrix( map(rhs,sys1), [x(t),y(t)], `-b` );

> b := evalm( -`-b` );

>

As seen previously, the general solution of this system is

> gen_soln := dsolve( sys1, {x(t),y(t)} );

>

which, when written in vector form and factored, appears as

> X_soln := eval( evalm(X), gen_soln );

>

Since this system is non-homogeneous, the particular solution is non-trivial. In this case, one particular solution is

> Xp := eval( evalm(X_soln), [_C1=0,_C2=0] );

>

and a basis for the solution of the corresponding homogeneous system is

> X1 := map( coeff, X_soln, _C1 );

> X2 := map( coeff, X_soln, _C2 );

>

This leads to the following vector form for the general solution of this system

> evalm(X) = _C1 * evalm(X1) + _C2 * evalm(X2) + evalm(Xp);

>

To understand the relationship between the general solution and the eigenvalue decomposition of the coefficient matrix, consider

> e_decomp := [ eigenvectors( A ) ];

>

The eigenvalues, and hence the eigenvalues, appear in complex conjugate pairs.

> for i from 1 to nops(e_decomp) do

> e_val[i] := op(1,e_decomp[i]);

> e_vec[i] := op(3,e_decomp[i]);

> e_sol[i] := exp( e_val[i]*t ) * evalm(op(e_vec[i]));

> print( lambda[i]=e_val[i], E[i]=evalm(e_vec[i]), XX[i]=evalm(e_sol[i]) );

> od:

>

The functions and are solutions to the system of ODEs, but they are complex-valued. Real-valued solutions, such as the ones returned by dsolve , would be more useful.

By Euler's formula, if and are real numbers, then

To apply this result in our case, it is necessary to tell Maple that the independent variable, , is real-valued. Then, the real and imaginary parts of one of the complex-valued solutions are

> assume( t, real );

> XX[1] := map( Re@evalc, evalm(e_sol[1]) );

> XX[2] := map( Im@evalc, evalm(e_sol[1]) );

> unassign('t');

>

These functions form a real-valued basis for the general solution of the homogeneous system. (These solutions may appear to differ from the ones reported by dsolve ; closer inspection reveals that these vectors are, at worst, parallel to the ones found by dsolve and so have the same linear span.)

>

Example 3: Repeated Eigenvalue

For a final example, consider the homogeneous system

> sys3 := [ diff( x(t), t ) = x(t) - 2*y(t), diff( y(t), t ) = y(t) ];

>

which has coefficient matrix and forcing vector

> A := genmatrix( map(rhs,sys3), [x(t),y(t)], `-b` );

> b := evalm( -`-b` );

>

The general solution to the system, as reported by dsolve , is

> infolevel[dsolve] := 3:

> gen_soln := dsolve( sys3, {x(t),y(t)} );

> infolevel[dsolve] := 0:

-> Solving ordering for ODE system:   [y(t), x(t)]

Methods for first order ODEs:

Trying to isolate the derivative dY/dX...

Successful isolation of dY/dX

-> Trying classification methods

trying a quadrature

trying 1st order linear

1st order linear successful

Methods for first order ODEs:

Trying to isolate the derivative dY/dX...

Successful isolation of dY/dX

-> Trying classification methods

trying a quadrature

trying 1st order linear

1st order linear successful

>

A particular solution to the system is the trivial solution

> X_soln := eval( evalm(X), gen_soln ):

> Xp := eval( evalm(X_soln), [_C1=0,_C2=0] );

>

Note that any values for _C1 and _C2 will yield a particular solution; using zero for all constants of integration is the easiest and most common choice.

A basis for the homogeneous solution is formed by the pair of functions

> X1 := map( coeff, X_soln, _C1 );

> X2 := map( coeff, X_soln, _C2 );

>

When these pieces are assembled, the general solution of the system can be written in the form

> evalm(X) = _C1 * evalm(X1) + _C2 * evalm(X2) + evalm(Xp);

>

Notice that both terms in the homogeneous solution involve the same exponential, , and, after factoring the exponential, one of the homogeneous terms has a coefficient that is not constant. These features will have to be explained during the eigenvalue decomposition.

> e_decomp := [ eigenvectors( A ) ];

>

> for i from 1 to nops(e_decomp) do

> e_val[i] := op(1,e_decomp[i]);

> e_vec[i] := op(3,e_decomp[i]);

> e_sol[i] := exp( e_val[i]*t ) * op(e_vec[i]);

> print( lambda[i]=e_val[i], E[i]=evalm(e_vec[i]), XX[i] = evalm(e_sol[i]) );

> od:

>

Because is an eigenvalue with algebraic multiplicity 2 and geometric multiplicity 1, the eigenvalue decomposition yields only one solution to the homogeneous equation. This solution

> XX[1] = evalm(e_sol[1]);

>

is referred to as a straight-line solution. The second solution for the basis of the homogeneous solution will be found in the form

> sol2_form := exp(e_val[1]*t) * (t*op(e_vec[1])+V2);

>

where

> V2 := vector( 2, [x2,y2] );

>

That is, assume

> r1 := equate( X, evalm(sol2_form) );

>

Now, substitution of the proposed solution into the (homogeneous) system of ODEs yields

> sol2_requires := eval( sys3, r1 );

>

Note that the second condition is trivially satisfied for all values of y2. However, the first condition is satisfied only when

> sol2_satisfied := solve( identity(sol2_requires[1],t), {x2,y2} );

>

This leads to the one-parameter family of solutions

> sol2_family := eval( evalm(sol2_form), sol2_satisfied );

>

Notice that the component of the solution that depends on the remaining parameter

> map( coeff, sol2_family, x2 );

>

is a multiple of the first solution to the homogeneous system

> evalm( e_sol[1] );

>

and so is not needed again in the basis. Therefore, the second basis solution for the homogeneous solution is chosen to be

> e_sol[2] := eval( evalm(sol2_family), x2=0 );

>

To summarize, the basis of homogeneous solutions found by this method is

> { evalm(e_sol[1]), evalm(e_sol[2]) };

>

and the basis of solutions found by inspection of the dsolve solution is

> { evalm(X1), evalm(X2) };

>

The spans of these bases are easily seen to be identical.

>

As a final note, observe that this system is decoupled. The ODE for is independent of x(t)

> ode_x,ode_y := selectremove( has, sys3, diff(x(t),t) );

>

and can be solved explicitly without knowledge of

> sol_y := dsolve( ode_y, y(t), [linear] );

>

Now, since is known, this result can be substituted into the ODE for

> ode_x2 := eval( ode_x, sol_y );

>

and solved for

> sol_x := dsolve( ode_x2, x(t) );

>

The result is

> eval( evalm(X), {sol_x,sol_y} );

>

which is equivalent to any of the solutions obtained above.

> evalm(X_soln);

>

The graphical solutions produced by DEplot are obtained using a numerical approximation to the solution. The numerical method used to compute the approximations can be specified via the method= argument (the default is method=classic[rk4] ). The dsolve command, with type=numeric , can be used to obtain direct access to a numerical solution to an initial value problem.

For the initial value problem

> IVP;

>

a procedure for the numerical approximation to the solution via Euler's method is obtained with

> soln_euler := dsolve( IVP, [x(t),y(t)], type=numeric, method=classical );

>

The approximate solution at is

> soln_euler(1);

>

The odeplot command can be used to create a phase-portrait for the solution

> odeplot( soln_euler, [x(t),y(t)], 0..4 );

>

or a plot of the individual components of the solution

> odeplot( soln_euler, [[t,x(t)],[t,y(t)]], 0..4,

> labels=[`t`,``], legend=[`x(t)`,`y(t)`] );

>