Carl Love

Carl Love

24663 Reputation

25 Badges

10 years, 56 days
Natick, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity

These are answers submitted by Carl Love

You say that expr is always a sum of terms. Then this should do what you want:

for item in [op](expr) do ... od

It may come as a suprise that the complex-valued function (-2)^x is a very smooth, continuous function of its real argument x, but here it is:

plot([(Re,Im)((-2.)^x), x= -2..2], scaling= constrained);

The fact that the integral signs are gray rather than blue suggests that there's a possibility of evaluating them with the value command. Try value(%) after pdsolve.

It's possible to automate this process so that starting with any univariate procedure f we can obtain a parameter-shifted g without needing to explicitly compute the inverse of the shift (that inverse is found automatically, and only once). Like this:

f:= (n::algebraic)-> `if`(n::complexcons, `if`(n=0,1,0), 'procname'(n))
shift:= proc(F::uneval, N::name)
    f:= eval(F,1), e:= eval(op(f)), 
    n:= `if`(nargs=2, N, indets(e, And(name, Not(constant)))[])
    subs(_f= op(0,f), _s= {solve}(e= _n, n), _n= n, _n-> _f~(_s)[])
end proc
g:= shift(f(n+1/3)):
g(1/3), g(1), g(x);
                         1, 0, f(x - 1/3)


Change the first argument of the eval to

:-`.`(I__b, diff(`ω_`(t), t))

This restores (temporarily, for this statement only) the default matrix multiplication operator, which has been overloaded by Physics

An alternative is

use `.`= :-`.` in
exactly what you had ...)
end use;

My personal choice for this very common situation is return with no argument, as in

set_r:= proc(r::integer) thismodule:-r:= r; return end proc;

I usually do this even if the result of the previous command is always NULL.

I was about to give you the patfunc answer, but you beat me to it. So, here's a shortcut syntax:

evalindets[2](expr, 'patfunc(identical(R0), anything)', 1= r0, subsop);

Yes, that help page is on my "speed dial", so to speak. Surely I've read it a few hundred times. So let me know if you have any questions about it.

Did you see the Reply that I made to your previous Question earlier today?
I wonder if you had any comment about that.

Yes, it's suprisingly easy:

evala(Norm(x - (2^(1/2)+2^(1/3))));

Or, equivalently,

evala(Minpoly(2^(1/2)+2^(1/3), x));

Both of these commands will accept extra arguments allowing you to specify a larger field of coefficients. The default field is the rationals, of course.

Yes, you can define your own type-conversion procedures to be used with convert, like this:

`convert/coefflist`:= (s::series)-> local k, S:= [op](s), L:= S[2]; 
    ([seq](Array(L..S[-3], {seq}(S[k+1]=S[k], k= 1..nops(S)-3, 2))), subs(0= (), L))


convert(series(cos(x), x, 10), coefflist);
convert(series(x^2*sin(x), x, 10), coefflist);
convert(series(x^(-2)*sin(x), x, 10), coefflist);

In the 2nd and 3rd examples, the integer returned after the list is the lowest exponent, which is returned when it's other than 0.

The procedure makes use of the fact that array entries that aren't explicitly set are automatically 0.

Try this:

GainQ1:= 8*wx^2*(m-1)/(Pi^2*sqrt((m*wx^2-1)^2+Q1^2*wx^2*(wx^2-1)^2*(m-1)^2));
params:= [m= 6.3, Q1= 1];
a:= eval(GainQ1, params);
Top_of_Q1:= maximize(a, wx= 0.1..2., location);


In the procedure, replace print(eq) with 

entries(eq, 'nolist')

The print command is never an appropriate way for a procedure to return its results.

I don't think that Preben or Tom explicitly said this, so I guess that it's worth saying: In this case, assumptions on the parameters can't help. That's a mathematical fact, not a Maple limitation, because the Fundamental Theorem of Algebra tells you that there'll always be 4 solutions for w(x) that are constant with respect to x.

My general personal opinion about Maple assumptions is that they shouldn't be used (even when they're mathematically true) in cases where it's known that they can't help. This is because they occasionally prevent you from getting valid solutions, which you'll then never know about.

An easy way to select only the nonconstant solution is to include an initial condition that includes a constant of integration, like this:

