Carl Love

Carl Love

28035 Reputation

25 Badges

12 years, 317 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity


These are answers submitted by Carl Love

For obvious reasons, the floating-point evaluation of frac(x) for large x is very prone to catasthropic cancellation[1]; indeed it's guaranteed to happen when abs(x) > 10^Digits and x is not an integer. Thus, it's amazing that this has not been corrected before. A correction is simple; here it is:

Frac:= (x::realcons)-> 
   `if`(hastype(x, float), frac(x), 'procname'(frac(x)))
:
`evalf/Frac`:= proc(x::realcons)
local oldDigits:= Digits, r;
   if x = 0 then return 0 fi;
   Digits:= Digits + 1 + max(ilog10~({1} union indets(x, rational)));
   r:= evalf(x);
   if ilog10(r) < -1 then
      while r=0 do
         Digits:= Digits + oldDigits;
         r:= evalf(x)
      od;
      Digits:= Digits + 2 - ilog10(r);
      r:= evalf(x)
   fi;
   evalf[oldDigits](r)
end proc:

Usage:

Frac(Pi^20);
evalf(%);

or
evalf(Frac(Pi^20));

I could overload frac so that it's diverted to Frac for the problematic arguments.

[1] Catasthropic cancellation is a phenomenon of floating-point evaluation (in all computer languages) which is a very common source of bugs. It happens when trying to compute the difference d:= x - y when the true value of abs(d) is much smaller than abs(x). In a base-b floating-point system with an m-digit mantissa, it's guaranteed to happen when 0 < abs(d/x) < b^(-m). (In the case of Maple, b=10 and m=Digits.)

The anomaly that you're experiencing is due to last name evaluation, the same issue that was at the root of the problem in your previous Question. It's not necessary for your procedure foo to return anything, but if you want to return something, it needs to be eval(r). You don't need to return anything because the statement r:-b:= 5 inside foo has already changed the global r. The same thing would happen with a Matrix.

Also, you're misusing the print command; you should be using eval instead.

It would be better if your procedure didn't return anything, because having it return something gives the false impression that it's acting on a local copy of the record when in fact it's acting on the global record. This remains true even if you use a different name for the local r.

Why do you have printlevel set to 0? The default value is 1, and that's what it should be if you want to see those plots. Also, why do you lower the values of prettyprint and verboseproc from their default values?

I think that it's not really considered a bug if you don't use enough digits of precision internally so that the overall computation returns with Digits precision. Maple's higher-level numeric evaluators for transcendental functions do internally adjust the value of Digits so that the returned result has the required precision. But frac is, I guess, too simplistic for that. Maybe I can easily correct that.

On the other hand, there is a library procedure `evalf/frac`, and that procedure does not adjust the value of Digits! There hardly seems to be any purpose for such a procedure if it doesn't accurately handle Digits. So, I think that that maybe should be considered a bug.

Assuming that your file contains the complete code, the result is identically 0. Indeed, your limit for every term is 0. Here is my analysis:

restart:

# I'm only interested in the integrand. I commented out
# the loop at the bottom of your code.

read "C:\\Users\\carlj\\Desktop\\MWE.txt":

type(integrand, `+`);

true

(1)

nops(integrand);

11360

(2)

indets(integrand, name);

{p, psi__p, t, t__1, t__2}

(3)

indets(integrand, function);

{_D175(t__2), _D677(t__1), _D678(t__1), _D680(t__1), _D681(t__1), cos(psi__p), exp(I*t), exp(I*2^(1/2)*t), exp(I*3^(1/2)*t), exp(p*(-cos(psi__p)^2)^(1/2)*t__1), exp(p*(-cos(psi__p)^2)^(1/2)*_D175(t__2)), exp(p*sin(psi__p)*t__1), exp(p*sin(psi__p)*_D175(t__2)), exp(-(1/2)*p*(-cos(psi__p)^2)^(1/2)*t__1), exp((1/2)*p*(-cos(psi__p)^2)^(1/2)*_D175(t__2)), exp(-(1/2)*p*sin(psi__p)*t__1), exp((1/2)*p*sin(psi__p)*_D175(t__2)), sin(psi__p)}

(4)

indets(integrand, function(dependent(t)));

{exp(I*t), exp(I*2^(1/2)*t), exp(I*3^(1/2)*t)}

(5)

indets(integrand, {dependent(t)^anything, anything^dependent(t)});

{1/(exp(I*t))^3, 1/exp(I*t), (exp(I*t))^3, 1/exp(I*2^(1/2)*t), 1/exp(I*3^(1/2)*t)}

(6)

I1:= eval(integrand, exp= 1):

C:= degree(I1, t);

0

(7)

#So, there are no terms containing any nonzero powers of t as factors.

No_T:= remove(depends, integrand, t):

evalb(No_T = 0);

true

