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

Jan, To reproduce Maple's computations for evalf( sin(Pi/8), 2 ); on a calculator you would have to know exactly how the sine function is being evaluated. Remember: in the above command, ALL numerical computations are being done with only two significant digits. The first comment about fixed-point computations is that numbers are represented in terms of two INTEGERS. For example, 123.45 = 12345x10^(-3) is represented by the pair [12345, -3]. All calculations are performed with this data structure. That is, all operations are done in terms of integers. For this reason it's easier to work with polynomials or rational functions than with transcendental functions. I think the following will give you an idea about what is happening. Approximate the sine function with it's two-term Maclaurin approximation:
F := x -> x - x^3/3!;
This gives
evalf[2]( F(Pi/8) );
                             0.36
evalf[3]( F(Pi/8) );
                             0.382
Before going further, I should also mention that the number of digits used in Maple's fixed-point calculations can be set via the Digits variable. (The default is 10 digits.)
Digits;
                             10
Digits := 2:
x := evalf(Pi/8);
                              0.37
x^3;
                             0.051
x^3/3!;
                             0.0085
x-x^3/3!;
                              0.36
The calculation of x^3 is basically, 0.37^3 = (37 x 10^(-2))^3 = 37^3 x 10^(-6) = 50563 / 1000000 = 51000 x 10^(-6) = 51 x 10^(-3) = 0.051. Notice that rounding (to 2 significant digits) is used in the integer division. Then, x^3/3! = (51 x 10^-3) / (6 x 10^0) = (51/6) x 10^(-3-0) = 85 x 10^(-4) = 0.0085. Now, x-x^3/3 = 37 x 10^(-2) - 85 x 10^(-4) = (3700 - 85) x 10^(-4) = 3615 x 10^(-4) = 3600 x 10^(-4) = 0.36. Here are Maple's results with 3 digits:
Digits := 3:
x := evalf(Pi/8);
                             0.392
x^3;
                             0.0602
x^3/3!;
                             0.0100
x-x^3/3!;
                             0.382
Fixed-point arithmetic is NOT an easy topic to understand. I am not an expert. I hope I have given a decent introduction to some of the basics. A google search should turn up some additional references. 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/~meade/
<\pre>
Jan, The "problem" you report is a very common misunderstanding. The second argument in evalf is the number of digits used in all floating-point arithmetic. It does NOT guarantee that the results are accurate to this many digits. Let me illustrate:
q := sin(Pi/8);
                              /1   \

                           sin|- Pi|

                              \8   /
qf := evalf( q );
                          0.3826834325
evalf( q , 2 );
                              0.36
evalf( qf, 2 );
                              0.38
evalf( q, 3 );
                             0.382
evalf( qf, 3 );
                             0.383
If your goal is to get answers accurate to a certain number of digits it is necessary to do the calculations with more digits. Exactly how many depends on the specific problem. Here is a quick procedure that shows how you could get the type of responses that you appear to be seeking:
MyRound := ( e, n::nonnegint )
           -> round( evalf[n+2]( 10^n*e ) )/10.^n;
Notice that MyRound does all calculations with 2 more digits than are reported in the result. If you do not like the trailing zeroes on the final result, you might want to look at using the printf command, e.g.,
printf( "%4.2f", MyRound(q,2) );
I hope this is useful, 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/~meade/
<\pre>
Emrah, You should look at the permute command in the combinat package. Here is a possible solution to your request:
P := ( n::integer, symbols::list )
      -> combinat[permute]( [seq( s$n, s=symbols )], n );
Then, your example would be produced with:
P( 3, [0,1] );
You could also do something like
P( 5, [heads,tails] );
The only real "trick" in this solution is creating the list of symbols with n copies of each symbol. This is what is done with the seq command above. If you want to use the permute command elsewhere in a worksheet I recommend loading the combinat package by issuing the command
with( combinat );
Then you can refer to permute without explicit reference to combinat[permute]. I hope this is useful. 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/~meade/
<\pre>
I do not believe I have an "answer" for you, but I do have more questions that might help you find what you seek. 1. I believe there is an error in the definition of pi[a1]. The last occurrence of e[q] is entered as e(q). 2. What do you know about the domains of e[q] and t? 3. What do you know about the other parameters in the problem? 4. Why do you believe this function has a maximum? The reason for the last question can be seen in the following limit:
pi[a1]:=Q*log(t)*(e[q]) - ( P*(1/log(t))*(1/(1- (e[q])))) - log(y)*log(t)*(e[q]*s[q]+(1-e[q])*s[p]):
limit( pi[a1], e[q]=1, left ) assuming t::RealRange(Open(0),Open(1));
                       signum(P) infinity
