It's not difficult to create a general zip3 that does not use zip itself. If you want to make a truly general zip3, see the source code for zip (use showstat( zip );). But, for many purposes, the following will suffice:

zip3 := proc( f, X,Y,Z )
  local n;
  n := nops(X);
  if n <> map(nops,{Y,Z})[] then
    error "arguments 2, 3, and 4 must all have the same number of elements"
  end if;
  return seq( f(X[i],Y[i],Z[i]), i=1..n )
end proc;

Then

> zip3( (x,y,z)->[x,y,z], [1,2,3], [2,3,4], [3,3,3] );
                [1, 2, 3], [2, 3, 3], [3, 4, 3]

> zip3( (x,y,z)->[x,y,z], [1,2,3], [2,4], [3,3,3] );
Error, (in zip3) arguments 2, 3, and 4 must all have the same number of elements

I hope this is useful. Doug ----------------------------------------------------------------- Prof. Douglas B. Meade Phone: (803) 777-6183 Department of Mathematics URL: http://www.math.sc.edu/~meade/ USC, Columbia, SC 29208 E-mail: mailto:meade@math.sc.edu

It's not difficult to create a general zip3 that does not use zip itself. If you want to make a truly general zip3, see the source code for zip (use showstat( zip );). But, for many purposes, the following will suffice:

zip3 := proc( f, X,Y,Z )
  local n;
  n := nops(X);
  if n <> map(nops,{Y,Z})[] then
    error "arguments 2, 3, and 4 must all have the same number of elements"
  end if;
  return seq( f(X[i],Y[i],Z[i]), i=1..n )
end proc;

Then

> zip3( (x,y,z)->[x,y,z], [1,2,3], [2,3,4], [3,3,3] );
                [1, 2, 3], [2, 3, 3], [3, 4, 3]

> zip3( (x,y,z)->[x,y,z], [1,2,3], [2,4], [3,3,3] );
Error, (in zip3) arguments 2, 3, and 4 must all have the same number of elements

I hope this is useful. Doug ----------------------------------------------------------------- Prof. Douglas B. Meade Phone: (803) 777-6183 Department of Mathematics URL: http://www.math.sc.edu/~meade/ USC, Columbia, SC 29208 E-mail: mailto:meade@math.sc.edu

While all of the previous comments are correct, I believe they miss the true explanation of the observed behavior. While Maple is a computer algebra system capable of working with expressions, functions, and many other mathematical objects, a Maple plot is nothing more than a finite collection of points. Maple uses an adaptive routine for selecting the points to plot. This is done so that when these points are connected the result is almost always a reasonable reproduction of the original function. There are exceptions. One of my favorite examples is to plot the tangent function on the interval [0,2*Pi].

plot( tan(x), x=0..2*Pi );

Maple detects the large derivatives near x=Pi/2 and x=3*Pi/2 and selects more points close to these two points.

plot( tan(x), x=0..2*Pi, s tyle=point ); # delete the space between s and t
                                         # (workaround for MaplePrimes editor)

To obtain a plot that looks like the typical tangent function it is necessary to explicitly restrict the vertical range.

plot( tan(x), x=0..2*Pi, view=[DEFAULT,-10..10] );

It looks as though Maple includes the vertical asymptotes. In reality, this is simply the (nearly vertical) line segment connecting the last point to the left of the asymptote with the first point to the right of the asymptote. The discont=true option suppresses these segments.

plot( tan(x), x=0..2*Pi, view=[DEFAULT,-10..10], discont=true );

Back to your situation. If Maple selects x=1 (exactly 1) then the point (1,2) will be plotted. I make not promises about how the discontinuity will be handled. The probability that Maple will choose a specific value of the independent variable is effectively 0. The real resolution to this problem would appear to be a modification to plot that looks for instances where a function is defined by alternate formula at a single point and to be sure to include these points in the collection of points that are plotted. (This is effectively what Paulina suggested, but it should be automated.) I hope this long-winded explanation is of some help - and maybe encourages someone to implement this idea. Doug ----------------------------------------------------------------- Prof. Douglas B. Meade Phone: (803) 777-6183 Department of Mathematics URL: http://www.math.sc.edu/~meade/ USC, Columbia, SC 29208 E-mail: mailto:meade@math.sc.edu

Jim, I am not sure why you see 0 as a result either. I do see two reasons to not use the particular implementations that you have shown.

1. when defining a function with an arrow operator, the rhs of the -> is not evaluated at execution

   This means Maple sees the N as a local variable to the function, to have N as a true parameter in this function, you should use unapply (and not define N until after the function definition)