dsolve({q, w(0)=C1})

Your command array should be Array. The lowercased array is a "deprecated" command, meaning that it's a long-obsolete command that shouldn't be used any more. Indeed, it was obsolete before the ExcelTools package even existed.

While that's certainly true, I don't think that it will fix your problem. It looks like Export automatically converts an array to an Array.

It looks like you're trying to learn Maple but that you have some experience with other computer languages. The code in my worksheet below may be a bit intense, but perhaps you're already aware of similar things in those other languages. Anyway, please look up what you can in Maple's help pages, and please ask any questions that you have about it. And don't spend too much time trying to understand a confusing help page; you can just ask about it here.

Given that your primary interest in your code is computing the maxes of Collatz sequences, these efficiency improvements are possible:

  1. (This one is easy, and obvious): You only need to check for a new max after doing a 3*i+1 step, not after an i/2 step.
  2. Memoization (this is subtle): Memoization refers to a procedure recording the results of its calls so that if called again with the same arguments, it doesn't need to redo the computation. Usually, this can be done very simply by putting option remember in the procedure's opening part. That doesn't work very well in this particular case, so my code uses an explicit table as a remember table.

Here's how the remember table works mathematically with the maxes of Collatz sequences (and this would work for any collection of finite sequences of reals sharing a common formula): Suppose that we have two positive integers a and b for which we want to compute the maxes Ma and Mb of their sequences A and B. First we compute A and Ma. We note the position in the sequence where Ma is attained. In the table, we record that Ma is the max for all sequences beginning with A[k] for k from 1 to K.

Next we compute and Mb. Suppose that at some point the sequence value in is one of  those A[k] for which we've already recorded a value. Then it's obvious that Mb is the larger of Ma and the largest value computed so far in B. So, we don't need to finish sequence B

Obviously that process (which I just detailed for 2 sequences) can be extended to any number of sequences. If we were only doing a small number of sequences, the extra work required to do this likely wouldn't be worth it. But it is worth it if you want the maxes of a large number of sequences. Below, I do it for all odd numbers upto 1 million in under 7 seconds. Without memoization, it takes over 60 seconds.


# This module requires 1-D input, aka plaintext or Maple Input.
# Some tools for working with Collatz sequences, especially for finding the maximum
# values in a large number of them.
CollatzMax:= module()
option `Author: Carl Love <>, 2022-Oct-5`;
uses LT= ListTools;
local A:= Array(1..1) #temp sequence storage; will be expanded as needed
    #remember table to map posints to the max of their Collatz seq:
    MaxCZ:= table('sparse'),
    Forget:= ()-> (MaxCZ:= table('sparse')), #table re-initializer

    #returns the max entry of the Collatz seq for a given posint:
    CollatzMax:= proc(n::posint)
    local i:= n, m, M:= n, k:= 1, K:= 0;
        A[k]:= n;
            M:= max((m:= MaxCZ[i]), i, M);
            if m<>0 then break fi; #Max found in table; use shortcut.
            if M=i then K:= k fi; #Current entry is max; record its position.  
            while i::even do i/= 2 od; #Compute next seq entry; repeat until odd.
            A(++k):= (i:= 3*i+1) #Record next entry.
        until i = 4;
        for k to K do MaxCZ[A[k]]:= M od;
    end proc,

    #applies CollatzMax to a seq of posints and returns the overall max:
    ModuleApply:= proc(b::posint) local a;
        max(seq(CollatzMax(a), a= 7..b, 2))
    end proc,

    #returns the set-valued functional inverse of table MaxCZ:
    InvertTable:= ()-> LT:-Classify(k-> MaxCZ[k], [indices](MaxCZ, 'nolist')),

    #returns just the Collatz seq of a posint, without any fancy stuff:
    Collatz:= proc(n::posint) local i:= n, k:= 1;
        A[k]:= n;
        while i<>1 do A(++k):= (i:= `if`(i::odd, 3*i+1, i/2)) od;
    end proc
end module

#The module's name used syntactically as a function calls the ModuleApply procedure:
M:= CodeTools:-Usage(CollatzMax(999999));

memory used=0.68GiB, alloc change=93.96MiB, cpu time=8.22s, real time=6.75s, gc time=2.16s


