Doug Meade

 

Doug

---------------------------------------------------------------------
Douglas B. Meade <><
Math, USC, Columbia, SC 29208 E-mail: mailto:meade@math.sc.edu
Phone: (803) 777-6183 URL: http://www.math.sc.edu

MaplePrimes Activity


These are answers submitted by Doug Meade

I don't think I quite understand what you are trying to do. So, if what I write here is not appropriate, please refine your explanation and I'll try again.

Let's start with your sample vector (as a Vector):

restart;
ev := <1,2,3>;
                                        [1]
                                        [ ]
                                  ev := [2]
                                        [ ]
                                        [3]

An explicit way to create the parameters you indicated would be:

param1 := [ u=ev[1], v=ev[2], w=ev[3] ];
                            [u = 1, v = 2, w = 3]

Here are two alternate ways to obtain the same list of equations:

param2 := [ (u,v,w) =~ convert(ev,list)[] ];
                            [u = 1, v = 2, w = 3]
param3 := (u,v,w) =~ [convert](ev,list)[];
                            [u = 1, v = 2, w = 3]

Note that the =~ is an elementwise operator. You might think about using (u,v,w) =~ ev , but this returns a Vector of equations that is not as useful for further work as the lists produced by param1, param2, and param3.

You can put all of this in a loop as follows:

for i from 1 to 4 do
 param[i] := [ (u[i],v[i],w[i]) =~ convert(ev,list)[] ];
end do;

Of course, you will need to modify ev to reflect the 4 different eigenvectors that you have.

The assign command can be used to convert these equations into assignments, if needed:

u,v,w;
                                   u, v, w
assign( param3 );
u,v,w;
                                   1, 2, 3

I hope this has been helpful,

Doug

---------------------------------------------------------------------
Douglas B. Meade  <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu

the difference is seen even earlier.

When you start the Plot Builder with an equation, such as y=x^(2/9), the Plot Type is "Plot" (with no other options) and the only choice is to produce a "2-D implicit plot". This means Maple is going to use its implicitplot command (from the plots package, ?implicitplot ) to create the plot. While it does make sense to talk about discontinuities for an implicitly defined function, Maple's implicitplot command does not provide support for this functionality.

On the other hand, if you start Plot Builder with an expression, such as x^(2/9), there are many different types of plots that Maple thinks it might be able to produce "2-D plot", "2-D polar plot", 3-D conformal plot of a complex-valued function", "2-D conformal plot of a complex-valued f unction", "2-D complex plot", "3-D complex plot". The default is the first choice: "2-D plot", i.e., the basic plot command. The detection of discontinuities is supported within the plot command, with the discontinuities=true optional argument.

This is the technical explanation for the different experiences. As noted above, that discontinuity detection is not supported within implicitplot appears to be a weakness. (This is not an easy problem to solve, but one that would be a nice addition to Maple.)

Doug

---------------------------------------------------------------------
Douglas B. Meade  <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu

If you want do not want to do the solve in advance, just put directly in the code:

restart;
eq := x^2=16*a;
                                   2       
                                  x  = 16 a

solve( eval( eq, a=2 ), x );
                                (1/2)      (1/2)
                             4 2     , -4 2     

solve( eval( eq, a=3 ), x )[2];
                                      (1/2)
                                  -4 3     

I think you'll agree that the code is easier to read in my first approach. I'd see what Maple can do with your equation before I just assume it won't be able to find the eight solutions you expect.

Also, don't panic if you see an answer that involves RootOf. Don't be afraid to ask for more help from MaplePrimes.

Doug

---------------------------------------------------------------------
Douglas B. Meade  <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu

Here's my Maple implementation of the solution of the same problem.

restart;
Params := {
  Dosis = 5000,
  tau = 1,
  k01 = 0.016029582,
  k12 = 0.132152668,
  k21 = 0.231583052,
  V1 = 5.2794,
  c10 = 0,
  c20 = 0,
  rho = 0
  }:
F := t -> piecewise(t<tau, Dosis/tau, rho );
              /         Dosis     \
t -> piecewise|t < tau, -----, rho|
              \          tau      /
ode := 
  diff( x1(t), t ) = - (k01+k21)*x1(t) + k12*x2(t) + piecewise(t<tau, Dosis/tau, rho ),
  diff( x2(t), t ) =        k21 *x1(t) - k12*x2(t);
 d                                                    /         Dosis     \  
