Part III. First Order Ordinary Differential Equations
Section 1. Entering, Solving, Plotting
Here begins a more formal introduction to Maple and differential equations. If you have only skimmed over Parts I and II, please review the discussion of differential equations found on the last three pages of Part II, Section 3. In this section, three first order equations are used to illustrate some of the features of the dsolve procedure, Maple's primary tool for generating solutions to ordinary differential equations.
Occasional references to "Ledder" refer to the textbook entitled
Differential Equations: A Modeling Approach
by Glenn Ledder, McGraw-Hill, 2005.
Solve this: dsolve
It is good practice to assign a name to a differential equation before beginning an analysis. Names should also be assigned to solutions. A differential equation for y as a function of t is entered using unevaluated derivatives. All appearances of the function y must be in the form y(t). There are two ways to do this.
1. Use the derivative operator D:
> | DE := D(y)(t) = y(t)/(t^2+1);
gen_soln := dsolve(DE); |
2. Use the derivative procedure diff:
> | DE := diff(y(t),t) = y(t)/(t^2+1);
gen_soln := dsolve(DE); |
The second method will be used in this manual because the output prettyprints in a form that closely resembles what appears in textbooks. The output is the same in both cases: A solution equation if Maple can find one, and nothing at all if it cannot.
Observe that the formula for the general solution displays the integration constant as _C1. This strange choice of notation is intentional to avoid conflicts with variables that may already be assigned in the worksheet.
Checking the solution: odetest
The solution can be checked by using substitution, subs, (and simplify).
> | subs(gen_soln,DE); |
> | simplify(%); |
The output is an identity, valid for all t, signaling that the solution formula is valid over the whole t axis.
Solutions can also be checked using a procedure called odetest. The syntax to check if "soln" solves "DE" is
odetest(soln,DE)
The solution must be in the form of an equation.
Maple generates the output to odetest by making the substitution, simplifying, and then subtracting the right side from the left. Thus a correct solution formula should produce the output 0.
> | odetest(gen_soln,DE); |
Plotting solution curves
Given the general solution, the sequence operator, $, can be used to plot several solution curves. The following plot shows the solutions to DE when _C1 = -3, -2, -1, 0, 1, 2, 3.
> | plot( [rhs(gen_soln)$_C1=-3..3], t=-2..2, color=red); |
Solving an initial value problem
If the differential equation is accompanied by an initial condition of the form
we have an initial value problem (IVP for short). Look for a solution by applying dsolve to a set containing DE and the initial condition:
dsolve( {DE, y(t0)=y0} )
For example, the solution to DE satisfying the initial condition y(1) = 3 is obtained as follows.
> | soln := dsolve( {DE, y(1)=3} ); |
Use eval to check the initial value.
> | eval(soln,t=1); |
The following picture shows the seven previous solutions, the IVP solution (dotted), and the initial point.
> | plot( [rhs(gen_soln)$_C1=-3..3, rhs(soln), [[1,3]]], t=-2..2, -4..4,
color=[red$7,blue,black], linestyle=[SOLID$7,DOT], style=[line$8,point]); |
Getting information out of a solution formula
Suppose, for example, we want to know the value of t for which the IVP solution has the value 1. fsolve should be able to find it.
> | fsolve( rhs(soln)=1, {t} ); |
Note that the unknown in fsolve, t, is placed inside set brackets. This generates the solution in the form of an equation that is used below in eval to verify that the value is correct. The input "check_value" is just to make a nice looking output.
> | check_value, eval(soln,%); |
This is not quite what we need. Apply evalf.
> | evalf(%); |
A solution formula in terms of the initial condition
A solution formula in terms of a generic initial condition, say, y(t0) = y0, can be obtained as follows.
> | IVP_soln := dsolve( {DE, y(t0)=y0} ); |
The unapply procedure turns the formula into a function phi of t, t0, and y0.
> | phi := unapply(rhs(IVP_soln),t,t0,y0); |
And the phi function can be used to create pictures of solutions satisfying a specific set of initial conditions.
Problem: Plot solutions starting at points in the t,y plane evenly spaced around the unit circle.
> | start_points := [[cos(n*Pi/4),sin(n*Pi/4)]$n=0..7]:
plot( [[t,phi(t,cos(k*Pi/4),sin(k*Pi/4)),t=cos(k*Pi/4)..3]$k=0..7, start_points], t=-3..3, -3..3, axes=boxed, color=[red$8,black$8], style=[line$8,point$8]); |
Now run time backwards from the same initial points.
> | plot( [[t,phi(t,cos(k*Pi/4),sin(k*Pi/4)),t=-3..cos(k*Pi/4)]$k=0..7,
start_points], t=-3..3, -3..3, axes=boxed, color=[red$8,black$8], style=[line$8,point$8]); |
Look before you leap
Before solving any first order differential equation it is a very good idea to put it into normal form
and think a bit about the function on the right hand side:
f (t, y ) .
That way unusual Maple output can be anticipated and pitfalls can be avoided. The essential facts to keep in mind are given in the following theorem.
Unique solution theorem Unique IVP solutions are guaranteed to exist inside of t,y rectangles where f (t, y) and
are both continuous. (Ledder, Section 2.4, Theorem 2.4.1)
Here is an example (very similar to DE above).
> | DE2 := diff(y(t),t) = y(t)/(t^2-1); |
The differential equation is already in normal form. Note that the function on the right side
is continuous and has a continuous y partial derivative in any rectangle that does not contain any points on the vertical line t = 1 or the vertical line t = -1.
Let's see Maple's general solution.
> | gen_soln2 := dsolve( DE2 ); |
Examination of the formula reveals that when _C1 is not 0, it is only valid for t between -1 and 1, unless the constant is allowed to take on imaginary values. Several solution curves are plotted below using an input similar to the one used to generate the first plot in this section.
> | plot( [rhs(gen_soln2)$_C1=-3..3], t=-2..2, color=red); |
It appears that solutions have a vertical asymptote at t = -1 and they disappear at t = 1. The next plot gets closer by restricting the horizontal range to -1..2 and the vertical range to -10..10.
> | plot( [rhs(gen_soln2)$_C1=-3..3], t=-1..2, -10..10, color=red); |
From this picture we conjecture that all solutions to DE2 that are valid on the open interval (-1,1) have a vertical asymptote at the left endpoint and will approach 0 as t approaches the right endpoint (except the zero solution which goes on forever in both directions).
The question is "How does Maple represent solutions when t is greater than 1 or when t is less than -1?" The answer, as illustrated below, is rather unusual (and probably not the one you might guess).
Question: What is the solution to DE2 satisfying the initial condition y(2) = 5?
Answer: See Maple's solution.
> | soln3 := dsolve( {DE2, y(2)=5} ); |
The same formula can be used, but with an imaginary constant:
= 3 i. See the graph below.
> | plot( rhs(soln3), t=0..4, -10..10, color=red); |
The next picture shows several solutions on the interval t > 1. The values for _C1 are listed first.
> | imag_const := [4*k*I$k=-3..3];
plot( [seq(rhs(gen_soln2), _C1=imag_const)], t=0..4, -10..10, color=red); |
Based upon this picture we conjecture that all solutions on the interval t > 1 start at (1,0) with a vertical tangent and have a horizontal asymptote (with the exception of the zero solution).
Note: The unique solution theorem (page 36) tells us that there is a solution starting at any point (a,b) where a > 1. Maple's solution formula, which is in terms of a complex constant, is not the one you would come up with "by hand". Try your algebra skills on Maple's solution (soln3) and see if you can make it into a formula that is real when t is greater than one. Hint. Factor
and then write the square root in the numerator as the product of two terms.
In the section exercises, you will find a problem asking for a conjecture about the nature of the solutions on the interval t < -1.
To obtain implicit solutions, add the keyword implicit
As a rule, dsolve attempts to return explicit solution formulas even when there are more than one of them.
> | DE3 := diff(y(t),t) = t/y(t);
gen_soln3 := dsolve( DE3 ); |
When an initial condition is added, dsolve picks the appropriate formula.
> | soln3 := dsolve( {DE3, y(0)=1} );
plot( rhs(soln3), t=-3..3, 0..4); |
This solution curve appears to be branch of a hyperbola. Several branches can be plotted in two ways.
Method 1 Plot solutions for a selected set of initial conditions.
A phi function can be constructed and you are asked to do so in the exercises. However, it is just as easy to simply make a list of the initial conditions and then use the sequence procedure to make solution formulas.
> | inits := [y(0)=1,y(0)=-1,y(1)=0.1,y(1)=-0.1,y(-1)=0.1,y(-1)=-0.1]:
solns := seq( rhs(dsolve( {DE3,IC} ) ), IC=inits): plot( [ solns], t=-3..3, y=-3..3, color=[red$2,blue$4], title="Six Solution Curves"); |
Method 2 Ask Maple for an implicit solution formula, replace y(t) with y, drop the constant of integration, convert it into a function F of t and y, and then apply contourplot to F to obtain an entire family of solution curves.
> | implicit_soln := dsolve( DE3, implicit); |
> | F := subs(y(t)=y,_C1=0,lhs(implicit_soln)); |
The contourplot procedure is in the plots package. The first picture below uses the optional equation
contours = [1,-1]
to obtain the same solution curves that are displayed above using plot. All the curves are colored blue.
> | with(plots):
contourplot( F, t=-3..3, y=-3..3, contours=[1,-1], color=blue, title="Solutions as Level Curves"); |
The last picture illustrates how to get a fancier picture with more contours (the default number) and background coloring varying from blue to red as the values for the contour lines increase from smallest to largest.
> | contourplot( F, t=-3..3, y=-3..3, filled=true, coloring=[blue,red],
title="Solutions as Level Curves"); |
Go to: Next Section Previous Section