How I got to this limit is another story. Basically, I noticed that the dependence on t was only through log(t). So, I made a substitution u=log(t). Because log is an increasing function, a critical point wrt (u,e[q]) is a critical point wrt (t,e[q]). While Maple still had trouble finding a critical point in terms of (u,e[q]), I was able to solve these equations step-by-step. What I found is that the critical point occurs when e[q]=1 and u^2=infinity. The worksheet I used to support my analysis can be downloaded from either of the following URLs: View 178_Max-Meade.mw on MapleNet or Download 178_Max-Meade.mw

View file details 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/~meade/
<\pre>
One way to achieve your desired result is to use ` ` to make make names out of the pieces of the expression. This preserves the structure but prevents Maple from thinking the arguments are numbers. For example: codegen[MathML](`5`^`4`); 5 4 5 4 `5`^`4` I hope this helps, 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/~meade/
Oops! The second argument should indicate the line of reflection. Here, I used:
plottools[reflect]( P1, [[0,0],[1,1]] );
<\pre>

You are correct that d and s have to be defined. You can use almost anything. Just to get a picture, try:

d := p->sin(p);
s := p->cos(p);
<\pre>

With apologies for my oversight.

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.ed
Jp, The simplest way to interchange the axes in a previously constructed plot is to use the reflect command in the plottools package:
P1 := plot( {d(p),s(p)}, p=0..50 ):
plottools[reflect]( P1 );
Note that the reflect command appears to suppress any axis labels that might have been created in the original plot. So, the plots[display] command would be needed to reinsert these labels:
with( plottools );
with( plots );
P1 := plot( {d(p),s(p)}, p=0..50, labels=["p","d,s"] ):
P1; # original plot, with labels
P2 := reflect( P1 ):
P2; # reflected plot, no labels
P3 := display( P2, labels=["d,s","p"] ):
P3; # reflected plot, with labels
Another approach would be to use a parametric representation of the curve:
P4 := plot( [ [d(p),p,p=0..50],
              [s(p),p,p=0..50] ], labels=["d,s","p"] ):
P4;
I hope these are 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
<\pre>
Here is how you can see the command used to prepare the plot: 1. start Plot Builder 2. add functions to be plotted 3. click Done 4. select type of plot to be produced 5. click Options (not Plot!) 6. select specific options for plot 7. click Command (not Plot!) 8. select and copy the plot command 9. paste this command at an input prompt 10. add a semicolon to the end of the command 11. execute the command 12. see your plot That's it. It's not as hard as it might seem. Steps 5, 7, and 9 are the critical ones. You should be warned that the plot command produced by this method is not always the simplest way to obtain the plot that you have created. I hope this helps, 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.ed
I think you will be happier with this approach. 1. Use a colon (:) at the end of the end if. This suppresses all output from the if .. then .. else .. end if regardless of whether you use a colon or semi-colon on the individual statements within the if block. 2. Put each result that you want to see in a print statement. For example:
with( plots ):
i:=1:
if (i=1) then
  p1:=plot(x,x=-1..1):
  p2:=plot(x^2,x=-1..1):
  print( display([p1,p2]) );
else
  ...
end if:
Note that I use square brackets instead of parentheses inside the display command. It makes no difference here, but it would if I were trying to set options for specific plots. (Sets are unordered, lists are ordered.) I hope this helps, 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
<\pre>