(8)

#So, there are no constant terms.

#Putting it all together, **every** term is of form
#  A*exp(B*I*t)
#where A is some junk that doesn't depend on t and B
#is a manifestly real constant. The limit for any term
#like that is...

int(A*exp(B*I*t)/L, t= -L..L);

-I*A*(exp((2*I)*L*B)-1)*exp(-I*L*B)/(B*L)

(9)

limit(%, L= infinity);

0

(10)

 

Download MWEanalysis.mw

You said (in a Reply far down from the main Question) that you were having trouble implementing von Mangoldt's function. It is a bit tricky to make an efficient implementation. Here's mine:

vonMangoldt:= proc(n::posint)
local k, exact, p;
   if isprime(n) or n=1 then return ln(n) fi;
   k:= prevprime(ilog2(n)+1);
   do
      p:= iroot(n,k,'exact');
      if exact then return procname(p) fi;
      if k=2 then return 0 fi;
      k:= prevprime(k)
   od
end proc:

I leave it as an exercise for you to verify that it's correct. In particular

  1. Why is the maximum possible exponent ilog2(n)?
  2. Why do we only need to check prime exponents?
  3. Why is recursion occasionally used (see procname(p))?

Having read all of your code that you've posted here, I've seen that you must learn to use the i- forms of commands when they are available. They are much more efficient. In particular, learn isqrt, iroot, irem, iquo, and ilog2.

This is a common problem that often perplexes new users, including me when I was new. I agree that one should never copy-and-paste Maple output to subsequent input in a worksheet that one intends to save. It makes worksheets unreadable and difficult to modify. Fortunately, there's usually an easy way to address the issue with coding rather than copy-and-pasting.

Here's a basic way to do it:

eval([P2,Q2], {ans}[1]);
#The {} make it invariant between subsequent runs with the same Maple version.

or
eval([P2,Q2], {ans}[2]);

That's a lot better than copy-and-paste, but it's still not totally satisfying: How do you decide between 1 or 2? To answer that, I need to know what your selection criterion is. For example, do you want the solution with the smaller value of P2?

It's also possible to do this without any indexing, which I find preferable:

add(diff~(r^2, coords));

I'll answer both of your questions, but I'll answer the second one first. All that you need is

solve('identity'(Eq, x));

      {a = 4, b = -1, c = 2}

