ORDINARY DIFFERENTIAL EQUATIONS POWERTOOLLesson 11 -- First-Order Linear SystemsProf. Douglas B. MeadeIndustrial Mathematics InstituteDepartment of MathematicsUniversity of South CarolinaColumbia, SC 29208 URL: http://www.math.sc.edu/~meade/E-mail: meade@math.sc.edu Copyright \251 2001 by Douglas B. MeadeAll rights reserved-------------------------------------------------------------------

11.A Equilibrium Analysis 11.A-1 Nullclines 11.A-2 Equilibrium Solutions11.B Graphical Analysis 11.B-1 Direction Fields 11.B-2 Phase Portraits 11.B-3 Solution Curves11.C Analytic Solutions 11.C-1 One-Step Solutions using dsolve 11.C-2 Eigenvalue Analysis Example 1 Real and Distinct Eigenvalues Example 2 Complex Eigenvalues Example 3 Repeated Eigenvalues11.D Numerical Solutions

restart;with( DEtools ):with( plots ):with( LinearAlgebra ):with( student ):

In the NiMlI3h5Rw==-plane, an isocline for the differential equation NiMvKiYlI2R5RyIiIiUjZHhHISIiLSUiZkc2JCUieEclInlH is a curve along which NiMqJiUjZHlHIiIiJSNkeEchIiI= has a constant value. Thus, an isocline is a level curve of the function NiMtJSJmRzYkJSJ4RyUieUc=.If the contstant value of NiMtJSJmRzYkJSJ4RyUieUc= on the isocline is zero, then the isocline is called a nullcline.A nullcline for a two-dimensional first-order system of differential equations is a curve along which at least one of the dependent variables (NiMlInhH or NiMlInlH) does not change, i.e., where NiMvKiYlI2R4RyIiIiUjZHRHISIiIiIh or NiMvKiYlI2R5RyIiIiUjZHRHISIiIiIh.The nullclines for the systemNiMvKiYlI2R4RyIiIiUjZHRHISIiLSUiRkc2JCUieEclInlH NiMvKiYlI2R5RyIiIiUjZHRHISIiLSUiR0c2JCUieEclInlH are found by graphing the curves NiMtJSJ5RzYjJSJ4Rw== defined implicitly by the equations NiMvLSUiRkc2JCUieEclInlHIiIh and NiMvLSUiR0c2JCUieEclInlHIiIh. While it can be difficult to give a general description of the nullclines for a general nonlinear system, the nullclines for a linear system are always straight lines that intersect at an equilibrium solution.Consider the general first-order linear system with constant coefficients, namely,NiMvKiYlI2R4RyIiIiUjZHRHISIiLCgqJiUiYUdGJi0lInhHNiMlInRHRiZGJiomJSJiR0YmLSUieUdGLkYmRiYlImZHRiY= NiMvKiYlI2R5RyIiIiUjZHRHISIiLCgqJiUiY0dGJi0lInhHNiMlInRHRiZGJiomJSJkR0YmLSUieUdGLkYmRiYlImdHRiY= which are rendered in Maple viaode1 := 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 the straight lines whose equations arenullcline := eval( sys, { diff(x(t),t)=0, diff(y(t),t)=0, x(t)=x, y(t)=y } );The change from NiMtJSJ4RzYjJSJ0Rw== to NiMlInhH and from NiMtJSJ5RzYjJSJ0Rw== to NiMlInlH is made to simplify the graphical display of the nullclines.For a concrete example, select the following values of the coefficients.value_of_constants := {a=1,b=-2,c=3,d=1,f=2,g=1};Then the system issys1 := eval( sys, value_of_constants );and the nullclines arenullcline1 := 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 );These nullclines are graphed in Figure 11.1, created from 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 ) ):display(plot_nullcline, title="Figure 11.1");

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, namely,sys1;with nullclinesnullcline1;the intersection of the nullclines can be approximated from Figure 11.1, the plot of the nullclines. The exact solution is found to beequil_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, as shown in Figure 11.2.plot_equil := pointplot( equil_point, symbol=CIRCLE, symbolsize=18, color=BLUE ):display([plot_nullcline,plot_equil], title="Figure 11.2");