MaxProducers:= CodeTools:-Usage(CollatzMax:-InvertTable()):

memory used=0.69GiB, alloc change=1.89MiB, cpu time=6.02s, real time=3.70s, gc time=2.95s


{704511, 2113534, 3170302, 4755454, 7133182, 10699774, 16049662, 24074494, 36111742, 54167614, 81251422, 121877134, 182815702, 231376126, 274223554, 308501500, 347064190, 411335332, 520596286, 780894430, 1171341646, 1757012470, 2635518706, 2964958546, 3335578366, 3953278060, 4447437820, 5003367550, 7505051326, 11257576990, 16886365486, 25329548230, 37994322346, 56991483520}


[704511, 2113534, 1056767, 3170302, 1585151, 4755454, 2377727, 7133182, 3566591, 10699774, 5349887, 16049662, 8024831, 24074494, 12037247, 36111742, 18055871, 54167614, 27083807, 81251422, 40625711, 121877134, 60938567, 182815702, 91407851, 274223554, 137111777, 411335332, 205667666, 102833833, 308501500, 154250750, 77125375, 231376126, 115688063, 347064190, 173532095, 520596286, 260298143, 780894430, 390447215, 1171341646, 585670823, 1757012470, 878506235, 2635518706, 1317759353, 3953278060, 1976639030, 988319515, 2964958546, 1482479273, 4447437820, 2223718910, 1111859455, 3335578366, 1667789183, 5003367550, 2501683775, 7505051326, 3752525663, 11257576990, 5628788495, 16886365486, 8443182743, 25329548230, 12664774115, 37994322346, 18997161173, 56991483520, 28495741760, 14247870880, 7123935440, 3561967720, 1780983860, 890491930, 445245965, 1335737896, 667868948, 333934474, 166967237, 500901712, 250450856, 125225428, 62612714, 31306357, 93919072, 46959536, 23479768, 11739884, 5869942, 2934971, 8804914, 4402457, 13207372, 6603686, 3301843, 9905530, 4952765, 14858296, 7429148, 3714574, 1857287, 5571862, 2785931, 8357794, 4178897, 12536692, 6268346, 3134173, 9402520, 4701260, 2350630, 1175315, 3525946, 1762973, 5288920, 2644460, 1322230, 661115, 1983346, 991673, 2975020, 1487510, 743755, 2231266, 1115633, 3346900, 1673450, 836725, 2510176, 1255088, 627544, 313772, 156886, 78443, 235330, 117665, 352996, 176498, 88249, 264748, 132374, 66187, 198562, 99281, 297844, 148922, 74461, 223384, 111692, 55846, 27923, 83770, 41885, 125656, 62828, 31414, 15707, 47122, 23561, 70684, 35342, 17671, 53014, 26507, 79522, 39761, 119284, 59642, 29821, 89464, 44732, 22366, 11183, 33550, 16775, 50326, 25163, 75490, 37745, 113236, 56618, 28309, 84928, 42464, 21232, 10616, 5308, 2654, 1327, 3982, 1991, 5974, 2987, 8962, 4481, 13444, 6722, 3361, 10084, 5042, 2521, 7564, 3782, 1891, 5674, 2837, 8512, 4256, 2128, 1064, 532, 266, 133, 400, 200, 100, 50, 25, 76, 38, 19, 58, 29, 88, 44, 22, 11, 34, 17, 52, 26, 13, 40, 20, 10, 5, 16, 8, 4, 2, 1]

(J,E):= ([lhs~, rhs~]@{indices})(MaxProducers, pairs)[]:   

nops(J); #number of distinct maxes:


nops(`union`(E[])); #number of distinct starting values for which max is now known:


#It works just as well for huge multi-word integers:
n:= 2^(2^10)+rand():

memory used=2.00MiB, alloc change=0 bytes, cpu time=31.00ms, real time=20.00ms, gc time=0ns



memory used=2.14MiB, alloc change=0 bytes, cpu time=16.00ms, real time=15.00ms, gc time=0ns




In your 2nd attempt, you are missing a comma after stetsel in the solve command.

In your followup, you didn't load package LinearAlgebra​​​​​. 

2 3 4 5 6 7 8 Last Page 4 of 361