Carl Love

Carl Love

28010 Reputation

25 Badges

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

MaplePrimes Activity


These are answers submitted by Carl Love

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.

This is the best that I could come up with so far. Its main drawback, to my mind, is that it creates a new DataFrame rather than just making a small change to the existing DataFrame. (In computer science lingo, we'd say that this doesn't operate inplace.)

data_2:= DataFrame(data_2, rows= [$1..upperbound(data_2, 1)])

My Answer is slightly different than @vv's because I assumed that you'd want to handle the more-general case of any number of rows in the new DataFrame. upperbound(..., 1) returns the number of rows, and the [$1.. ...] creates the list [1, 2, ..., nr], where nr is the number of rows.

I could easily have missed something easier. Although there are many help pages for DataFrame, they are not well organized IMO. In particular, the overview page ?DataFrame should at least list the name of every command specific to DataFrames, like the overview pages for most packages do.

Here's a way to find discrete local maxima. Suppose we have a list of (x,y)-points, sorted with respect to x, and we want the local maxima of the y-values.

#Generate test data by extracting a point-data matrix from a plot:
plot(sin(x), x= 0..5*Pi):
A:= indets(%, rtable)[]:

#Filter out consecutive equal y-values:
A1:= A[[seq](`if`(A[i,2]<>A[i+1,2], i, NULL), i= 1..upperbound(A)[1]-1)]:

#Find discrete local maxima:
convert(
    A1[[seq](
        `if`(A1[i,2] > max(A1[[i-1,i+1],2], i, NULL),
        i= 2..upperbound(A1)[1]-1
    )],
    listlist
);
        [[1.59043127800000, 1.], [7.87361657400000, 1.], [14.1568018700000, 1.]]

 

 

restart:

local D:

(A,B,C,D,S1,S2):= (
    [0,0,0], [1,0,0], [1,1,0], [0,1,0], [0,0,1], [0,0,sqrt(2)/2]
):

use Student:-MultivariateCalculus in
    A1:= Angle(Plane(S1,B,C),Plane(S1,C,D))*180/Pi;
    A2:= Angle(Plane(S2,B,C),Plane(S1,C,D))*180/Pi
end use;

60

180*arccos((1/6)*6^(1/2)*2^(1/2))/Pi

evalf(A2)*degrees;

54.73561030*degrees

 

Download PyramidAngle.mw

Just by blind guessing, I added a sqrt(x) term to the argument of the exponential, and I got a near-perfect fit! Both this new model and your original can be linearized, which gives more-accurate results and which can be probabilistically quantified.

MaplePrimes won't let me display my worksheet, but you can download it.

Download BatteryDrainage.mw

Here's the plaintext code from the worksheet:

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>:
f := x-> a*exp(1-x/b+c*sqrt(x));

#The above model f(x) can be linearized:
ln(f(x)) = (ln(a)+1) + (-1/b)*x + c*sqrt(x);
LinearParms:= <ln(a)+1, -1/b, c>:
P:= Statistics:-LinearFit(
    [1, x, sqrt(x)], A, ln~(T), x, output= parametervector, summarize
);

Summary:
----------------
Model: 8.2359404+.12972487e-1*x-.58328210*x^(1/2)
----------------
Coefficients:
              Estimate  Std. Error  t-value   P(>|t|)
Parameter 1    8.2359    0.2723      30.2450   0.0000
Parameter 2    0.0130    0.0027      4.7533    0.0089
Parameter 3   -0.5833    0.0647     -9.0137   0.0008
----------------
R-squared: 0.9920, Adjusted R-squared: 0.9880

Model:= eval(f(x), solve({seq}(P =~ LinearParms), {a,b,c}));
plot([<A|T>, Model], x= 0..444); 

 

5 6 7 8 9 10 11 Last Page 7 of 394