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
| > | 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
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 ); |
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:
| > | 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); |
| > | tydatalist := [ seq( [tdata[i],ydata[i]], i=0..n ) ]: plot( tydatalist, `t`=0..1.2, `y`=-0.5..1.6, color='green'); |
and phase plots like this:
| > | xydatalist := [ seq( [xdata[i],ydata[i]], i=0..n ) ]: plot( xydatalist , color='blue'); |
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); |
| > |
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.