--- x1(t) = -(k01 + k21) x1(t) + k12 x2(t) + piecewise|t < tau, -----, rho|, 
 dt                                                   \          tau      /  

   d                               
  --- x2(t) = k21 x1(t) - k12 x2(t)
   dt                              
ic :=
  x1(0) = c10*V1,
  x2(0) = c20*V1 * k21/k12;
                                             c20 V1 k21
                     x1(0) = c10 V1, x2(0) = ----------
                                                k12    
sol := dsolve( eval( { ode, ic }, Params ), {x1(t),x2(t)} ):

The "exact" solution is, in simplified form:

simplify( sol );
 /                 /       
{ x1(t) = piecewise|t < 1, 
 \                 \       
  1250000000000    /      1      /                              (1/2)\  \
- ------------- exp|- ---------- \189882651 + 5 1357482764902521     / t|
     8014791       \  1000000000                                        /

     1250000000000    /    1      /                               (1/2)\  \
   - ------------- exp|---------- \-189882651 + 5 1357482764902521     / t|
        8014791       \1000000000                                         /

     2500000000000    14487755750000000000                  (1/2)    /
   + ------------- - ---------------------- 1357482764902521      exp|
        8014791      3626646882265280396037                          \

      1      /                               (1/2)\  \   
  ---------- \-189882651 + 5 1357482764902521     / t| + 
  1000000000                                         /   

   14487755750000000000                  (1/2)    /      1      
  ---------------------- 1357482764902521      exp|- ---------- 
  3626646882265280396037                          \  1000000000 

  /                              (1/2)\  \          
  \189882651 + 5 1357482764902521     / t|, 1 <= t, 
                                         /          
   14487755750000000000     /      1              /         
- ---------------------- exp|- ---------- (t - 1) \189882651
  3626646882265280396037    \  1000000000                   

                       (1/2)\\                 (1/2)    14487755750000000000  
   + 5 1357482764902521     /| 1357482764902521      + ---------------------- 
                             /                         3626646882265280396037 

                  (1/2)    /    1              /          
  1357482764902521      exp|---------- (t - 1) \-189882651
                           \1000000000                    

                       (1/2)\\
   + 5 1357482764902521     /|
                             /

     1250000000000    /      1              /         
   + ------------- exp|- ---------- (t - 1) \189882651
        8014791       \  1000000000                   

                       (1/2)\\
   + 5 1357482764902521     /|
                             /

     1250000000000    /    1              /          
   + ------------- exp|---------- (t - 1) \-189882651
        8014791       \1000000000                    

                       (1/2)\\    14487755750000000000                  (1/2) 
   + 5 1357482764902521     /| - ---------------------- 1357482764902521      
                             /   3626646882265280396037                       

     /    1      /                               (1/2)\  \   
  exp|---------- \-189882651 + 5 1357482764902521     / t| + 
     \1000000000                                         /   

   14487755750000000000                  (1/2)    /      1      
  ---------------------- 1357482764902521      exp|- ---------- 
  3626646882265280396037                          \  1000000000 

  /                              (1/2)\  \
  \189882651 + 5 1357482764902521     / t|
                                         /

     1250000000000    /      1      /                              (1/2)\  \
   - ------------- exp|- ---------- \189882651 + 5 1357482764902521     / t|
        8014791       \  1000000000                                        /

     1250000000000    /    1      /                               (1/2)\  \\  
   - ------------- exp|---------- \-189882651 + 5 1357482764902521     / t||, 
        8014791       \1000000000                                         //  

                   /       144739407500000000000   
  x2(t) = piecewise|t < 1, --------------------- - 
                   \          264794003528097      

   916116746675642750000000000                   (1/2)    /    1      /
  ------------------------------ 1357482764902521      exp|---------- \
  119817765346309672026096544179                          \1000000000  
                               (1/2)\  \    916116746675642750000000000   
-189882651 + 5 1357482764902521     / t| + ------------------------------ 
                                       /   119817765346309672026096544179 

                  (1/2)    /      1      /                              (1/2)\  
  1357482764902521      exp|- ---------- \189882651 + 5 1357482764902521     / t
                           \  1000000000                                        

  \
  |
  /

     72369703750000000000    /    1      /                               (1/2)
   - -------------------- exp|---------- \-189882651 + 5 1357482764902521     
       264794003528097       \1000000000                                      

  \  \   72369703750000000000    /      1      /         
  / t| - -------------------- exp|- ---------- \189882651
     /     264794003528097       \  1000000000           

                       (1/2)\  \           916116746675642750000000000   
   + 5 1357482764902521     / t|, 1 <= t, ------------------------------ 
                               /          119817765346309672026096544179 

                  (1/2)    /    1              /          
  1357482764902521      exp|---------- (t - 1) \-189882651
                           \1000000000                    

                       (1/2)\\    916116746675642750000000000      /  
   + 5 1357482764902521     /| - ------------------------------ exp|- 
                             /   119817765346309672026096544179    \  

      1              /                              (1/2)\\ 
  ---------- (t - 1) \189882651 + 5 1357482764902521     /| 
  1000000000                                              / 

                  (1/2)   72369703750000000000    /      1              
  1357482764902521      + -------------------- exp|- ---------- (t - 1) 
                            264794003528097       \  1000000000         

  /                              (1/2)\\    916116746675642750000000000   
  \189882651 + 5 1357482764902521     /| - ------------------------------ 
                                       /   119817765346309672026096544179 

                  (1/2)    /    1      /                               (1/2)\  
  1357482764902521      exp|---------- \-189882651 + 5 1357482764902521     / t
                           \1000000000                                         

  \    916116746675642750000000000                   (1/2)    /      1      
  | + ------------------------------ 1357482764902521      exp|- ---------- 
  /   119817765346309672026096544179                          \  1000000000 

  /                              (1/2)\  \   72369703750000000000    /
  \189882651 + 5 1357482764902521     / t| + -------------------- exp|
                                         /     264794003528097       \

      1              /                               (1/2)\\
  ---------- (t - 1) \-189882651 + 5 1357482764902521     /|
  1000000000                                               /

     72369703750000000000    /    1      /                               (1/2)
   - -------------------- exp|---------- \-189882651 + 5 1357482764902521     
       264794003528097       \1000000000                                      

  \  \   72369703750000000000    /      1      /         
  / t| - -------------------- exp|- ---------- \189882651
     /     264794003528097       \  1000000000           

                       (1/2)\  \\\ 
   + 5 1357482764902521     / t|| }
                               /// 

