306ch1.html

Chapter 1: First-order equations

Direction fields

Here's how to make a direction field for the first-order differential equation
dy/dt = 2*sin(t) - y .

> with(DEtools):

> dfieldplot( diff(y(t),t) = 2*sin(t)-y(t), y(t),
t=0..5, y=-6..6, arrows=line );

If you don't put in "arrows=line", you get arrow-heads on your mini-tangents.

>

A numerical method: Euler's

Euler's method in Maple

>    h:=0.1:          # This is our step-size.
n:=ceil(4.0/h): # This calculates how many steps
                 # of size h we need to take, to get to 4.0.
t:=0.0:
y:=0.05:         # The initial condition.
data[0]:=[t,y]: # We store the initial point so we can
                 # plot it later.
for i from 1 to n do    # This is how to make a block of
                        # code get executed n times.
   mk := y + sin(y*t):
   y := y + h*mk:
   t := t + h:
   data[i] := [t,y]: # Store the current point so we can
                      # plot it later.
od:              # This ("do" backwards) marks the end
                 # of the "do" block.
t:='t': y:='y': h:='h': # Clean up

datalist1 := [ seq( data[i], i=0..n ) ];
                 # This puts all the data into a list
                 # suitable for plotting.
plot( datalist1, 0..4, 0..8, style=POINT, color=red);
                 # This plots the (approximation to the)
                 # solution, with a horizontal (t) range
                 # from 0 to 4, and a vertical (y) range
                 # of 0 to 8. The POINT style means each
                 # data point is marked with a symbol of
                 # some kind. If "style=POINT" were omitted,
                 # the data points would be joined with
                 # straight line segments (try it!).

>

Estimating the error by repeating the computation with a smaller step size

Euler's method gives only an approximation to the solution.
To test the accuracy, we can do the computation over again,
but with a value of h that is half of the old value. Then

- the new approximate solution should be more accurate than the old approximation, and

- the difference between the two approximations can be used to estimate the error in the new approximation.

>    h:= 0.1/2:       # This is our step-size, half of the value above.
n:=ceil(4.0/h): # This calculates how many steps
                 # of size h we need to take, to get to 4.0.
t:=0.0:
y:=0.05:         # The initial condition.
data[0]:=[t,y]: # We store the initial point so we can
                 # plot it later.
for i from 1 to n do    # This is how to make a block of
                        # code get executed n times.
   mk := y + sin(y*t):
   y := y + h*mk:
   t := t + h:
   data[i] := [t,y]: # Store the current point so we can
                      # plot it later.
od:              # This ("do" backwards) marks the end
                 # of the "do" block.
t:='t': y:='y': h:='h': # Clean up

datalist2 := [ seq( data[i], i=0..n ) ];
                 # This puts all the data into a list
                 # suitable for plotting.
plot( datalist2, 0..4, 0..8, style=POINT, color=blue);
                 # This plots the (approximation to the)
                 # solution, with a horizontal (t) range
                 # from 0 to 4, and a vertical (y) range
                 # of 0 to 8. The POINT style means each
                 # data point is marked with a symbol of
                 # some kind. If "style=POINT" were omitted,
                 # the data points would be joined with
                 # straight line segments (try it!).

>    p1:=plot( datalist1 ,style=point, color=red ):
p2:=plot( datalist2 ,style=point, color=blue ):
with(plots):        # (See Section below on combining plots)
display([p1,p2]);   # (See Section below on combining plots)

Warning, the name changecoords has been redefined

>

>    errorbarset := {}:
for i from 1 to nops(datalist1) do
  xvalue := op(1,op(i,datalist1)):
  yendvalue := op(2,op(i,datalist1)):
  ycenter := op(2,op(2*i-1,datalist2)):
  width := abs(ycenter-yendvalue):
  errorbarset := errorbarset union
      {[ [xvalue,ycenter-width],
         [xvalue,ycenter+width] ]}:
od:
errorbarplot := plot(errorbarset,color=black):
with(plots):
display({p1,p2,errorbarplot});

We expect that the exact solution will pass through each of the vertical error
bars, shown here as black line segments.

Making sure the step-size h is small enough