(I don't know if the identity option was available yet in Maple 13.)

So, you really don't need to extract the coefficients, which is why I answered the second question first. But, if you still want to extract the coefficients, it can be done by

C:= coeffs(lhs(Eq), indets(Eq, function));
solve({C});

which'll produce the same result as the solve(identity(...)).  For some more-general cases, the coefficient extraction may need to be expanded to:

T:= {x} union indets(Eq, function(dependent(x)));
C:= coeffs(collect((lhs-rhs)(Eq), T), T);

 

I reduced your M-animating program to this. I believe that this'll work in Maple 7. If not, let me know.

restart;
Col2:= ()-> 'color'= 'COLOR'('RGB', 'evalf[8](rand()/10.^12)' $ 3):

makeM:= proc(x0, y0, h, b, w)
local Hm;
   Hm:= h/3;
   plots[display](
      map(
         L-> plottools[polygon](map(`+`, L, [x0,y0]), Col2()),
         [   [[0,0], [0,h], [w,h], [w,0]],
             [[b/2,Hm], [w,.75*h], [w,h], [b/2,Hm+2*w*(h-Hm)/b]],
             [[b/2,Hm+2*w*(h-Hm)/b], [b-w,h], [b-w,.75*h], [b/2,Hm]],
             [[b-w,h], [b,h], [b,0], [b-w,0]]
         ]
      )
   ) 
end proc:

nframes:=47:
x0:=0: y0:=0: b:=7.8: w:=1.2:
plots[display](
   [seq(makeM(x0, y0, 12-0.37*i*(1-i/nframes), b, w), i= 1..nframes)],
   scaling= constrained, linestyle= 1, thickness= 0,
   insequence= true
);

 

Use

expr:= (-(exp(j*z))^2*Ei(1, j*z)+ln(-z)-ln(z)+Ei(1, -j*z))/(2*j*exp(j*z)):
evalc(expr) assuming j > 0, z > 0;

The resulting expression contains only Ei terms of a single real argument:

simplify(%);
             exp(j z) Ei(-j z) - Ei(j z) exp(-j z)
             -------------------------------------
                              2 j                 

 

You want the color to apply to the individual plots, not to the display as a whole. So the color option needs to go inside the parentheses of the individual plots, like this:

display(
     polygon([[x0,y0], [x0,y0+h], [x0+w,y0+h],[x0+w, y0]], color=red),
     polygon([[x0+b/2, y0+Hm],  [x0+w, y0+3*h/4],[x0+w,y0+h],[x0+b/2, y0+Hm+2*w*(h-Hm)/b]], color=white),
     polygon([[x0+b/2, y0+Hm+2*w*(h-Hm)/b],  [x0+b-w,y0+ h],[x0+b-w,y0+3*h/4],[x0+b/2, y0+Hm]], color=blue),
     polygon([ [x0+b-w,y0+ h],  [x0+b,y0+h], [x0+b,y0], [x0+b-w, y0]], color=turquoise),
     opts
):

And here's an unrelated problem that I spotted in your code: When you're making an animation each frame of which is made from several separate plots, then each frame should use display to bind it into a single plot. And then all the frames are put together into an animation with display(..., insequence= true). So, your code defining lM[i] should be:

lM[i]:= display(
     polygon([[x0,y0], [x0,y0+h], [x0+w,y0+h], [x0+w, y0]], color=red),
     polygon([[x0+b/2, y0+Hm], [x0+w, y0+3*h/4],[ x0+w,y0+h], [x0+b/2, y0+Hm+2*w*(h-Hm)/b]], color=white),
     polygon([[x0+b/2, y0+Hm+2*w*(h-Hm)/b],  [x0+b-w,y0+ h],[x0+b-w,y0+3*h/4],[x0+b/2, y0+Hm]], color=blue),
     polygon([[x0+b-w,y0+ h], [x0+b,y0+h], [x0+b,y0], [x0+b-w, y0]], color=turquoise),
     opts
):

So, there's no limit to how many levels at which you can use display.

I'm not sure if this example (below) shows the situation that you're talking about. Part of my confusion is that there are many meanings in Maple for what's colloquially called a "table". 

In the example, I make a column of random data called X to simulate the data that you've imported from Excel. Then I define a function F and apply it to each member of X to produce a new column that I display to the right of the original column, calling the resulting two-column table XY. Finally, I show one of many ways that a table can be displayed, a DataFrame.

If this example isn't what you're talking about, then you'll need to upload your worksheet to this site by using the green up-arrow on the toolbar of the editor for this site (MaplePrimes).
 

R:= rand(-99..99):  X:= Vector(8, R);

X := Vector(8, {(1) = -82, (2) = 80, (3) = -44, (4) = 71, (5) = -17, (6) = -75, (7) = -10, (8) = -7})

(1)

F:= x-> 2*x+3;

proc (x) options operator, arrow; 2*x+3 end proc

(2)

XY:= <X | F~(X)>;

XY := Matrix(8, 2, {(1, 1) = -82, (1, 2) = -161, (2, 1) = 80, (2, 2) = 163, (3, 1) = -44, (3, 2) = -85, (4, 1) = 71, (4, 2) = 145, (5, 1) = -17, (5, 2) = -31, (6, 1) = -75, (6, 2) = -147, (7, 1) = -10, (7, 2) = -17, (8, 1) = -7, (8, 2) = -11})

(3)

DataFrame(XY, columns= [x, 'f'(x)]);

Matrix([[``, x, f(x)], [1, -82, -161], [2, 80, 163], [3, -44, -85], [4, 71, 145], [5, -17, -31], [6, -75, -147], [7, -10, -17], [8, -7, -11]])

(4)

 


 

Download Columns.mw

A little bit of plotting shows that with the given parameter values, the first positive root is always between 0.7 and 0.9. So, I use this fsolve command (here embedded in a procedure):

restart:
F[lambda[i]]:= 
   (sin(lambda[i]) - tan(y*lambda[i])*cos(lambda[i]))*((1-1/B)*lambda[i]^2-k) + 
   (lambda[i]^3/B+lambda[i])*(cos(lambda[i]) - tan(y*lambda[i])*sin(lambda[i]))
: 
SolveF:= subs(
   _F= F[lambda[i]],
   proc(
      B::And(realcons, satisfies(B-> 1 <= evalf(B) and evalf(B) <= 30)),
      y::complexcons:= 0.9, 
      k::complexcons:= 0.003
   )
      fsolve(_F, lambda[i]= 0.7..0.9)
   end proc
):
SolveF(15);
plot(SolveF, 1..30);

The value of i is irrelevant, because i doesn't appear algebraically in the expression.

The above procedure is "quick and dirty": In particular, if you change the parameter values y and k from their defaults, the range 0.7 to 0.9 will likely need to be changed also.

I have two suggestions. I can't test these as I'm posting this from my phone at the moment. The first is to replace :- with :: in the stopat command. If that works, then that is the way to do it. If that doesn't work, set 

kernelopts(opaquemodules= false);

and then use the same stopat command that you originally tried. I'm almost certain that this will work, but the former option is better (if it works) than the latter because the latter causes a slight change to your testing environment that theoretically could cause different behavior than you'd have in your actual-use environment.

First 172 173 174 175 176 177 178 Last Page 174 of 395