Carl Love

Carl Love

27286 Reputation

25 Badges

11 years, 359 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity


These are answers submitted by Carl Love

The command Worksheet:-WorksheetToMapleText should do most of the work of converting a worksheet to a procedure. You'll still need to do a little clean-up by hand. Save the worksheet without output (Ctrl-D) before doing this.

Yes, the not-publically-documented command for doing this is `tools/genglobal`. Most commands that generate arbitrary constant names use it. I believe that dsolve used to use it but doesn't anymore, instead always using _C1_C2, ..., possibly aliased to c__1c__2, .... So, we just put a little wrapper around commands to change the constants:

restart:
SeqConst:= proc(e, prefix::symbol:= :-_C, new_prefix::symbol:= :-_C)
    subsindets['flat', 2](e, suffixed(prefix, nonnegint), new_prefix, `tools/genglobal`[1]) 
end proc
:
SeqConst(dsolve(diff(y(x),x)=2*y(x)));

                             y(x) = _C2 exp(2 x)

SeqConst(dsolve(diff(y(x),x)=3*y(x)));
                             y(x) = _C3 exp(3 x)
 

Here is my Maple implementation of the Mathematica command Cases that you showed. If that command allows more-complicated patterns or wildcards that you'd like implemented, let me know. The objects in the pattern are either types or the triple-underscore zero-or-more-of-anything wildcard. (All explicit numbers are also considered types in Maple.)

WildcardPatternPredicate:= proc(Pat::list({type, identical(:-___)}))
local P:= Array(1..0, datatype= integer[1]), t, k, np:= nops(Pat);
    for k,t in Pat do
        P,= `if`(k=1, "proc(", `if`(k<np or t::type, ", ", ""));
        if t::type then P,= sprintf("%a::(%a)", _||k, t);
        elif k=np then break
        elif (t:= Pat[k+1]) = :-___ then next
        else P,= sprintf("%a::seq(Not(%a))", _||k, t)
        fi
    od;
    if t::type then P,= ", $" else np-- fi;
    P,= sprintf(") option overload; %a; true end proc", _||np);
    overload([parse(String(P)), ()-> false])@op
end proc
:
Cases:= (L, P::list({type, identical(:-___)}))->
    select(WildcardPatternPredicate(P), L)
:
Cases(combinat:-permute(4), [___,2,___,3,___]);
   [[1, 2, 3, 4], [1, 2, 4, 3], [1, 4, 2, 3], [2, 1, 3, 4], 
     [2, 1, 4, 3], [2, 3, 1, 4], [2, 3, 4, 1], [2, 4, 1, 3], 
     [2, 4, 3, 1], [4, 1, 2, 3], [4, 2, 1, 3], [4, 2, 3, 1]]

#Large-case time test:
CodeTools:-Usage(Cases(combinat:-permute(9), [___,2,___,3,___])):
nops(%) = 9!/2;
memory used=1.99GiB, alloc change=169.75MiB, 
cpu time=3.16s, real time=6.99s, gc time=2.55s

                        181440 = 181440


The brilliant idea of using overload with the seq parameter modifier to effectively allow type-checking of expression sequences is due to @sursumCorda . Turning the patterns into predicates would be much more complicated without this idea.

This method is somewhat slower than most of the other methods presented in this thread, but it's not horribly slow. However, the patterns that it can handle are much more complicated than just permutations; in particular, the non-wildcards can be any type; they needn't be specific entities such as 2 or 3

I don't know the extent of all your cases, but in this case the following does it:

primpart(factor(expr));
                        40 x - 24 y - 3

primpart returns what's left after all the constants have been factored out.

Your error was mixing up the argument order for has. The object being searched is its first argument; the things being searched for (disjunctively) are its second argument.

I know that you like rule-based and pattern-based constructs. Unfortunately, the command applyrule (and its ilk, such as patmatch) has very limited capability, was written decades ago, the writer likely left Maple decades ago, and likely will never be significantly improved. So, you might as well stop torturing yourself by trying to use these commands, because it's usually impossible. (There might be a way to use applyrule in this case, but I'm not going to take the time to find out.)

Here's a way with subsindets:

subsindets(
    eq, `*`,   
    proc(e) 
    local (S,R):= selectremove(type, e, anything^(And(fraction, 2 &under denom)));
        `if`(denom(R)=1 or S=1, e, R %* S)
    end proc
);    