The direction field (see Lesson 1) can be a source of much information about a system of ODEs, even when an explicit formula for the solution is not known. For starters, the direction field for a system shows the direction in which the solution moves at any point in space.Again consider the system sys1;first seen in Lesson 11, Section A . The aim here is to construct and study the direction field for this system.Even though this is an autonomous system (recall from Lesson 1 that the independent variable does not explicitly appear), Maple's DEplot command requires an interval for the independent variableNiMtJSFHNiMlInRH. The specific interval chosen is irrelevant, provided the endpoints are numeric. Hence, pick the domainDOMAIN := t = 0 .. 1:The direction field for this system appears in Figure 11.3.plot_direction_field := DEplot( sys1, [x(t),y(t)], DOMAIN, WINDOW, arrows = MEDIUM ):display(plot_direction_field, title="Figure 11.3");From this picture it is evident that solutions spiral outward from a point in the second quadrant. To see that the center of these spirals is the equilibrium solution, superimpose the plot of the equilibrium solution and the nullclines, as done in Figure 11.4.display( [plot_direction_field, plot_equil, plot_nullcline], title="Figure 11.4" );The nullclines show exactly where the individual components pass through critical points and, hence, where they transition between increasing and decreasing functions (of NiMlInRH). In practice, if a computer is not available to create the direction field, the constancy of the sign of the derivative of each component in each region created by the nullclines provides insight into the qualitative behaviour of solutions to a linear, autonomous system of ODEs.

A phase portrait for an autonomous system of ODEs displays solutions to the system in "phase space" where solutions NiMtJSJ4RzYjJSJ0Rw== and NiMtJSJ5RzYjJSJ0Rw== are treated as parametric equations for the curve NiMtJSJ5RzYjJSJ4Rw==. (See Lesson 1.) 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 NiMtJSJ4RzYjJSJ0Rw== and NiMtJSJ5RzYjJSJ0Rw==, each initial condition can be specified either as triples of numbers NiM3JSYlInRHNiMiIiEtJSJ4RzYjRiQtJSJ5R0Yq or as pairs of equations NiM3JC8tJSJ4RzYjJiUidEc2IyIiISZGJkYqLy0lInlHRicmRi9GKg==. For example, initial conditions at the points with integer coordinates along the axes can be specified asICy := [0,0,i] $ i=-5..5;for the NiMlInlH-axis,ICx := [x(0)=i,y(0)=0] $ i=-5..5;for the NiMlInhH-axis, andIC := ICx, ICy:as a composit list for both axes.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.Using the domainDOMAIN;for the independent variable, the phase portrait for this system and these initial conditions appears in Figure 11.5.DEplot( sys1, [x(t),y(t)], DOMAIN, [ IC ], arrows=NONE, scene=[ x, y ], linecolor=BLUE, title="Figure 11.5" );To restrict the solutions to a pre-determined "window", include a viewing window such asWINDOW;This gives Figure 11.6.plot_phase_portrait := DEplot( sys1, [x(t),y(t)], DOMAIN, [ IC ], WINDOW, arrows=NONE, scene=[ x, y ], linecolor=BLUE ):display(plot_phase_portrait, title="Figure 11.6");A simple change to the arrows= option includes the direction field in the phase portrait, as seen in Figure 11.7.DEplot( sys1, [x(t),y(t)], DOMAIN, [ IC ], WINDOW, arrows=SMALL, scene=[ x, y ], linecolor=BLUE, title="Figure 11.7" );In Figure 11.8, this plot is superimposed on the direction field, nullclines, and equilibrium solution.display( [plot_direction_field, plot_nullcline, plot_equil, plot_phase_portrait], title="Figure 11.8" );Even though the the independent variableNiMtJSFHNiMlInRH 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, this is done in Figure 11.9.DEplot( sys1, [x(t),y(t)], DOMAIN, [ IC ], arrows=NONE, scene=[ x, y ], linecolor=t, title="Figure 11.9" ); 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].Figure 11.9 animates the tracing of the trajectories in the phase portrait.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, title="Figure 11.9" );Note that to provide a consistent resolution for all solution curves in the animation, the stepsize= option has been used.

