Chapter 2: Systems of two first-order equations

Direction fields and vector fields for 2D systems

Here is how to plot the vector field of the system

diff(x(t),t) = 2*x*(1-1/2*x)-x*y

diff(y(t),t) = 3*y*(1-1/3*y)-2*x*y

>    with(plots):
fieldplot( [2*x*(1-x/2)-x*y, 3*y*(1-y/3)-2*x*y],
            x=-0.10..0.5, y=2.5..3.5);

Warning, the name changecoords has been redefined

[Maple Plot]

Since it's hard to see what's going on near an equilibrium point (because the arrows

get so small), sometimes it's better to scrap the speed information, and make all the

arrows the same length. That's what we get with a direction field:

>    with(DEtools):
dfieldplot( [diff(x(t),t)=2*x(t)*(1-x(t)/2)-x(t)*y(t),
             diff(y(t),t)=3*y(t)*(1-y(t)/3)-2*x(t)*y(t)], [x,y],
             t=0..1, x=0..0.5, y=2.5..3.5 );

[Maple Plot]

Euler's method for 2D systems, combining plots

Here is how use Maple to perform Euler's method for a 2D system.
It's a straightfoward extension of the 1D version.
This example applies Euler's method to the Van der Pol initial value problem:

diff(x(t),t) = y

diff(y(t),t) = -x+(1-x^2)*y

x(0) = 1

y(0) = 1

>    deltat:=0.25:
tend:=2.0:
n:=ceil(tend/deltat): # number of steps of size deltat
t:=0.0:  
x:=1.0:  
y:=1.0:  
tdata[0] := t:  
xdata[0] := x:  
ydata[0] := y:  
for k from 0 to n-1 do     
   mk := y:     
   nk := -x + (1-x^2)*y:     
   x := x + deltat*mk:     
   y := y + deltat*nk:     
   t := t + deltat:     
   tdata[k+1] := t:   
   xdata[k+1] := x:  
   ydata[k+1] := y:  
od:  
t:='t': x:='x': y:='y':                      # (clean-up)
printf("%s   %s   %s   %s\n",
       "i","tdata[i]","xdata[i]","ydata[i]");
for i from 0 to n do
  printf("%d %8.3f %12g %12g\n",
          i,tdata[i],xdata[i],ydata[i])
od:

i   tdata[i]   xdata[i]   ydata[i]

0    0.000     1.0          1.0     

1     .250     1.250         .750   

2     .500     1.43750       .332031

3     .750     1.520508     -.115864

4    1.000     1.491542     -.457989

5    1.250     1.377045     -.69065 

6    1.500     1.204382     -.880162

7    1.750      .984342    -1.082121

8    2.000      .713811    -1.336613

>   

You can make time-plots of the dependent variables like this:

>    txdatalist := [ seq( [tdata[i],xdata[i]], i=0..n ) ];
plot( txdatalist, `t`=0..1.2, `x`=-0.5..1.6);

txdatalist := [[0., 1.0], [.25, 1.250], [.50, 1.43750], [.75, 1.520507812], [1.00, 1.491541862], [1.25, 1.377044600], [1.50, 1.204382088], [1.75, .9843415620], [2.00, .7138112118]]
txdatalist := [[0., 1.0], [.25, 1.250], [.50, 1.43750], [.75, 1.520507812], [1.00, 1.491541862], [1.25, 1.377044600], [1.50, 1.204382088], [1.75, .9843415620], [2.00, .7138112118]]

[Maple Plot]

>    tydatalist := [ seq( [tdata[i],ydata[i]], i=0..n ) ]:
plot( tydatalist, `t`=0..1.2, `y`=-0.5..1.6, color='green');

[Maple Plot]

and phase plots like this:

>    xydatalist := [ seq( [xdata[i],ydata[i]], i=0..n ) ]:
plot( xydatalist , color='blue');

[Maple Plot]

You can combine the Euler x,y data plot with a direction field by using "display" from the "plots" package:
plot1 will be a data structure containing the Euler x,y data plot, and
plot2 will be a be a data structure containing a direction field plot for the same Van der Pol system.

>    plot1 := plot( xydatalist, color='blue'):
with(DEtools):
plot2 := dfieldplot( [diff(x(t),t)=y(t),
             diff(y(t),t)=-x(t) + (1-x(t)^2)*y(t)], [x,y],
             t=0..1, x=-1.6..1.6, y=-1.6..1.6):
with(plots):
display(plot1,plot2);

[Maple Plot]

>   

The Euler's method solution with stepsize 0.25 doesn't  follow the direction field as closely as it could:

to improve the Euler's method solution,  change the line

"deltat := 0.25" above to read "deltat := 0.05", and re-run the blocks of code from Euler's method on down.