2. sum (and Sum) are used for symbolic sums, or sums with variable numbes of terms

   Use add for adding a finite collection of terms. In your original posts, N is fixed BEFORE defining the function so add can be used. In my recommendation, below, N is not defined at the function definition so you MUST use sum (or Sum).

Here is my recommendation for implementing your function:

restart;
U := unapply( sum(a[n](t)*sin(n*Pi*x),n=1..N), (t,x) );
N := 2;
U(t,x);
diff( U(t,x), t );
D[1](U)(t,x);

Note that it makes no difference here if you use sum or add - so long as N is finite. Doug Doug ----------------------------------------------------------------- Prof. Douglas B. Meade Phone: (803) 777-6183 Department of Mathematics URL: http://www.math.sc.edu/~meade/ USC, Columbia, SC 29208 E-mail: mailto:meade@math.sc.edu

Victor, You missed the emphasis in my previous post. The second argument to animate needs to be a LIST. This means it is delimted with square brackets. Change the parentheses () around the second argument to square brackets [] and it works.

with(plots);
with(DEtools);
animate(
phaseportrait, [ [(D(x))(t) = z(t), (D(z))(t) = 2*x(t)*(k-x(t)^2), (D(y))(t) = k-x(t)^2],
[x(t), y(t), z(t)], t = -6 .. 6, [[x(0) = (1+k)^(1/2), y(0) = 0, z(0) = 0]],
stepsize = 0.1, scene = ([x(t), y(t)]), scaling = constrained, linecolor = blue], k=0..2);

Interesting plots. Doug ----------------------------------------------------------------- Prof. Douglas B. Meade Phone: (803) 777-6183 Department of Mathematics URL: http://www.math.sc.edu/~meade/ USC, Columbia, SC 29208 E-mail: mailto:meade@math.sc.edu

Victor, You missed the emphasis in my previous post. The second argument to animate needs to be a LIST. This means it is delimted with square brackets. Change the parentheses () around the second argument to square brackets [] and it works.

with(plots);
with(DEtools);
animate(
phaseportrait, [ [(D(x))(t) = z(t), (D(z))(t) = 2*x(t)*(k-x(t)^2), (D(y))(t) = k-x(t)^2],
[x(t), y(t), z(t)], t = -6 .. 6, [[x(0) = (1+k)^(1/2), y(0) = 0, z(0) = 0]],
stepsize = 0.1, scene = ([x(t), y(t)]), scaling = constrained, linecolor = blue], k=0..2);

Interesting plots. Doug ----------------------------------------------------------------- Prof. Douglas B. Meade Phone: (803) 777-6183 Department of Mathematics URL: http://www.math.sc.edu/~meade/ USC, Columbia, SC 29208 E-mail: mailto:meade@math.sc.edu

Here is how you could use the plots[animate] command to do what you want. I will illustrate using the more general DEplot command, but this should work for phaseportrait in the same way.

animate(
 DEplot,
 [
  [ diff(x(t),t)=x(t)*(1-y(t)), diff(y(t),t)=a*y(t)*(x(t)-1) ],
  [x(t),y(t)], t=-7..7,
  [[x(0)=1.2,y(0)=1.2],[x(0)=1,y(0)=.7]],
  stepsize=.2,
  title=`Lotka-Volterra model`,
  color=[a*y(t)*(x(t)-1),x(t)*(1-y(t)),.1],
  linecolor=t/2,
  arrows=MEDIUM, method=rkf45
 ], a=0.5..3 );

Notice that the first argument is a plot command (DEplot or phaseportrait) and the second argument is a LIST containing all of the arguments to the plot command. The animation parameter should appear somewhere in the plot commands. The third argument to animate identifies the plot parameter and the range of parameter values. The animation contains 25 frames. These plots can take a while to create, so be patient. Good luck! Doug ----------------------------------------------------------------- Prof. Douglas B. Meade Phone: (803) 777-6183 Department of Mathematics URL: http://www.math.sc.edu/~meade/ USC, Columbia, SC 29208 E-mail: mailto:meade@math.sc.edu

Here is how you could use the plots[animate] command to do what you want. I will illustrate using the more general DEplot command, but this should work for phaseportrait in the same way.

animate(
 DEplot,
 [
  [ diff(x(t),t)=x(t)*(1-y(t)), diff(y(t),t)=a*y(t)*(x(t)-1) ],
  [x(t),y(t)], t=-7..7,
  [[x(0)=1.2,y(0)=1.2],[x(0)=1,y(0)=.7]],
  stepsize=.2,
  title=`Lotka-Volterra model`,
  color=[a*y(t)*(x(t)-1),x(t)*(1-y(t)),.1],
  linecolor=t/2,
  arrows=MEDIUM, method=rkf45
 ], a=0.5..3 );

