• Overview
• Plotting flows in 2 diminsions
• Flow Lines in 3 dimensions
• Numerical approaches
  • Euler's Method
  • The dsolve command

Flows in Vector Fields

Worksheet by Mike May, S.J.- maymk@slu.edu

> restart: with(DEtools): with(plots):

Warning, the name changecoords has been redefined

Overview

Creating flow lines through vector fields is really solving a system of differential equations. We let the vector be the velocity vectory at a point, and want to solve for a path so that the velocity vector is the derviative of the parameterized position vector with respect to time.

This can sound intimidating, since it is really solving a system of differential equations, but it is not hard to do graphically. For flows of vector fields in two dimensions, we will use other software, but we want to show Maple for completeness. In particular, we can use Maple when we want to plot flows of vector fields in 3 dimensions, or when we want to do more complicated systems. (We used these tools again when we looked at planetary motion.)

Since we are solvng differential equations we loaded the DEtools package when we started the worksheet

>

Plotting flows in 2 diminsions

We start with a simple example. Consider the vector field [-y, x].

> fieldplot([-y, x], x=-5..5, y=-5..5);

It is clear that a flow path must swirl around the center. If we try to trace a path that is always tangent to the field lines we get circles. We would like Maple to trace out the curve for us. We need to do some set-up.

Given the vector field, [f(x,y), g(x,y)], saying that that field gives the derivative of motion would turn into x' = f(x,y) and y'= g(x, y). Since we will do paths in time, both x and y are functions of time. We now record the two differential equations.

> ODE1 := diff(x(t),t) = -y(t);
ODE2 := diff(y(t),t) = x(t);

Next we need to specify a starting position. This should be a tell us where our path starts at a given time.

> initpos1 := [x(0) = 3, y(0) =4];

Now we try to plot the vector field and curve.

> DEplot({ODE1, ODE2}, [x, y], t=0..5, [initpos1], x=-6..6, y=-6..6);

The path is turning out to be a circle, but we need more time to get all the way around. We could also add other initial positions.

> DEplot({ODE1, ODE2}, [x, y], t=0..7, [[x(0)=3, y(0)=4], [x(2)=1, y(2)=2]], x=-6..6, y=-6..6);

Two tidbits to clean things up if needed. We can suppress the arrows with the arrows=none option. We can also specify how smooth the curve is by giving the stepsize of the numerical solution. The default stepsize leads to 20 steps.

> ODE1 := diff(x(t),t) = -3*y(t);
ODE2 := diff(y(t),t) = x(t);
DEplot({ODE1, ODE2}, [x, y], t=0..7, [[x(0)=3, y(0)=4]], x=-7..7, y=-6..6, linecolor=red);
DEplot({ODE1, ODE2}, [x, y], t=0..7, [[x(0)=3, y(0)=4]], x=-7..7, y=-6..6, arrows=none, linecolor=blue, stepsize=.05);

As noted in the last worksheet, when using a new command, it is worthwhile to look at the manual page for that command.

> ?DEplot

Exercises:

1) For the vector field [y^2, 2x^2], find a flow line through (1,2).

>

2) Describe typical flow lines for the following vector fields.

(a) F(x, y) = [y, x].

>

(b) F(x, y) = [x, y].

>

(c) F(x, y) = -yi + (x + y/10)j.

>

(d) F(x, y) = (x - y) i + (x - y) j.

>

Flow Lines in 3 dimensions

If we want a flow line in a vector field of three dimensions, most things follow the obvious generalizations. The sticking point is that the DEplot3d command won't show the field lines. We can use fieldplot3d to plot the field marks, and display3d to put the two together.

>

> ODE1 := diff(x(t),t) = -3*y(t);
ODE2 := diff(y(t),t) = x(t);
ODE3 := diff(z(t),t) = 1;
vfield := [-3*y, x, 1];
initvals1 := [x(0)=3, y(0)=2, z(0)=-4];