Dave, The command you seek is nops. If L is a list, then nops(L) gives the number of elements in the list. (In general, nops gives the number of operands of an expression.) This can be found on the help page for op or nops. I do agree that this could have been explicitly explained on the help page for lists. I hope this helps, 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 problem is that RootOf is used for symbolic/analytic roots. Since you have complex-valued floating-point values for the coefficients, this is problematic. Either solve or fsolve will return the full set of (approximate) solutions. You can then extract the roots that you want/need or work through the full set of 4 solutions: nz1 := solve(E1*NZ^4 - NZ^2*(2*E1^2) + E1^3 - E1*E2^2=0); -2.374580431 + 2.675986179 I, 2.374580431 - 2.675986179 I, -2.325160311 + 2.623756975 I, 2.325160311 - 2.623756975 I nz1 := fsolve(E1*NZ^4 - NZ^2*(2*E1^2) + E1^3 - E1*E2^2=0); -2.374580431 + 2.675986179 I, -2.325160311 + 2.623756975 I, 2.325160311 - 2.623756975 I, 2.374580431 - 2.675986179 I 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
The basic problem, as I see it is that Maple's limit command implements a continuous limit. Just a couple weeks ago I talked about this limitation with people from Maplesoft at the Maple Conference. My recommendation is to implement an extension to the limit command that would fill this void. The syntax I recommend is:
       limit( f(n), n=infinity, discrete );
<\pre>

While this is of no immediate use to you, I can also show you a way to get the results that you seek.

Your approach is close, but does not force Maple to make the assumptions early enough. What you need to do is form each subsequence, then (expand and) simplify under the assumption BEFORE giving the result to limit. Try:

restart;

X := n->(-1)^n*(2-1/n)*sin(n*Pi/3);

X0 := simplify(expand( X(3*n  ) )) assuming n::posint;
X1 := simplify(expand( X(3*n+1) )) assuming n::posint;
X2 := simplify(expand( X(3*n+2) )) assuming n::posint;

seq( X0, n=1..5 );
seq( X1, n=1..5 );
seq( X2, n=1..5 );

limit( X0, n=infinity );
limit( X1, n=infinity );
limit( X2, n=infinity );
<\pre>

The expand is needed to get Maple to expand the sine functions; simplify is needed to remove instances of (-1)^even and (-1)^odd. The goal is to reduce the expression to one that coicides with a subsequence formed from a function defined for all real numbers (so that Maple's continuous limit command can find the correct result).

I hope some of 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
This is pretty straight-forward except that Maple does not do well if the full integral is entered all in one command. The problem is that the bounds on a definite integral are not automatically applied as ``assumptions'' when the integrand is being evaluated. To see what this means, try evaluating J1 and J2 below without the assuming phrase.
restart;
qx := piecewise(x>20 and x<40,1/21,0);
I1 := Int( qx(x1),x1=x-10..40 );           
I2 := Int( qx(x1),x1=x-20..40 );                    
J1 := value( I1 ) assuming x>20 and x<30;
J2 := value( I2 ) assuming x>20 and x<30;
I3 := Int( J1*J2, x=20..30 );
J3 := value( I3 );
                             26500 / 1323 
------------------------------------------------------------- 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 not sure why you have a problem with this. I just tried your example and it works for me. I see that the extension (.exe) is not needed, but does not cause any problems either. I assume you can launch notepad from outside the procedure. 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 am not aware of a specific Maple command for finding the dominant eigenvalue of a matrix. But, this is not too difficult to implement. The approach I will take is to sort the eigenvalues according to their complex modulus. Then, the dominant eigenvalue (if it exists) will be the first element of the resulting list. To begin, find the eigenvalues of a floating-point matrix: > restart; > with( LinearAlgebra ): > A := convert( RandomMatrix( 5, 5 ), float ): > Lambda := Eigenvalues( A ): Note now that the modulus of the vector of eigenvalues can be found with the abs command: > abs(Lambda); To re-order the eigenvalues by decreasing modulus, use the sort command with a second argument that provides the order function: > Lambda1 := abs(l)>abs(m) )[]>; Note that the first argument of sort needs to be a vector. Thus, it is necessary to convert the vector of eigenvalues into a list. In the end, the list is coverted back to a vector. All of this can be put together in a single procedure as follows: > EigenvalueSort := A -> > (a,b)->abs(a)>abs(b) > )[]>:

This procedure can be used, as follows:


> EigenvalueSort( RandomMatrix(9,9) );



-------------------------------------------------------------------
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
First 40 41 42 43 44 Page 42 of 44