>    print([`step-size h`,`steps n`,`approx. y(4)`]): # table header.
h:=0.1:
for j from 1 to 7 do # We will try 7 different step sizes.

   n:=ceil(4.0/h): # Compute the needed number of steps
                     # for the current step-size.
   t:=0.0:
   y:=0.05:
   for i from 1 to n do
      mk := y + sin(y*t):
      y := y + h*mk:
      t := t + h:
   od:
   print([h,n,y]):
   h := h/2.0:     # Halve the step-size.
od:

> 9.584681599-9.545526386;

> [9.584681599 - .039155213,
9.584681599 + .039155213];

Euler's method with two different step sizes,
for an example problem with an exact solution

>

and approximate the solution on the interval t = T0 = 0 to t = T = 2, by means of Euler's method.

This is done twice, for h = 0.1 and for h = 0.05.
Plots of the two approximate solutions are combined with
a plot of error bars based on comparing the approximate solutions,
and a plot of the exact solution.

>    h:=0.1:    # step size
T0:=0.0:   # initial time
y0:=0.2:   # initial y value
T:=2.0:    # time we want to get to
datalist1 := 'datalist1': datalist2 := 'datalist2':
for itry from 1 to 2 do
   t:=T0:
   y:=y0:
   data[0]:=[t,y]:
   n:=ceil((T-t)/h);
   for i from 1 to n do
      mk := y + sin(t):
      y := y + h*mk:
      t := t + h:
      data[i]:=[t,y]:
   od:
   if(itry=1) then datalist1:= [ seq( data[i], i=0..n ) ] fi;
   if(itry=2) then datalist2:= [ seq( data[i], i=0..n ) ] fi;

   h := h/2: # reduce the step size by a factor of 2
             # and start over again
od:
p1:=plot( datalist1 ,style=point, color=red ):
p2:=plot( datalist2 ,style=point, color=blue ):
errorbarset := {}:
for i from 1 to nops(datalist1) do
  xvalue := op(1,op(i,datalist1)):
  yendvalue := op(2,op(i,datalist1)):
  ycenter := op(2,op(2*i-1,datalist2)):
  width := abs(ycenter-yendvalue):
  errorbar := [[xvalue,ycenter-width],           [xvalue,ycenter+width]]:
  errorbarset := errorbarset union {errorbar}:
od:
errorbarplot := plot(errorbarset,color=black):
t:='t':      # exact solution obtained by the technique of                                                            # integrating factors; see the text
yexact := exp(t)*(y0 + int(exp(-s)*sin(s),s=0..t)):
exactsolnplot:=plot( yexact, t=0..2,color=black):
with(plots): # (See Section below on combining plots)
display([p1,p2,errorbarplot,exactsolnplot]);
h:='h': T0:='T0': t:='t': T:='T': y:='y': y0:='y0':

>

Plotting a function (such as an analytical solution)

Say you have found that

is a solution of your differential equation. You can plot the solution like this
(the same way you plot any function):

> y:= 'y': t:='t':
y := 1/(1/4 + 15*exp(-4*t)):
plot(y,t=0..3);

>

What is Maple doing to plot a function?

When asked to plot a function as above, Maple simply selects a large number of values of the

independent variable (t in the above example), evaluates the function at each of them, and then joins

the corresponding points with straight lines. You can see the points it's choosing by selecting

the "point" style:

> plot(y,t=0..3,style=point);

Usually it chooses enough points so that when it joins them with straight lines, the resulting polygonal arc

looks like a smooth curve. If it looks angular, you can increase the number of points used like this:

> plot(y,t=0..3,style=point,numpoints=200);

Ranges and labels

You can dictate a range for the vertical axis like this:

> plot(y,t=0..3,2..4);

You can add a label to the vertical axis, and a title, like this:

> plot(y,t=0..3,`hi ho`=2..4,title=`JR's plot`);

Combining plots

Frequently you'll want to superimpose several plots.

If you simply want to plot several functions for which you have formulas, you simply put the

formulas in a comma-separated list enclosed by curly braces ({}):

> plot({x^2,x^4},x=-1..1);

To combine plots which are of differing types, such as a direction field and a plot of an

analytical solution formula, you give each plot a name, and then use the "display" function

in the "plots" package.

> restart:
plot1:=plot( exp(-t),t=0..2):
with(DEtools):
plot2:=dfieldplot(diff(y(t),t)=-y(t),y(t),t=0..2,y=0..1):
with(plots):
display( [ plot1, plot2 ] );

Warning, the name changecoords has been redefined

>