• Outline of Lesson 11
• Initialization
• 11.A Equilibrium Analysis
  • 11.A-1 Nullclines
  • 11.A-2 Equilibrium Solutions
• 11.B Graphical Analysis
  • 11.B-1 Direction Fields
  • 11.B-2 Phase Portrait
  • 11.B-3 Solution Curves
• 11.C Analytic Solutions
  • 11.C-1 One-Step Solutions using
  • 11.C-2 Eigenvalue Analysis
    • Example 1: Real and Distinct Eigenvalues
    • Example 2: Complex Eigenvalues
    • Example 3: Repeated Eigenvalues
• 11.D Numerical Solutions

ORDINARY DIFFERENTIAL EQUATIONS POWERTOOL

Lesson 11 -- First-Order Linear Systems

Prof. Douglas B. Meade

Industrial Mathematics Institute

Department of Mathematics

University of South Carolina

Columbia, SC 29208

URL:   http://www.math.sc.edu/~meade/

E-mail:  meade@math.sc.edu

Copyright  2001  by Douglas B. Meade

All rights reserved

-------------------------------------------------------------------

Outline of Lesson 11

11.A Equilibrium Analysis

                11.A-1 Nullclines

                11.A-2 Equilibrium Solutions

11.B Graphical Analysis

                11.B-1 Direction Fields

                11.B-2 Phase Portraits

                11.B-3 Solution Curves

11.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 Eigenvalues

11.D Numerical Solutions

>

Initialization

>

restart;

>	with( DEtools ):

>	with( plots ):

>	with( LinearAlgebra ):

>	with( student ):

Warning, the name changecoords has been redefined

>

11.A Equilibrium Analysis

11.A-1 Nullclines

In the -plane, an isocline for the differential equation is a curve along which has a constant value.  Thus, an isocline is a level curve of the function .

If the contstant value of 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 ( or ) does not change, i.e., where or .

The nullclines for the system

are found by graphing the curves defined implicitly by the equations and .  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,

which are rendered in Maple via

>	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 the straight lines whose equations 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, 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 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 );

>

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");

>

11.A-2 Equilibrium Solutions

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 nullclines

>	nullcline1;

>

the intersection of the nullclines can be approximated from Figure 11.1, 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, 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");

>

11.B Graphical Analysis

11.B-1 Direction Fields

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 variable .  The specific interval chosen is irrelevant, provided the endpoints are numeric.  Hence, pick the domain

>	DOMAIN := 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 ). 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.

>

11.B-2 Phase Portrait

A phase portrait for an autonomous system of ODEs displays solutions to the system in "phase space" where solutions and are treated as parametric equations for the curve . (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 and , 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;

>

for the -axis,

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

>

for the -axis, and

>	IC := 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 domain

>

DOMAIN;

>

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 as

>

WINDOW;

>

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

>

11.B-3 Solution Curves

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 is

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

>

Figure 11.10, containing graphs of the - and -components of the particular solution satisfying this initial condition, is created and displayed with

>	plot_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 in

>	IC2 := [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 with

>	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] ):

>

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 whereas the blue curve on the right is the graph of both of which correspond to the first initial condition.  The pair of green curves corresponds to the second initial condition.

>

11.C Analytic Solutions

11.C-1 One-Step Solutions using dsolve

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)} );

>

more easily grasped if written as

>	xt = eval(x(t), gen_soln); yt = eval(y(t), gen_soln);

>

The solution satisfying an initial condition, say

>

IC;

>

can be found by substituting the initial condition into the general solution, thereby producing the equations

>	q1 := 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 is

>	part_soln := eval( gen_soln, q2 );

>

again more easily grasped if written as

>	xt = eval(x(t), part_soln); yt = eval(y(t), part_soln);

>

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.  The equivalence of is established with

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

>	evalb(q3);

>

whereas that of is established with

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

>	evalb(q4);

>

That these solutions satisfy the system of differential equations is established with

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

>	map( evalb, q5 );

>

and that they satisfy the initial condition is established with

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

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