For the two-dimensional linear system, 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 isIC := [x(0)=0,y(0)=1];Figure 11.10, containing graphs of the NiMlInhH- and NiMlInlH-components of the particular solution satisfying this initial condition, is created and displayed withplot_soln_x := DEplot( sys1, [x(t),y(t)], DOMAIN, [ IC ], scene=[ t, x ], linecolor=black ):plot_soln_y := DEplot( sys1, [x(t),y(t)], DOMAIN, [ IC ], scene=[ t, y ], linecolor=red ):display( [ plot_soln_x, plot_soln_y ], labels=[`t`,``], title="Figure 11.10" );If multiple initial conditions are specified as inIC2 := [x(0)=0,y(0)=1], [x(0)=0,y(0)=-1];one solution curve is produced for each initial condition. The individual plots are created withplot_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] ):and can be displayed together, as in Figure 11.11.display( [ plot_soln_x2, plot_soln_y2 ] , labels=[`t`,``], title="Figure 11.11" );To eliminate some of the clutter from the multiple plots, it is sometimes advantageous to display each pair of solutions in side-by-side plots, as seen in the following figure.display( array([plot_soln_x2,plot_soln_y2]), scaling=constrained );Note how the colors are used to indicate pairs of solutions corresponding to the same initial condition. The blue curve on the left is the graph of NiMtJSJ4RzYjJSJ0Rw== whereas the blue curve on the right is the graph of NiMtJSJ5RzYjJSJ0Rw== both of which correspond to the first initial condition. The pair of green curves corresponds to the second 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 ODEssys1;the general solution isgen_soln := dsolve( sys1, {x(t),y(t)} );more easily grasped if written asxt = eval(x(t), gen_soln); yt = eval(y(t), gen_soln);The solution satisfying an initial condition, sayIC;can be found by substituting the initial condition into the general solution, thereby producing the equationsq1 := eval( eval( gen_soln, t=0 ), IC );which can then be solved for the constants of integration, namely,q2 := solve( q1, {_C1,_C2} );The corresponding solution of this system of differential equations ispart_soln := eval( gen_soln, q2 );again more easily grasped if written asxt = eval(x(t), part_soln); yt = eval(y(t), part_soln);Alternatively, the solution to the initial value problemIVP := sys1 union convert(IC,set);could be found with the single commandivp_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. The equivalence of NiMtJSJ4RzYjJSJ0Rw== is established withq3 := eval( x(t), part_soln ) = eval( x(t), ivp_soln );evalb(q3);whereas that of NiMtJSJ5RzYjJSJ0Rw== is established withq4 := eval( y(t), part_soln ) = eval( y(t), ivp_soln );evalb(q4);That these solutions satisfy the system of differential equations is established withq5 := simplify( eval( convert(sys1,list), ivp_soln ) );map( evalb, q5 );and that they satisfy the initial condition is established withq6 := eval( ivp_soln, t=0 );map( evalb, eval( IC, q6 ) );