Or, with floating-point coefficients,

sol2 := evalf( simplify( sol ) );
 /                 /                                       
{ x1(t) = piecewise\t < 1., -8776.8314 exp(-0.3741028158 t)
 \                                                         

                   5                                        5                 
   - 3.031464614 10  exp(-0.005662486200 t) + 3.119232928 10 , 1. <= t, 8776.\

  8314 exp(0.3741028158 - 0.3741028158 t)

                   5                                       
   + 3.031464614 10  exp(0.005662486200 - 0.005662486200 t)

                   5                                                        \  
   - 3.031464614 10  exp(-0.005662486200 t) - 8776.8314 exp(-0.3741028158 t)/, 

                   /                      5
  x2(t) = piecewise\t < 1., 5.466113491 10 

                   5                                                          
   - 5.550121102 10  exp(-0.005662486200 t) + 8400.7612 exp(-0.3741028158 t), 

                         5                                       
  1. <= t, 5.550121102 10  exp(0.005662486200 - 0.005662486200 t)

   - 8400.7612 exp(0.3741028158 - 0.3741028158 t)

                   5                                                        \
   - 5.550121102 10  exp(-0.005662486200 t) + 8400.7612 exp(-0.3741028158 t)/

  \ 
   }
  / 

Note that Maple's solution contains terms of the form exp(a+b*t) that Mathematica wrote in the form exp(a)*exp(b*t). Note also that the MMA solution shows the two pieces of the solution in the opposite order: t>1 comes first and t<1 corresponds to the part labeled "True".

The easiest confirmation that Maple's solution is consistent with MMA's is to look at the plot of the two solutions found by Maple:

plot( eval( [x1(t)/V1,x2(t)/V1/k21*k12], sol2 union Params ), t=0..25,
      color=[blue,red], legend=[x1(t)/V1,x2(t)/(V1*k21/k12)] );

There are other ways to get the same information. If all you really need is the plot, I'd suggest avoiding the symbolic solution altogether and using the odeplot command (from the plots package, ?odeplot ).

I hope this has been helpful.

Doug

---------------------------------------------------------------------
Douglas B. Meade  <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu

 

Maple gives you only what you ask for.

Note that in your example it is more efficient to find the general solution first, and then to substitute values of a afterwards. So, I would have written your original code as follows:

eq := x^2=16*a;
                                   2       
                                  x  = 16 a
sol := solve( eq, x );
                                (1/2)      (1/2)
                             4 a     , -4 a     
for a from 1 to 3 do
  sol;
end do;
                                    4, -4
                                (1/2)      (1/2)
                             4 2     , -4 2     
                                (1/2)      (1/2)
                             4 3     , -4 3     

If all you want are the solutions when a=2, then don't ask for a=1 or a=3:

eval( sol, a=2 );
                                (1/2)      (1/2)
                             2 2     , -4 2     

And, the "second" solution when a=3:

eval( sol[2], a=3 );
                                      (1/2)
                                  -4 3   

I hope this is helpful.

Doug

---------------------------------------------------------------------
Douglas B. Meade  <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu

Here's a screen shot showing what I get with the Plot Builder in Maple 13.02.

The Find Discontinuities check box is, as others have described, in the right-hand column.

Doug

---------------------------------------------------------------------
Douglas B. Meade  <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu

I would encourage you to take advantage of the ability to use negative indices to refer to elements from the end of a list. That is, L[-1] is the last entry in L and L[-3] is the third entry from the end of L.

Here is one way to write a procedure to do what you want:

LgSm := proc( L )
  local n, S;
  S := sort(L);
  n := nops(L);
  [ seq( [S[i],S[-i]][], i=1..floor(n/2) ),
    `if`( n::odd, S[ceil(n/2)], NULL ) ];
end proc:
LgSm( [2,3,5,2,7,9] );
                              [ 2, 9, 2, 7, 3, 5 ]

The only special handling is for lists with an odd number of elements, assuming you don't want to duplicate the middle entry. (Change floor to ceil and see what changes - when the list has an odd number of elements.)

If you are not familiar with the seq command, it can be replaced by a for .. do .. end do loop. I'll leave those detail to you.

Doug

---------------------------------------------------------------------
Douglas B. Meade  <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu

Nice pictures. I think you're going to like all the options you have to make a really effective presentation of these systems.

My one word suggestion is "DEplot". This command, in the DEtools package, is very flexible - and powerful. It really is preferrable to using odeplot in this case.

Let me demonstrate with some code:

restart:
with( DEtools ):
eq1:=diff(x(t),t)=-x(t)*(1-x(t)^2-y(t)^2)^2-y(t),diff(y(t),t)=-y(t)*(1-x(t)^2-y(t)^2)^2+x(t):
ic1 := [ [0,0.68,0.68], [0,-0.68,-0.68], [0,1.5,1.5], [0,-1.5,1.5], [0,1.5,-1.5], [0,-1.5,-1.5] ]:
PP := DEplot( {eq1}, {x(t),y(t)}, t=0..20, ic1, arrows=thin, stepsize=0.1, animatecurves=true, linecolor=t ):
Warning, numpoints adjusted from 201 to 217 for animation
PP;


eq2:=diff(x(t),t)=x(t)*(1-x(t)^2-y(t)^2)-y(t),diff(y(t),t)=y(t)*(1-x(t)^2-y(t)^2)+x(t):
ic2:=[ [0,0.01,0.01], [0,-0.01,-0.01], [0,1.5,1.5], [0,-1.5,1.5], [0,1.5,-1.5], [0,-1.5,-1.5] ]:
QQ := DEplot( {eq2}, {x(t),y(t)}, t=0..20, ic2, arrows=thin, stepsize=0.1, animatecurves=true, linecolor=t ):
Warning, numpoints adjusted from 201 to 217 for animation
QQ;

I was not successful in my attempts to include an animated GIF file in my MaplePrimes post. I thought I'd seen them in other posts, but they don't seem to be working here. Bummer. So, all you get is a static image of the final frame of the animation.

While these do not use arrows to indicate the direction of flow, this information is visible in a static image through the use of t in the coloring of the solution curves. You could get fancier than this if desired, but I think this clearly shows that the convergence to the unit circle happens faster than the convergence to the origin (in the first example). Of course, maybe you can have your students (or audience) come to this conclusion on their own.

The animations I've included in this response are slightly different from the ones obtained with the above commands. I had to make some adjustments to get the filesize below the 800K maximum allowed on MaplePrimes. To be specific, I shortened the time interval to 0..10 and lowered the number of frames to 10 (from 25, the default).

You can put these animations side by side as follows:

plots:-display( Array( [PP,QQ] ) ); 