>

11.C-2 Eigenvalue Analysis

If the constant-coefficient, first-order linear system of ODEs is written in the matrix form

X' = A X + b

where X = X(t) is the vector of unknown functions, X' = X is the vector of derivatives of the unknown functions, 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 .  For example, the matrix form of the system

would be

X' = =

An eigenpair for a matrix consists of a scalar and a direction (expressed by an eigenvector V) satisfying the defining equation

A V = V

Thus, the eigenvector represents a direction that is invariant under multiplication by the matrix .  At most, a vector in this direction has its length changed by the scale factor , called the eigenvalue.

If this defining equation is rewritten in the form

V = 0

where is the 2 x 2 identity matrix ,  then the eigenvalues are found as the roots of the characteristic equation, that is, the equation .

If the eigenpairs are ( , V1) and ( , V2), then the general solution of the system is given by the sum

= V1 + V2

where , are arbitrary constants.

The examples below detail how the solution of the linear system can be constructed from the eigenpairs of the system matrix .

>

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 := < x(t), y(t) >;

>

The coefficient matrix , 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 as

>	map(diff, X, t ) = A.X + 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( 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( 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 eigenvalues

>	lambda[1] = 0; lambda[2] = 2;

>

The eigenvalues are the roots of the characteristic equation

>	CharacteristicPolynomial(A,lambda) = 0;

>

which is Maple's built-in command for obtaining the equivalent of the equation .

That each eigenpair satisfied an equation of the form

A V = V

is established with the following computations:

A V1 = = = V1 = 0 =   

A V2 = = = V2 = 2 =  

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 yields

>	L,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 ( ) and corresponding eigenvectors ( ).  These eigenpairs, as well as the eigensolutions , are obtained and printed via

>	for 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 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 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 is

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

>

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

>	X_soln := eval( X, gen_soln );

>

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

>	Xp := eval( 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

>	L,V := Eigenvectors( A );

>

The eigenvalues, and hence the eigenvalues, appear in complex conjugate pairs.  These complex conjugate eigenpairs and the products are

>	for 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 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

 

This is how the complex exponentials in and can be simplified.  However, note that the sum and difference of the complex-conjugate pair

are respectively,

Hence, the real and imaginary parts of are respectively

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 = U + V is also a solution, as the followng calculation verifies.

W' = ( U + V)' = U' + V' =   A U + A V = A ( U + 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 can be obtained via

>	U := 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.)

>

Example 3: Repeated Eigenvalues

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

>	(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 , is

>	infolevel[dsolve] := 3:

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

>	infolevel[dsolve] := 0:

`-> Solving each unknown as a function of the next ones using the order:\
[x(t), y(t)]`
`Methods for first order ODEs:`
`--- Trying classification methods ---`
`trying a quadrature`
`trying 1st order linear`
`<- 1st order linear successful`

`Methods for first order ODEs:`
`--- 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( 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 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

>	X = _C1 * X1 + _C2 * X2 + Xp;

>

Notice that both terms in the homogeneous solution involve the same exponential, .  After factoring the exponential, one of the homogeneous terms has a coefficient that is not constant.  It is a function of .  These features will have to be explained during the eigenvalue decomposition which starts as usual

>	L,V := Eigenvectors( A );

>

It is not difficult to see that is an eigenvalue of the coefficient matrix.  The fact that this eigenvalue appears twice in the list of eigenvalues means its algebraic multiplicity of 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 via

>	for 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 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 form

= +

and obtained in Maple as

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

>

where

>	V2 := <x2,y2>;

>

That is, assume

>	r1 := equate( X, 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( 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

>

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( 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 is independent of , as seen from

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

>

Thus, the second equation can be solved for without knowledge of .  This gives

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

>

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

>	ode_x2 := eval( ode_x, sol_y );

>

and solved for , giving

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

>

The result is these calculations is

>	subs(sol_x,sol_y, X);

>

which is equivalent to any of the solutions obtained above.

>

X_soln;

>

11.D Numerical Solutions

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