There are a vast number of features of the Maple language that do not work in 2D Input. Using `?` in symbols (names) is one of them. For many of the features that don't work in 2D, that fact is mentioned on the help page of the feature.

I've always found the depends parameter modifier to be undependable. (Note for other readers: This has nothing to do with the depends command.) I'd guess that your goal is to make the 2 essentially a parameter of a procedure-valued procedure so that such procedures can be constructed on-the-fly for pattern matching. Here's a way to do it:

f3:= subs(_ty= 2, proc(x::seq(Not(_ty)), y::_ty) [[x], [y]] end proc):
f3(0, 1, 2);

                         [[0, 1], [2]]

The repeated use of abs(...)^2 in the equation strongly suggests that a complex solution is expected. To get complex solutions from fsolve, include complex as the last argument. Doing that, I get

             tg := -3.401652391 - 3.401652391*I

There may be other complex solutions.

Factoring the non-constant terms of the polynomial is sufficient to control the wild oscillations:

restart:
p0:= expand(2*sin(45*arcsin(x/2)));
p:= factor(p0) - sqrt(7-sqrt(5)-sqrt(30-6*sqrt(5)))/2;
plot(p, x= -2.1..2.1, numpoints= 1000, size= [1500,300]);

The stated task can be done in Maple by

{seq}([seq](p[]), p= Iterator:-TopologicalSorts(4, {2<3}));

{[1, 2, 3, 4], [1, 2, 4, 3], [1, 4, 2, 3], [2, 1, 3, 4], [2, 1, 4, 3], [2, 3, 1, 4],
 [2, 3, 4, 1], [2, 4, 1, 3], [2, 4, 3, 1], [4, 1, 2, 3], [4, 2, 1, 3], [4, 2, 3, 1]}

You can put any number of inequalities in the second argument.

This command does not generate all the permutations and then filter them to get the ones that match the pattern. Rather, it generates only the ones that match in the first place.

You wrote:

L:=combinat:-permute([1,2,3,4]);
map(X->`if`(has(X,2) and has(X,3) and ListTools:-Search(2,X)<ListTools:-Search(3,X),X,NULL) ,L);

The commands has and member are not the same thing! Using has for a job that can be done by member can be very inefficient and sometimes incorrect. Also, one shouldn't use map for a job that can be done by select (however, that's only slightly inefficient). Your code can be replaced by

select(P-> member(2,P,'p') and member(3,P,'q') and p<q, combinat:-permute(4));


Here's a little example that highlights the crucial difference between has and member:

L:= [3, x+2];
has(L,2);
                              true
member(2,L)
                             false

 

All your for loops could be replaced by a single very short elementwise operation:

Test:= Q =~ Multi

See help page ?elementwise.

It looks like you're trying to invert a symbolic 5x5 matrix by solving a system of equations. You could instead simply ask Maple for the inverse:

MAinv:= 1/MA

There is a tremendous amount known about symbolic matrix inverses, but perhaps you want to see what you discover on your own. That's fine. I won't spoil it by disclosing any findings. But I'll hint to you that you'll discover the same things easier if you make MA a 3x3 matrix. And you can make its entries such that you can tell immediately which row and column they come from:

M:= <a1,a2,a3; b1,b2,b3; c1,c2,c3>

I've made the following syntax changes, and now it runs:

  1. Replaced all "dot" multiplication by multiplication
  2. Put after (1+F[s])
  3. Changed D[a] to D__a

I also clarified the input style (indentation, line breaks, removal of superfluous parentheses, etc.), mostly to help me see where syntax changes might be needed.

eqn1:= {
    (1+1/beta)*diff(f(eta), eta, eta, eta) 
    - (1+F[S])*(diff(f(eta), eta)^2)
    + diff(f(eta), eta, eta)*f(eta) 
    + M*(A-diff(f(eta), eta))
    - 1/(R*D__a)*diff(f(eta), eta)
    = 0,
 
    diff(theta(eta), eta, eta)
    + 4/3*N*diff((1+(K-1)*theta(eta))^3*diff(theta(eta), eta), eta)
    + Pr*f(eta)*diff(theta(eta), eta)
    + (1+1/beta)*Pr*Ec*diff(f(eta), eta, eta)^2
    + M*Ec*(diff(f(eta), eta)-A)^2 
    = 0,

    #BCs at eta=0:
    f(0) = 0, theta(0) = 1+alpha*D(theta)(0), 
    D(f)(0) = 1+(1+1/beta)*lambda*(D@@2)(f)(0), 

    #BCs at eta=10:
    D(f)(10) = 0, theta(10) = 0
}: 
#Supply numeric values for parameters:
sys1:= eval(
    eqn1, {
        A = 1, Ec = .2, K = 2.5, M = .5, N = .5, Pr = 6.5, R = 1,
        alpha = .5, beta = 2, lambda = .5, D__a = .3, F[S] = 1
    }
):