If the constant-coefficient, first-order linear system of ODEs is written in the matrix formX' = A X + b where X = X(t) is the vector of unknown functions, X' = NiMqJiUiZEciIiIlI2R0RyEiIg== X is the vector of derivatives of the unknown functions, NiMlIkFH is the coefficient matrix, and b is the non-homogeneous ("forcing") vector, the solution can be constructed from the eigenvalues and eigenvectors (eigenpairs) of the matrix NiMlIkFH. For example, the matrix form of the systemNiMvLSUjeCdHNiMlInRHLCYqJiIiIyIiIiUieEdGK0YrKiYiIiRGKyUieUdGK0Yr NiMvLSUjeSdHNiMlInRHLCYqJiIiJiIiIiUieEdGK0YrKiYiIihGKyUieUdGKyEiIg== would beX' = NiMtJSdNQVRSSVhHNiM3JDcjJSN4J0c3IyUjeSdH = NiMtJSdNQVRSSVhHNiM3JDckIiIjIiIkNyQiIiYsJCIiKCEiIg== NiMtJSdNQVRSSVhHNiM3JDcjJSJ4RzcjJSJ5Rw== An eigenpair for a matrix NiMlIkFH consists of a scalar NiMlJ2xhbWJkYUc= and a direction (expressed by an eigenvector V) satisfying the defining equationA V = NiMlJ2xhbWJkYUc= V Thus, the eigenvector represents a direction that is invariant under multiplication by the matrix NiMlIkFH. At most, a vector in this direction has its length changed by the scale factor NiMlJ2xhbWJkYUc=, called the eigenvalue.If this defining equation is rewritten in the formNiMtJSFHNiMsJiUiQUciIiIqJiUnbGFtYmRhR0YoJSJJR0YoISIi V = 0 where NiMlIklH is the 2 x 2 identity matrix NiMtJSdNQVRSSVhHNiM3JDckIiIiIiIhNyRGKUYo , then the eigenvalues are found as the roots of the characteristic equation, that is, the equation NiMvLSUkZGV0RzYjLCYlIkFHIiIiKiYlJ2xhbWJkYUdGKSUiSUdGKSEiIiIiIQ==.If the eigenpairs are (NiMmJSdsYW1iZGFHNiMiIiI=, V1) and (NiMmJSdsYW1iZGFHNiMiIiM=, V2), then the general solution of the system is given by the sumNiMmJSJYRzYjJSJnRw== = NiMqJiYlImNHNiMiIiJGJy0lJGV4cEc2IyomJiUnbGFtYmRhR0YmRiclInRHRidGJw== V1 + NiMqJiYlImNHNiMiIiMiIiItJSRleHBHNiMqJiYlJ2xhbWJkYUdGJkYoJSJ0R0YoRig= V2 where NiUmJSJjRzYjJSJrRy9GJiIiIiIiIw==, are arbitrary constants.The examples below detail how the solution of the linear system can be constructed from the eigenpairs of the system matrix NiMlIkFH.