Notice that the first argument is a plot command (DEplot or phaseportrait) and the second argument is a LIST containing all of the arguments to the plot command. The animation parameter should appear somewhere in the plot commands. The third argument to animate identifies the plot parameter and the range of parameter values. The animation contains 25 frames. These plots can take a while to create, so be patient. Good luck! Doug ----------------------------------------------------------------- Prof. Douglas B. Meade Phone: (803) 777-6183 Department of Mathematics URL: http://www.math.sc.edu/~meade/ USC, Columbia, SC 29208 E-mail: mailto:meade@math.sc.edu

My initial experiences with the MapleNet server through MaplePrimes have been a little problematic. While it has - eventually - worked, the initial download have taken quite a long time. (I have used MapleNet before, and the delays I have experienced with this server are much longer than from previous uses.) This delay seems to have been related to the initial updating of Java. At the least, users should be given estimates of the time needed to do this download. An ongoing problem is the appearance of popup windows with warnings about a "General exception". The text of the warning is "Optional package installation aborted." It's not uncommon to get several of these warnings for a single MapleNet download. Clicking OK to make these disappear does not appear to have any long-term consequences - the applications do appear to work. I agree with Fred that the new file manager is a nice addition to MaplePrimes. I also agree that it's unfortunate that there is no support for deleting files. I understand how this has the potential to create dead links in posts but I also see it as a serious limitation. For example, I have several items that I would like to post but I would like to be able to do some initial experiments to find the best way to work with this new feature before I make something publicly available. Thanks for listening, Doug ----------------------------------------------------------------- Prof. Douglas B. Meade Phone: (803) 777-6183 Department of Mathematics URL: http://www.math.sc.edu/~meade/ USC, Columbia, SC 29208 E-mail: mailto:meade@math.sc.edu

