306ch2.mws

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) = x-y

Diff(y,t) = x+y

> with(plots):
fieldplot( [x-y,x+y], x=-1..1, y=-1..1);

[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. TO do this, we just divide the vector field by its own length:
That's what we get with a direction field:

> f1:=x-y:
f2:=x+y:
len:=sqrt(f1^2+f2^2):
fieldplot( [f1/len,f2/len], x=-1..1,y=-1..1);

[Maple Plot]

Euler's method for 2D systems

Here is how use Maple to perform Euler's method for a 2D system.
It's a straightfoward extension of the 1D version.
I'm applying it 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

> dt:=0.25:
n:=ceil(1.0/dt):
t:=0.0:
x:=1.0:
y:=1.0:
txdata[0] := [t,x]:
tydata[0] := [t,y]:
xydata[0] := [x,y]:
for i from 1 to n do
mk := y:
nk := -x + (1-x^2)*y:
x := x + dt*mk:
y := y + dt*nk:
t := t + dt:
txdata[i] := [t,x]:
tydata[i] := [t,y]:
xydata[i] := [x,y]:
od:
t:='t': x:='x': y:='y': # (clean-up)
txdatalist := [ seq( txdata[i], i=0..n ) ]:
tydatalist := [ seq( tydata[i], i=0..n ) ]:
xydatalist := [ seq( xydata[i], i=0..n ) ]:

If we don't suppress the output on the last three lines, we'll see

txdatalist := [[0, 1.0], [.25, 1.250], [.50, 1.4375...

tydatalist := [[0, 1.0], [.25, .750], [.50, .332031...

xydatalist := [[1.0, 1.0], [1.250, .750], [1.43750,...

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

> plot( {txdatalist,tydatalist},
`t`=0..1.2, `x and y`=-0.5..1.6);

and phase plots like this:

> plot( xydatalist );

You can combine the Euler x,y data plot with a direction field using "display" from the "plots" package, as before

(see Ch1 guide).

Plotting parametrized curves.

Here is the syntax for plotting a parametrized curve:

> x:=exp(-3*t);
y:=exp(-t);
plot( [x,y,t=0..2] );

x := exp(-3*t)

y := exp(-t)

This form of the plot command is also useful for plotting a function y(t) when you actually have t(y):

> restart:
t:=9*y+8*cos(y);
plot([t,y,y=0..3]);

t := 9*y+8*cos(y)

Note that y is on the vertical axis - where we want it!

>