Carl Love

Carl Love

28085 Reputation

25 Badges

13 years, 99 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity


These are answers submitted by Carl Love

Assignments to sets and lists (whether by := or by subsop) are extremely inefficient. It is only begrudingly that Maple allows it for lists, and even that's only for small lists (upto size 100 I think). I think that that's allowed for the benefit of newbies doing homework problems (where efficiency doesn't matter much) who want to use lists as vectors or arrays. But nobody (in their right mind) wants to use a set as an array.

Assigning to set positions is even more inefficient than it is for lists because the new element needs to be put in the correct position (since sets are always stored sorted). So, even if S[k]:= x were allowed, after doing it, it wouldn't necessarily be true that S[k] = x or that nops(S) was the same as it was before. Those two facts alone make S[k]:= x too weird to consider using, regardless of efficiency concerns.

So, whenever you feel the urge to do S[k]:= x, you need to rethink that and come up with a better way to code it. There's always a better way. MaplePrimes can offer help in any situation like that.

If you change solve to eliminate, you'll get a less-trivial solution: a = (k^2-1)/3, b= 0, c= 0.

The Adams-Bashforth method is implemented in dsolve via

dsolve({...}, numeric, method= classical[adambash]);

See ?dsolve,numeric,classical

Use

int(1/(1+sqrt(x)), x= 0..1, method= ftoc);

To aid in remembering that, ftoc stands for fundamental theorem of calculus.

For the real-rooted polynomial, use the characteristic polynomial of a very sparse symmetric matrix all of whose entries are in {-1, 0, 1} of order 600. Then multiply that by any integer polynomial of degree 200 and expand. This only takes a few seconds (in modern Maple):

P:= expand(
   LinearAlgebra:-CharacteristicPolynomial(
      LinearAlgebra:-RandomMatrix(
         600, shape= symmetric, density= 1/200, generator= rand(-1..1)
      ), x
   )*
   randpoly(x, degree= 200)
);

 

Another way is to make the 1/2 a symbol, do the integration, then make the symbol equal 1/2. Thus:

restart:
f:= ((r-I*t)^(-s)-(r+I*t)^(-s))/(I/r):
int(f, t= 0..infinity) assuming s > 1, r > 0:
simplify(eval(%, r= 1/2));

 

On the second pass through the final loop (thus, when i=1), the dsolve doesn't find any solutions and it returns NULL. Thus S[1] is the NULL sequence. Thus S[1][1] has an "invalid subscript selector" because you're asking for the first member of an empty sequence.

I did manage to get a numeric solution from dsolve with no trouble.

If you want to apply the operation to column j of matrix A and put the results back into column j (thereby overwriting the original column j), you can do

A[..,j]:= evalf(convert~(A[..,j], temperature, Celsius, kelvin)):

Note that kelvin is with a lowercase k (because Kelvin means something else in Maple, completely unrelated to temperature).

In Maple, negative reals raised to fractional powers return the principal complex value. This is the best general answer for a system that generally handles complex numbers. But Maple also provides the command surd to handle the real case. So,

f:= x-> surd(4*x^2-4, 5)^4;

That means the 5th root to 4th power.

And here's another solution that I don't like as much, but I'll mention it because if I don't someone else likely will:

f:= proc(x) uses RealDomain; (4*x^2-4)^(4/5) end proc:

There are many ways to do it. I'd do

R:= unapply(index~(<P1, P2, P3, P4, P1 | Q1, Q2, Q3, Q4, Q1>, i), i):
plot([seq(R(i), i= 1..21)], color= plots:-setcolors()[1]);

That'll plot them all in the same color. If you'd like different colors, just remove the color option completely. The index command in the first line requires Maple version at least "Maple 2017". If you have an earlier version, I can very easily give you another command. Let me know.

 

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

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