Consider the system of ODEssys2 := [ diff( x(t), t ) = x(t) + y(t), diff( y(t), t ) = x(t) + y(t) ];Let the vector of unknown functions beX := < x(t), y(t) >;The coefficient matrix NiMlIkFH, and forcing vector b can be extracted via the GenerateMatrix command from the LinearAlgebra package:(Aneg,b) := GenerateMatrix( eval(sys2,[x(t)=_x,y(t)=_y]), [_x,_y] ):A := -Aneg;Thus, the system can be expressed in vector form asmap(diff, X, t ) = A.X + b;The general solution to this system of ODEs isgen_soln := dsolve( sys2, {x(t),y(t)} );which, when written in vector formX_soln := eval( X, gen_soln );and separated into terms involving at most one of the two constants of integration, can be written asX1 := map( coeff, X_soln, _C1 ):X2 := map( coeff, X_soln, _C2 ):Xp := eval( X_soln, [_C1=0,_C2=0] ):X = _C1 * X1 + _C2 * X2 + Xp;The two vectors V1 := X1; V2 := eval(X2, t=0);are the eigenvectors corresponding, respectively, to the eigenvalueslambda[1] = 0; lambda[2] = 2;The eigenvalues are the roots of the characteristic equationCharacteristicPolynomial(A,lambda) = 0;which is Maple's built-in command for obtaining the equivalent of the equation NiMvLSUkZGV0RzYjLCYlIkFHIiIiKiYlJ2xhbWJkYUdGKSUiSUdGKSEiIiIiIQ==.That each eigenpair satisfied an equation of the form A V = NiMlJ2xhbWJkYUc= V is established with the following computations:A V1 = NiMtJSdNQVRSSVhHNiM3JDckIiIiRihGJw== NiMtJSdNQVRSSVhHNiM3JDcjIiIiNyMsJEYoISIi = NiMtJSdNQVRSSVhHNiM3JDcjIiIhRic= = NiMmJSdsYW1iZGFHNiMiIiI= V1 = 0 NiMtJSdNQVRSSVhHNiM3JDcjIiIiNyMsJEYoISIi = NiMtJSdNQVRSSVhHNiM3JDcjIiIhRic= A V2 = NiMtJSdNQVRSSVhHNiM3JDckIiIiRihGJw== NiMtJSdNQVRSSVhHNiM3JDcjIiIiRic= = NiMtJSdNQVRSSVhHNiM3JDcjIiIjRic= = NiMmJSdsYW1iZGFHNiMiIiM= V2 = 2 NiMtJSdNQVRSSVhHNiM3JDcjIiIiRic= = NiMtJSdNQVRSSVhHNiM3JDcjIiIjRic= The eigenpairs of a 2 x 2 matrix can generally be found by manual computation. However, the calculations can become tedious, and are often prone to arithmetic errors. Instead, we will rely on Maple's Eigenvectors command, the output of which must be carefully deconstructed to extract the eigenvalues and eigenvectors.Application of that command here yieldsL,V := Eigenvectors(A);Observe that this assignment creates both a vector containing the eigenvalues and the corresponding matrix of eigenvectors. It is not difficult to extract the eigenvalues (NiMmJSdsYW1iZGFHNiMlImlH) and corresponding eigenvectors (NiMmJSJFRzYjJSJpRw==). These eigenpairs, as well as the eigensolutions NiMtJSRleHBHNiMqJiYlJ2xhbWJkYUc2IyUiaUciIiIlInRHRis= NiMmJSJFRzYjJSJpRw==, are obtained and printed viafor i from 1 to Dimension(L) do e_val[i] := L[i]; e_vec[i] := Column(V,i); e_sol[i] := exp( e_val[i]*t ) * e_vec[i]; print( lambda[i]=e_val[i], E[i]=e_vec[i], XX[i] = e_sol[i] );end do:Note that the vectors NiMmJSNYWEc2IyIiIg== and NiMmJSNYWEc2IyIiIw== are the basis vectors for the linear combination that is the general solution of this homogeneous system.