I'm afraid there are no simpler ways - at present - to use Maple to construct the types of figures you desire. Why do you want to use Maple to create these figures? What do you want to do once you have the figures? Depending on your responses to these questions, you might find geometry software more suitable. Geometer's Sketchpad and Cabri are two well-established products in this area. A newer entry to this market is Geometry Expressions (http://www.geometryexpressions.com/). This is still in beta testing, but that means it is FREE. I have found the software pretty easy to use. One unique feature of Geometry Expressions is that it can produce symbolic results (GSP and Cabri are limited to numeric results.) I hope this is helpful. Doug ----------------------------------------------------------------- Prof. Douglas B. Meade Phone: (803) 777-6183 Department of Mathematics URL: http://www.math.sc.edu/~meade/ USC, Columbia, SC 29208 E-mail: mailto:meade@math.sc.edu

I'm afraid there are no simpler ways - at present - to use Maple to construct the types of figures you desire. Why do you want to use Maple to create these figures? What do you want to do once you have the figures? Depending on your responses to these questions, you might find geometry software more suitable. Geometer's Sketchpad and Cabri are two well-established products in this area. A newer entry to this market is Geometry Expressions (http://www.geometryexpressions.com/). This is still in beta testing, but that means it is FREE. I have found the software pretty easy to use. One unique feature of Geometry Expressions is that it can produce symbolic results (GSP and Cabri are limited to numeric results.) I hope this is helpful. Doug ----------------------------------------------------------------- Prof. Douglas B. Meade Phone: (803) 777-6183 Department of Mathematics URL: http://www.math.sc.edu/~meade/ USC, Columbia, SC 29208 E-mail: mailto:meade@math.sc.edu

The behavior you observed is probably explained by numerical affects with an ill-conditioned matrix. Can you tell us EXACTLY how you generated the system that gave you the bad results? Do you see this with one specific system or in general?

Here is your code modified to display the condition number (infinity-norm) and the largest component of the residual vector (||b-Ax||_infinity). I ran this 100 times for 140x140 matrices. You can see that the condition numbers are fairly large, but that the residual vectors are very reasonable. I suspect that if I created an ill-conditioned matrix I could force the computed solution to have a large component in the residual.

Notice, in particular that the ratio of the condition number to the largest element in the residual is about 10^(-20).

> restart;

> with( LinearAlgebra ):

> interface( rtablesize=150 );

> my_solve := proc(A::Matrix)

> local sz, local_A, B, sol, check;

> sz := Dimension(A);

> local_A := Matrix(A, datatype=float);

> B := Vector(1...sz[1], 1);

> sol := LinearSolve(local_A, B);

>

> check := local_A . sol - B;

> print(ConditionNumber(local_A,infinity),Norm(check,infinity));

> end proc:

> for i to 100 do
  my_solve( RandomMatrix(140, 140, generator=rand(0..1000)/1000.0) )
end do;

>

This post generated using the online HTML conversion tool
Download the original worksheet

The behavior you observed is probably explained by numerical affects with an ill-conditioned matrix. Can you tell us EXACTLY how you generated the system that gave you the bad results? Do you see this with one specific system or in general?

Here is your code modified to display the condition number (infinity-norm) and the largest component of the residual vector (||b-Ax||_infinity). I ran this 100 times for 140x140 matrices. You can see that the condition numbers are fairly large, but that the residual vectors are very reasonable. I suspect that if I created an ill-conditioned matrix I could force the computed solution to have a large component in the residual.

Notice, in particular that the ratio of the condition number to the largest element in the residual is about 10^(-20).

> restart;

> with( LinearAlgebra ):

> interface( rtablesize=150 );

> my_solve := proc(A::Matrix)

> local sz, local_A, B, sol, check;

> sz := Dimension(A);

> local_A := Matrix(A, datatype=float);

> B := Vector(1...sz[1], 1);

> sol := LinearSolve(local_A, B);

>

> check := local_A . sol - B;

> print(ConditionNumber(local_A,infinity),Norm(check,infinity));

> end proc:

> for i to 100 do
  my_solve( RandomMatrix(140, 140, generator=rand(0..1000)/1000.0) )
end do;

>

This post generated using the online HTML conversion tool
Download the original worksheet

Do the following examples illustrate the problems you are seeing?

The matrix M is the matrix in the first example in the online help for LinearSolve. Note that I have modified my_solve so that the use of LinearAlgebra is local to the procedure (just good programming style) and that the matrix is not converted to floating-point. The second change allows me to send the matrix in different formats and compare the results.

> my_solve := proc(A::Matrix)

>   local sz, local_A, B, sol;

>   use LinearAlgebra in

>     sz := Dimension(A);

> #    local_A := Matrix(A, datatype=float);

>     B := Vector(1...sz[1], 1);

>     sol := LinearSolve(A, B);

>   end use;

> end proc;

> ?LinearSolve

Exact Solution

> M := <<1,1,1,4>|<1,1,-2,1>|<3,1,1,8>|<-1,1,-1,-1>|<0,1,1,0>>:

> x := my_solve( M );

> M.x;

>

Floating-Point Solution (Digits=10)

> x_float := my_solve( evalf(M) );

> M.x_float;

> map(fnormal,%);

>

Floating-Point Solution (Digits=20)

> Digits := 20:

> x_float20 := my_solve( evalf(M) );

> M.x_float20;

> map(fnormal,%);

>

This post generated using the online HTML conversion tool
Download the original worksheet

I see nothing surprising in this. The solutions appear different because they are using different variables as the free variable. All other differences are in the use of floating-point numbers and the roundoff that this introduces.

Doug

Do the following examples illustrate the problems you are seeing?

The matrix M is the matrix in the first example in the online help for LinearSolve. Note that I have modified my_solve so that the use of LinearAlgebra is local to the procedure (just good programming style) and that the matrix is not converted to floating-point. The second change allows me to send the matrix in different formats and compare the results.

> my_solve := proc(A::Matrix)

>   local sz, local_A, B, sol;

>   use LinearAlgebra in

>     sz := Dimension(A);

> #    local_A := Matrix(A, datatype=float);

>     B := Vector(1...sz[1], 1);

>     sol := LinearSolve(A, B);

>   end use;

> end proc;

> ?LinearSolve

Exact Solution

> M := <<1,1,1,4>|<1,1,-2,1>|<3,1,1,8>|<-1,1,-1,-1>|<0,1,1,0>>:

> x := my_solve( M );

> M.x;

>

Floating-Point Solution (Digits=10)

> x_float := my_solve( evalf(M) );

> M.x_float;

> map(fnormal,%);

>

Floating-Point Solution (Digits=20)

> Digits := 20:

> x_float20 := my_solve( evalf(M) );

> M.x_float20;

> map(fnormal,%);

>

This post generated using the online HTML conversion tool
Download the original worksheet

I see nothing surprising in this. The solutions appear different because they are using different variables as the free variable. All other differences are in the use of floating-point numbers and the roundoff that this introduces.

Doug