sol1:= dsolve(sys1, numeric);

 

I was skeptical about whether there was any theoretical justification for the square root term that I included in yesterday's Answer. Although it worked amazingly well for fitting your data, intuitive dimensional analysis tells me that sqrt(amps) is a bogus concept. I am not a battery expert, although I have long been curious about the thing that you've measured.

After lengthy web research (mostly on the technical pages of high-tech battery companies), I found that there is a model for precisely the phenomenon that you are measuring: Peukert's law, which says that t = C/i^k, where t is the time to drain the battery to a pre-set cut-off voltage at constant current i, and C and k are empirical constants specific to the battery. The exponent is called the Peukert constant for the battery, and it's a function of the battery chemistry and age. The range k = 1.05..1.2 is typical for lead-acid batteries for example. The is a function of the battery's size or "capacity"; its dimensions are (current x time).

This model fits your data extremely well (R-squared = 99.98%)! Assuming that you're just a battery hobbyist, the accuracy of your measurements is impressive.

Regarding plots combining data points and a function fitted to them: It's best to plot the data as discrete points. If they're plotted in connect-the-dots form, that tends to visually exaggerate the residuals.

Peukert's Law

A model for the time needed to discharge a battery to a specific voltage at constant current

 

Maple author: Carl Love <carl.j.love@gmail.com> 2024-September-6

Data collected by Thomas Dean

restart:

# Amps
A := <444, 218, 74, 312/10, 209/10, 119/10, 625/100>
# Time until battery discharges to 10.5 volts
T := <5, 15, 60, 3*60, 5*60, 10*60, 20*60>:

Model_Peukert:= unapply(Statistics:-PowerFit(A,T,i,summarize), i);

Summary:
----------------
Model: 14182.414/i^1.28408026766874
----------------
Coefficients:
              Estimate  Std. Error  t-value  P(>|t|)
Parameter 1    9.5598    0.0975      98.0834   0.0000
Parameter 2   -1.2841    0.0240     -53.4568   0.0000
----------------
R-squared: 0.9998, Adjusted R-squared: 0.9997

proc (i) options operator, arrow; HFloat(14182.413917653834)/i^HFloat(1.2840802676687426) end proc

plot(
    [<A|T>, Model_Peukert], (min..max)(A), style= [point, line],
    color= [COLOR(RGB, 1, 0, 0), COLOR(RGB, 0, 0.7, 0)],
    symbol= diagonalcross, symbolsize= 20, thickness= 2,
    view= [-max(A)/20..max(A), -max(T)/20..max(T)],
    labels= [
        `#mn("current")`*``(Unit(ampere)),
        `#mn("time")`*``(Unit(minute))
    ],
    axes= frame,
    labeldirections= [horizontal, vertical], labelfont= [helvetica, 16],
    legend= ["measured     ", "Peukert"],
    title= cat(
        "Time to discharge to 10.5 V at constant current\n"
        "Peukert model:  ",
        t = subsindets(Model_Peukert(i), float, evalf[3]),
        "\n"
    ),        
    titlefont= [times, 16], gridlines= false
);

 

Download PeukertsLaw.mw

Like this:

FlattenSet:= (s::set)-> subsindets(s, set, (x-> `if`(x::set, x[], x))~)

The crucial difference between this and @vv 's Answer is that this only unpackages sets that are members of other sets.

The minimum shown is, of course, a bug. The correct value can be obtained by using the location option:

Min:= minimize(G(x,y),x=0..1,y=0..1, location);
                        9  
   Min := 5.838093910 10 , {[{x = 0.8949243201, y = 0.9715796788}, -0.7102080449]}

eval(G(x,y), Min[2][1][1]);
                          -0.710208045

Now, the erroneous value is still shown as the first number, but the true minimum is shown as the last number. Its correctness can be confirmed by direct evaluation and your plot.

1 2 3 4 5 6 7 Last Page 2 of 390