> DEplot3d({ODE1, ODE2, ODE3}, [x,y,z], t=0..10, [initvals1],
x=-7..7, y=-7..7, z=-7..7, shading=z, stepsize=.1 );
fieldplot3d(vfield, x=-7..7, y=-7..7, z=-7..7, shading=z,
grid=[5,5,5], shading=z);
patha :=DEplot3d({ODE1, ODE2, ODE3}, [x,y,z], t=0..10, [initvals1],
x=-7..7, y=-7..7, z=-7..7, shading=z, stepsize=.1 ):
vfielda :=fieldplot3d(vfield, x=-7..7, y=-7..7, z=-7..7, shading=z,
grid=[5,5,5], shading=z):
display3d({patha, vfielda});

We see that the flow line tries to form a helix

>

Exercise:

3) Describe a flow line grenerated by the field [y-z, z-x, x-y] through the point [1, 2, 3].

>

Numerical approaches

Since the book looks at numerical computation of flow lines, we want to look at them a little bit. Maple uses advanced techniques that are beyond the scope of this course. We can however do Euler's method to find specified points.

Euler's Method

Consider our favorite vector field, [-y, x]. The flow lines should be circles around the origin.

The vector field corresponds to the system x' = -y, y' =x, where both x and y are understood to be functions of t. Euler's method uses linear approximation to say that
x[t+delt] = x[t] +(delt)*(-y[t]);
y[t+delt] = y[t] + (delt)*(x[t]);

Then we put the instructions in a loop to compute to our desired time.

> pt[0.0] := [3,4]; delt := .1;
for i from 0 to 63 do
pt[i*delt+delt] := [pt[i*delt][1]-delt*pt[i*delt][2], pt[i*delt][2]+delt*pt[i*delt][1]];
od;

>

We would like to plot the points we computed.

> pointplot([seq(pt[i*delt],i=0..63)], connect=true, view=[-7..7, -7..7]);

It becomes clear that our plot wanders off the true path. In fact Euler's method is related to integration by the left hand rule and is not very accurate unless we make delt quite small.

> pta[0.0] := [3,4]; delta := .1:
clump := 5:
truedel := delta/clump;
for i from 0 to 63 do
ptb[0] := pta[i*delt]:
for j from 1 to clump do
ptb[j] := [ptb[j-1][1]-(delt/clump)*ptb[j-1][2],
ptb[j-1][2]+(delt/clump)*ptb[j-1][1]]:
od:
pta[(i+1)*delt] := ptb[clump]:
od:
pointplot([seq(pta[i*delt],i=0..63)], connect=true, view=[-7..7, -7..7]);

This method lets us evaluate points. The point pta[3.1] should be halfway around the circle.

> pta[3.1];

Exercise:

4) How far is the point pta[6.2] from its expected value, which is back at the starting point? How does this compare to pt[6.2]?

>

The dsolve command

(Optional material - this is for completeness for the curious. It will not be on the test.)

You should be suspicious that Maple has a cleaner way of doing this. Maple uses the command dsolve.

We start by defining our DEs and initial conditions.

> DE1 := diff(x(t),t)=-y(t);
DE2 := diff(y(t),t)=x(t);
IC1 := x(0)=3;
IC2 := y(0)=4;

We now use the dsolve command to create a function, named solcurve), that is our solution curve.

> solcurve := dsolve({DE1, DE2, IC1, IC2},[x(t), y(t)], numeric);

To test the function we can evaluate at 2Pi, which should be back to our starting point.

> solcurve(2*Pi);

The answer is accurate to 8 decimal places which will be good enough for us.

You should notice that solcurve is a vector valued function giving 3 values, each of which is an equation in t, x, or y. We can use the second and third values to create a plot of the solution curve.

> pointplot([seq(
[op(2,solcurve(t*.1)[2]),op(2,solcurve(t*.1)[3])],
t=0..80)], connect=true);

>

>

This is the technique that was used in the worksheet on planetary motion to find plots of orbits starting only with the force of gravity, initial position, and initial velocity.

>