Returning to the system introduced at the beginning of this lesson, namely,EQ1 := select(has,sys1,diff(x(t),t))[1]: EQ2 := select(has,sys1,diff(y(t),t))[1]: Sys1 := [EQ1,EQ2];recall that this system is non-homogeneous. In fact, the coefficient matrix and forcing term are(Aneg,b) := GenerateMatrix( eval(Sys1,[x(t)=_x,y(t)=_y]), [_x,_y] ):A := -Aneg;As seen previously, the general solution of this system isgen_soln := dsolve( sys1, {x(t),y(t)} );which, when written in vector form and factored, appears asX_soln := eval( X, gen_soln );Since this system is non-homogeneous, the particular solution is non-trivial. In this case, one particular solution isXp := eval( X_soln, [_C1=0,_C2=0] );and a basis for the solution of the corresponding homogeneous system isX1 := map( coeff, X_soln, _C1 );X2 := map( coeff, X_soln, _C2 );This leads to the following vector form for the general solution of this systemevalm(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, considerL,V := Eigenvectors( A );The eigenvalues, and hence the eigenvalues, appear in complex conjugate pairs. These complex conjugate eigenpairs and the products NiMtJSRleHBHNiMqJiYlJ2xhbWJkYUc2IyUiaUciIiIlInRHRis= NiMmJSJFRzYjJSJpRw== arefor i from 1 to Dimension(L) do e_val[i] := L[i]; e_vec[i] := Column(V,i); e_sol[i] := exp( e_val[i]*t ) * e_vec[i]; print( lambda[i]=e_val[i], E[i]=e_vec[i], XX[i]=e_sol[i] );od:The vectors NiMmJSNYWEc2IyIiIg== and NiMmJSNYWEc2IyIiIw== 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 NiMlJmFscGhhRw== and NiMlJWJldGFH are real numbers, then NiMvLSUkZXhwRzYjLCYlJmFscGhhRyIiIiomJSViZXRhR0YpJSJJR0YpRikqJi1GJTYjRihGKSwmLSUkY29zRzYjRitGKSomRixGKS0lJHNpbkdGM0YpRilGKQ==This is how the complex exponentials in NiMmJSNYWEc2IyIiIg== and NiMmJSNYWEc2IyIiIw== can be simplified. However, note that the sum and difference of the complex-conjugate pairNiMvJiUiekc2IyIiIiwmJSJhR0YnKiYlImJHRiclImlHRidGJw== NiMvJiUiekc2IyIiIywmJSJhRyIiIiomJSJiR0YqJSJpR0YqISIi are respectively, NiMvLCYmJSJ6RzYjIiIiRigmRiY2IyIiI0YoKiZGK0YoJSJhR0Yo NiMvLCYmJSJ6RzYjIiIiRigmRiY2IyIiIyEiIiooRitGKCUiaUdGKCUiYkdGKA== Hence, the real and imaginary parts of NiMmJSJ6RzYjIiIi are respectivelyNiMvJSJhRyomLCYmJSJ6RzYjIiIiRiomRig2IyIiI0YqRipGLSEiIg== NiMvJSJiRyooLCYmJSJ6RzYjIiIiRiomRig2IyIiIyEiIkYqRi1GLiUiaUdGLg== Linear combinations of solutions to a linear equation are again solutions of the equation. This is the principle of superposition, which says nothing more than that if U and V are solutions of X' = A X, then W = NiMlImFH U + NiMlImJH V is also a solution, as the followng calculation verifies.W' = (NiMlImFH U + NiMlImJH V)' = NiMlImFHU' + NiMlImJHV' = NiMlImFH A U + NiMlImJH A V = A (NiMlImFH U + NiMlImJH V) = A W (See also Section 19.B in Lesson 19.)Thus, the real and imaginary parts of a single complex solution are themselves distinct (real) solutions. The real and imaginary parts of NiMmJSNYWEc2IyIiIg== can be obtained viaU := map( evalc@Re, e_sol[1] );V := map( evalc@Im, e_sol[1] );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.)

For a final example, consider the homogeneous systemsys3 := [ diff( x(t), t ) = x(t) - 2*y(t), diff( y(t), t ) = y(t) ];which has coefficient matrix and forcing vector(Aneg,b) := GenerateMatrix( eval(sys3,[x(t)=_x,y(t)=_y]), [_x,_y] ):A := -Aneg;The general solution to the system, as reported by dsolve, isinfolevel[dsolve] := 3:gen_soln := dsolve( sys3, {x(t),y(t)} );infolevel[dsolve] := 0:A particular solution to the system is the trivial solutionX_soln := eval( X, gen_soln ):Xp := eval( 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 functionsX1 := 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 formX = _C1 * X1 + _C2 * X2 + Xp;Notice that both terms in the homogeneous solution involve the same exponential, NiMtJSRleHBHNiMlInRH. After factoring the exponential, one of the homogeneous terms has a coefficient that is not constant. It is a function of NiMlInRH. These features will have to be explained during the eigenvalue decomposition which starts as usualL,V := Eigenvectors( A );It is not difficult to see that NiMvSSdsYW1iZGFHNiIiIiI= is an eigenvalue of the coefficient matrix. The fact that this eigenvalue appears twice in the list of eigenvalues means its algebraic multiplicity of NiMvSSdsYW1iZGFHNiIiIiI= is 2. The geometric multiplicity is only 1 because there is only one linearly dependent eigenvector.In more complicated situations it is often easier to work with this information in a different form.e_decomp := Eigenvectors( A, output=list );The new form collects all information for each eigenvalue in a separate list. Each sublist containing three entries. The first entry in each sublist is an eigenvalue. The second entry is the algebraic multiplicity of the corresponding eigenvalue. The algebraic multiplicity is the number of times the eigenvalue is a root of the characteristic equation. The third entry in each sublist is a set. This set contains all eigenvectors that belong to the eigenvalue at the beginning of the corresponding list. The number of eigenvectors that correspond to an eigenvalue is called the geometric multiplicity of the eigenvalue.Extraction of the one eigenpair is implemented viafor i from 1 to nops(e_decomp) do e_val[i] := e_decomp[i][1]; e_vec[i] := e_decomp[i][3][1]; e_sol[i] := exp( e_val[i]*t ) * e_vec[i]; print( lambda[i]=e_val[i], E[i]=e_vec[i], XX[i]=e_sol[i] );od:Because NiMvJSdsYW1iZGFHIiIi is an eigenvalue with algebraic multiplicity 2 and geometric multiplicity 1, the eigenvalue decomposition yields only one solution to the homogeneous equation. This solution, namely,XX[1] = 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 formNiMmJSNYWEc2IyIiIw== = NiMlInRH NiMmJSJYRzYjIiIi + NiMtJSRleHBHNiMqJiYlJ2xhbWJkYUc2IyIiIkYqJSJ0R0Yq NiMmJSJWRzYjIiIj and obtained in Maple asunassign('V2'); sol2_form := exp(e_val[1]*t) * (t*e_vec[1]+V2);whereV2 := <x2,y2>;That is, assumer1 := equate( X, sol2_form );Now, substitution of the proposed solution into the (homogeneous) system of ODEs yieldssol2_requires := eval( sys3, r1 );Note that the second condition is trivially satisfied for all values of y2. However, the first condition is satisfied only whensol2_satisfied := solve( identity(sol2_requires[1],t), {x2,y2} );This leads to the one-parameter family of solutionssol2_family := eval( sol2_form, sol2_satisfied );Notice that the component of the solution that depends on the remaining parametermap( coeff, sol2_family, x2 );is a multiple of the first solution to the homogeneous systeme_sol[1];and so is not needed again in the basis. Therefore, the second basis solution for the homogeneous solution is chosen to bee_sol[2] := eval( sol2_family, x2=0 );To summarize, the basis of homogeneous solutions found by this method is{ e_sol[1], e_sol[2] };and the basis of solutions found by inspection of the dsolve solution is{ X1, 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 NiMtJSJ5RzYjJSJ0Rw== is independent of NiMtJSJ4RzYjJSJ0Rw==, as seen from ode_x,ode_y := selectremove( has, sys3, diff(x(t),t) );Thus, the second equation can be solved for NiMtJSJ5RzYjJSJ0Rw== without knowledge of NiMtJSJ4RzYjJSJ0Rw==. This gives sol_y := dsolve( ode_y[1], y(t), [linear] );Now, since NiMtJSJ5RzYjJSJ0Rw== is known, this result can be substituted into the ODE for NiMtJSJ4RzYjJSJ0Rw==, givingode_x2 := eval( ode_x, sol_y );and solved for NiMtJSJ4RzYjJSJ0Rw==, givingsol_x := dsolve( ode_x2[1], x(t) );The result is these calculations issubs(sol_x,sol_y, X);which is equivalent to any of the solutions obtained above.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 problemIVP;a procedure for the numerical approximation to the solution via Euler's method is obtained withsoln_euler := dsolve( IVP, [x(t),y(t)], type=numeric, method=classical );The approximate solution at NiMvJSJ0RyIiIg== issoln_euler(1);The odeplot command can be used to create a phase portrait for the solution, and results in Figure 11.12.odeplot( soln_euler, [x(t),y(t)], 0..4, title="Figure 11.12" );Alternatively, a plot of the individual components of the solution can also be obtained with the odeplot command, as shown in Figure 11.13.odeplot( soln_euler, [[t,x(t)],[t,y(t)]], 0..4, labels=[`t`,``], legend=[`x(t)`,`y(t)`], title="Figure 11.13" );See also the discussion in Section 1.C for additional details and options for working with Maple-generated numerical solutions to an IVP.

[Back to ODE Powertool Table of Contents]