What's particularly nice about this presentation is that when you start the animation, it starts it for both frames. You can watch the solutions side by side. If you move back and forth manually, you do so in both frames. Pretty nice. (I just wish I could export that as a single graphic to upload on MaplePrimes.)

I'm not done. When I was checking the specific syntax for the DEplot command and its options, I discovered a new option: animatefield. This would be even simpler, and gives exactly the message you are trying to communicate (I think):

PP := DEplot( {eq1}, {x(t),y(t)}, t=0..10, x=-1.5..1.5, y=-1.5..1.5, arrows=thin, stepsize=0.1, animatefield=true, scaling=constrained ):
PP2;
QQ := DEplot( {eq2}, {x(t),y(t)}, t=0..10, x=-1.5..1.5, y=-1.5..1.5, arrows=thin, stepsize=0.1, animatefield=true, scaling=constrained ):
QQ2;
plots:-display( Array( [PP,QQ] ) ); 

Note the use of scaling=constrained to ensure the circular nature is evident in the animations, particularly when they are viewed side-by-side.

You might want to explore the online help for DEplot ( ?DEplot ) to see if there are other options that you can use to improve your presentation.

I've enjoyed putting this together. I wish I was teaching differential equations this semester!

Doug

---------------------------------------------------------------------
Douglas B. Meade  <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu

 

While Thomas' suggestion works, it would seem to be advisable to avoid the warning message altogether by using:

Explore( plot( a*x^2+b, x=-2..2, y=-5..5 ) );

I do note that without the explicit vertical range, the exploration is not visually appealing. All of the plots look the same, it's only the labels on the vertical axis that change. This is something that I hope would be improved in future releases of Maple.

The Explore command is an underused gem. It does have its limitations, but it's still quite powerful.

Doug

---------------------------------------------------------------------
Douglas B. Meade  <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu

Or, open with a XML/text editor and directly modify the SHIFT-RETURNs. Of course, be sure to save a copy of the original in case you mess something up with your first attempt.

Doug

---------------------------------------------------------------------
Douglas B. Meade  <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu

And a good one at that. Hopefully, the Documentation group will make note of this and address it in a (near) future release.

Doug

---------------------------------------------------------------------
Douglas B. Meade  <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu

Sure,

Just use display to combine the two plots, like this:

SYS := [ diff( x(t), t ) = 2*y(t),
         diff( y(t), t ) = 2*x(t)-y(t),
         diff( z(t), t ) = 1 ];
         [ d                  d                         d          ]
         [--- x(t) = 2 y(t), --- y(t) = 2 x(t) - y(t), --- z(t) = 1]
         [ dt                 dt                        dt         ]
IC := [[x(0) = 0, y(0) = .5, z(0)=0],
       [x(0) = 0, y(0) = 1, z(0)=0],
       [x(0) = 0, y(0) = -.5, z(0)=0],
       [x(0) = 0, y(0) = -1, z(0)=0]];
    [[x(0) = 0, y(0) = 0.5, z(0) = 0], [x(0) = 0, y(0) = 1, z(0) = 0], 

      [x(0) = 0, y(0) = -0.5, z(0) = 0], [x(0) = 0, y(0) = -1, z(0) = 0]]
RHS := eval( map( rhs, SYS ), [x(t)=x,y(t)=y,z(t)=z] );
                              [2 y, 2 x - y, 1]
with( plots ):
p1 := fieldplot3d( RHS, x=-1..1, y=-1..1, z=-1..1, grid=[5,5,5], axes=boxed ):
p2 := DEplot3d(SYS, [x(t), y(t),z(t)], t = -1 .. 1, IC, stepsize = .2, linecolor=blue ):
display( p1, p2, view=[-1..1,-1..1,-1..1] );

It would be nicer if this was more convenient, but this is not so bad.

Doug

---------------------------------------------------------------------
Douglas B. Meade  <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu

OK. Continuing with my previous example, here's how to access the eigenvalues from the vector:

Lambda[1]-Lambda[2];
                                       (1/2)
                                 I 1223     

You will note that when I called the Eigenvectors command I used a "multiple assignment" (because I knew the result would be two objects (eigenvalues and eigenvectors).

If I had not done this and, instead, given this compound object a single name, then I would have disected the result as follows:

EV := Eigenvectors( A );
          [  5   1       (1/2)]  [        87                  87        ]
          [- - + - I 1223     ]  [------------------  ------------------]
          [  2   2            ]  [71   1       (1/2)  71   1       (1/2)]
    EV := [                   ], [-- + - I 1223       -- - - I 1223     ]
          [  5   1       (1/2)]  [2    2              2    2            ]
          [- - - - I 1223     ]  [                                      ]
          [  2   2            ]  [        1                   1         ]

lambda[1] := EV[1][1];
lambda[2] := EV[1][2];
v[1] := Column( EV[2], 1 );
v[2] := Column( EV[2], 2 );
                                     5   1       (1/2)
                      lambda[1] := - - + - I 1223     
                                     2   2            
                                     5   1       (1/2)
                      lambda[2] := - - - - I 1223     
                                     2   2            
                                [        87        ]
                                [------------------]
                                [71   1       (1/2)]
                        v[1] := [-- + - I 1223     ]
                                [2    2            ]
                                [                  ]
                                [        1         ]

                                [        87        ]
                                [------------------]
                                [71   1       (1/2)]
                        v[2] := [-- - - I 1223     ]
                                [2    2            ]
                                [                  ]
                                [        1         ]
lambda[1]-lambda[2];
                                       (1/2)
                                 I 1223     

Does this get to the root of your problem?

Doug

---------------------------------------------------------------------
Douglas B. Meade  <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu

 

DEplot3d does not provide direct support for plotting the direction field. (It seems like it should, and maybe this can be looked at for a future release.)

But, the fieldplot3d command (in the plots package), see ?fieldplot3d , can be used to plot the direction field for a system of 3 ODEs.

SYS := [ diff( x(t), t ) = 2*y(t),
         diff( y(t), t ) = 2*x(t)-y(t),
         diff( z(t), t ) = 1 ];

         [ d                  d                         d          ]
         [--- x(t) = 2 y(t), --- y(t) = 2 x(t) - y(t), --- z(t) = 1]
         [ dt                 dt                        dt         ]

RHS := eval( map( rhs, SYS ), [x(t)=x,y(t)=y,z(t)=z] );

                              [2 y, 2 x - y, 1]

with( plots ):
fieldplot3d( RHS, x=-1..1, y=-1..1, z=-1..1, grid=[5,5,5], axes=boxed );

Does this do what you want?

Doug

---------------------------------------------------------------------
Douglas B. Meade  <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu

This can be done with straightforward vector addition / subtraction. Here's a random example to show what I mean.

restart;
interface( prettyprint=1 ):   # to force Maple to output matrix results in a format that can be copied to MaplePrimes
                                      
with( LinearAlgebra ):
A := RandomMatrix( 2, 2 );
                                    [-38  87]
                               A := [       ]
                                    [-18  33]
Lambda,EV := Eigenvectors( A );
              [  5   1       (1/2)]  [        87                  87        ]
              [- - + - I 1223     ]  [------------------  ------------------]
              [  2   2            ]  [71   1       (1/2)  71   1       (1/2)]
Lambda, EV := [                   ], [-- + - I 1223       -- - - I 1223     ]
              [  5   1       (1/2)]  [2    2              2    2            ]
              [- - - - I 1223     ]  [                                      ]
              [  2   2            ]  [        1                   1         ]

Note that Lambda is a vector that contains the two eigenvalues of the matrix; EV is a 2x2 matrix whose columns are eigenvectors for the corresponding entry in the vector of eigenvalues.


ev1 := EV[1..-1,1];
ev2 := EV[1..-1,2];
                                [        87        ]
                                [------------------]
                                [71   1       (1/2)]
                         ev1 := [-- + - I 1223     ]
                                [2    2            ]
                                [                  ]
                                [        1         ]

                                [        87        ]
                                [------------------]
                                [71   1       (1/2)]
                         ev2 := [-- - - I 1223     ]
                                [2    2            ]
                                [                  ]
                                [        1         ]

This is one way to extract the columns of EV as the eigenvectors of A; you could also use, e.g., the Column command (Column(EV,1) ).

Now,  you just subtract or do whatever you'd normally do with these vectors.


ev2-ev1;
                  [        87                   87        ]
                  [------------------ - ------------------]
                  [71   1       (1/2)   71   1       (1/2)]
                  [-- - - I 1223        -- + - I 1223     ]
                  [2    2               2    2            ]
                  [                                       ]
                  [                   0                   ]

I hope this is helpful,

Doug

---------------------------------------------------------------------
Douglas B. Meade  <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu
First 13 14 15 16 17 18 19 Last Page 15 of 44