Carl Love

Carl Love

23646 Reputation

25 Badges

9 years, 266 days
Natick, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity


These are answers submitted by Carl Love

There are a few reasons why it doesn't work, and several ways to fix it. I'm guessing that you intend for x1x2x3 to be the roots of eq? But you never substituted them in for x in eq. I think this is the simplest fix that does what I think you really want:

restart:
eq := x-> a*x^3 + b*x^2 + c*x + d;
 s := x1 + x2 + x3 = -b/a;
 sp := x1*x2 + x1*x3 + x2*x3 = c/a;
 p := x1*x2*x3 = -d/a;
sol:= ()-> solve({eq~([x1,x2,x3])[], p, s, sp}, {x1, x2, x3},explicit);
sol(); #no solution, yet                       
 a := 1;
b := -5;
c := 6;
d := -1;
 sol(); #very long solution quickly returned   

 

It seems to be a bug such that the columns of the selector matrix A are used in reverse order. You can use this instead:

eval~(`if`~(A, B, 'NULL'));

If I ignore the typeset formulas in your Question---which necessarily have errors as alreadyb pointed out---and I focus instead on the text of your Question, then I think that the command BesselJZeros is the major part of the solution. See help page ?BesselZeros.

The number of terms isn't a mathematically invariant property of expressions. In other words, it's not a property that's necessarily preserved by transformations under which expressions remain mathematically equal. Rather, it's a structural property. So, expand shouldn't be used. This is sufficient:

nTerms:= (e::algebraic)-> `if`(e::`+`, nops(e), 1)

Here's a procedure that returns the values of the four formulas in your Question. This code likely requires 1D Input. I can modify it if need be, but I strongly recommend learning 1D Input.

ANOVA_3D:= (X::And(rtable, 3 &under rtable_num_dims))->
local
    N:= numelems(X), 
    Add:= ArrayTools:-AddAlongDimension,
    Xij:= Add(X, 3),
    SS:= A-> add(A^~2)*numelems(A)/N,
    SSi, SSj, all    
;
    Record(
        "SS__P"=  (SSi:= SS(Add(Xij, 2))) - (all:= SS(<add(Xij)>)),
        "SS__A"=  (SSj:= SS(Add(Xij, 1))) - all,
        "TSS"=    SS(X) - all,
        "SS__AP"= SS(Xij) - SSi - SSj + all
    )
:
#Example usage:
(n, k, r):= (90, 100, 110): #There'll be nearly a million entries.
#essentially the same as Array, but allows the random initializer:
X:= rtable((1..n, 1..k, 1..r), frandom(0..2, 1), datatype= hfloat):
R:= ANOVA_3D(X);
                                                               
R := Record(SS__P = 35.6539282805752, SS__A = 35.7316001050640, 
  TSS = 330279.334091332, SS__AP = 2965.63736086327)

#Access values like this:
R:-SS__P;
                        35.6539282805752

Let me know when you want to access values from the F-distribution.

You'll need 3 initial conditions. I used 3 random values in interval 0..1. This system of ODEs is known to be particularly sensitive to changes in initial conditions; the term "butterfly effect" was essentially invented because of this system.

For the plot, I used thickness= 0.1. A very fine line is ideal for this. But you'll need Maple 2022 for that. If you have earlier Maple, change 0.1 to 0. Maple has numerous schemes to color the plot.

restart:
Digits:= 15:
V:= [x,y,z]:
eqs:= diff~(V(t), t)=~ [35*(y-x), -x*z-7*x+28*y, x*y-3*z](t);
sol:= dsolve({eqs[], (V(0)=~ V||~0)[]}, numeric, maxfun= -1, parameters= V||~0);
sol(parameters= ['rand(0.0..1.)()' $ 3]);
plots:-odeplot(sol, V(t), t= 0..50, numpoints= 5000, thickness= 0.1);

To get the first entry from each pair (where "first" is with respect to the order displayed, which is not necessarily the order that they were entered) do

op~(1, ({entries}@op~@ListTools:-Classify)(lhs, AA));

This should work in Maple 12:

c:= plots:-spacecurve([cos(t), sin(t), t], t= -Pi..Pi);
map2(op, 1, indets(c, specfunc(anything, CURVES)))[1];

 

Like this:

restart:
oldGAMMA:= eval(GAMMA):
unprotect(GAMMA):
GAMMA:= overload([
    proc(x::And(integer, Not(0)), y::And(algebraic, Not(complexcons)), $)
    option overload;
        'procname'(x,y)
    end proc,
    oldGAMMA
]):
protect(GAMMA, oldGAMMA)
:
#Test:
GAMMA(-1,x), oldGAMMA(-1,x), GAMMA(3);

 

Like this (I've divided all the quantities by 1000 to reduce visual clutter):
 

Network Flow

 

There are 3 regions in this network: a, b and c.

Region a contains 2 nodes, b has 4, and c has 5

 

restart:

(a,b,c):= (2,4,5):

Let "x[i,j] be a quantity of interest that moves from region a to b, where 1<=i<=a, and 1<=j<=b."

Similarly, "y[j,k]  is a quantity from region b to c,  where 1<=j<=b, and 1<=k<=c."

"z[i,k]  is a quantity from region a to c, where 1<=i<=a, and 1<=k<=c."

 

I wish to determine the flow required from from regions a to c, a to b and b to c, that produce the minimum cost to each of the 5 nodes in region c.

The x-array represents a quantity from  i to  j; the y-array corresponds to  j to k, and the z-array concerns i to k

X:= Matrix((a,b), symbol= x):
Y:= Matrix((b,c), symbol= y):
Z:= Matrix((a,c), symbol= z):

Each node in region c requires a quantity

RegionC:= <5, 15, 8, 10, 15>:

The capacities at region a 

RegionA:= <90, 75>:

The capacities at region b

RegionB:= <35, 20, 30, 15>:

Cost from each node in region a to each node in region b

Cost1:= <
    2, 1, 3/2,   3;
  5/2, 2, 7/2, 3/2
>:

Cost from each node in region b to each node in region c

Cost2:= <
    3/2, 4/5, 1/2, 3/2,   3;
      1, 1/2, 1/2,   1, 1/2;
      1, 3/2,   2,   2, 1/2;
    5/2, 3/2, 3/5, 3/2, 1/2
>:

Cost from nodes in region a to c

Cost3:= <
    11/4, 7/2, 5/2, 3,   5/2;
       3, 7/2, 7/2, 5/2, 2
>:

The objective function

Cost__Total:= (add@(add@`*`~)~)([Cost||(1..3)], [X,Y,Z]):

 

I wish to impose 4 constraints

CapB:= add(X[i], i= 1..a) <=~ RegionB:

CapA:= add(<X|Z>[..,j], j= 1..b+c) <=~ RegionA:

ReqC:= add(<Y,Z>[i], i= 1..a+b) >=~ RegionC:

InEqOutB:= add(<X, -Y^%T>[i], i= 1..a+c) =~ 0:

Cons:= seq~({CapA, CapB, ReqC, InEqOutB}):

The flow solution

Sol:= Optimization:-LPSolve(Cost__Total, Cons, assume= nonnegative):

 

The solution matrices are

(X,Y,Z):= (round~)~(eval([X,Y,Z], Sol[2]))[]:

(AtoB, BtoC, AtoC):= DataFrame~([X,Y,Z])[];

"AtoB,BtoC,AtoC:=[[[,1,2,3,4],[1,0,20,8,0],[2,0,0,0,15]]],[[[,1,2,3,4,5],[1,0,0,0,0,0],[2,0,15,5,0,0],[3,5,0,0,0,3],[4,0,0,3,0,12]]],[[[,1,2,3,4,5],[1,0,0,0,0,0],[2,0,0,0,10,0]]]"

MinCost:= round(Sol[1]*1000);

103800

N:= GraphTheory:-Graph({
    seq(seq(`if`(X[i,j]=0, NULL, [[A__||i, B__||j], X[i,j]]), i= 1..a), j= 1..b),
    seq(seq(`if`(Y[j,k]=0, NULL, [[B__||j, C__||k], Y[j,k]]), j= 1..b), k= 1..c),
    seq(seq(`if`(Z[i,k]=0, NULL, [[A__||i, C__||k], Z[i,k]]), i= 1..a), k= 1..c)
});

GRAPHLN(directed, weighted, [A__1, A__2, B__2, B__3, B__4, C__1, C__2, C__3, C__4, C__5], Array(1..10, {(1) = {3, 4}, (2) = {5, 9}, (3) = {7, 8}, (4) = {6, 10}, (5) = {8, 10}, (6) = {}, (7) = {}, (8) = {}, (9) = {}, (10) = {}}), `GRAPHLN/table/1`, Matrix(10, 10, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 20, (1, 4) = 8, (1, 5) = 0, (1, 6) = 0, (1, 7) = 0, (1, 8) = 0, (1, 9) = 0, (1, 10) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 0, (2, 4) = 0, (2, 5) = 15, (2, 6) = 0, (2, 7) = 0, (2, 8) = 0, (2, 9) = 10, (2, 10) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 0, (3, 4) = 0, (3, 5) = 0, (3, 6) = 0, (3, 7) = 15, (3, 8) = 5, (3, 9) = 0, (3, 10) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = 0, (4, 5) = 0, (4, 6) = 5, (4, 7) = 0, (4, 8) = 0, (4, 9) = 0, (4, 10) = 3, (5, 1) = 0, (5, 2) = 0, (5, 3) = 0, (5, 4) = 0, (5, 5) = 0, (5, 6) = 0, (5, 7) = 0, (5, 8) = 3, (5, 9) = 0, (5, 10) = 12, (6, 1) = 0, (6, 2) = 0, (6, 3) = 0, (6, 4) = 0, (6, 5) = 0, (6, 6) = 0, (6, 7) = 0, (6, 8) = 0, (6, 9) = 0, (6, 10) = 0, (7, 1) = 0, (7, 2) = 0, (7, 3) = 0, (7, 4) = 0, (7, 5) = 0, (7, 6) = 0, (7, 7) = 0, (7, 8) = 0, (7, 9) = 0, (7, 10) = 0, (8, 1) = 0, (8, 2) = 0, (8, 3) = 0, (8, 4) = 0, (8, 5) = 0, (8, 6) = 0, (8, 7) = 0, (8, 8) = 0, (8, 9) = 0, (8, 10) = 0, (9, 1) = 0, (9, 2) = 0, (9, 3) = 0, (9, 4) = 0, (9, 5) = 0, (9, 6) = 0, (9, 7) = 0, (9, 8) = 0, (9, 9) = 0, (9, 10) = 0, (10, 1) = 0, (10, 2) = 0, (10, 3) = 0, (10, 4) = 0, (10, 5) = 0, (10, 6) = 0, (10, 7) = 0, (10, 8) = 0, (10, 9) = 0, (10, 10) = 0}, order = C_order))

GraphTheory:-DrawNetwork(N);

 


 

Download LPTransNetw.mw

@Earl (Perhaps that title should be "A get-Maple-to-do-more-algebra approach.") I should've shown that simplification of W. I think that that would've avoided that version difference of implicitplot. And it'll make your attempt to solve W for y work better. Indeed, by solving for x and separately and plotting both, you can get a perfect plot (shown below).

We simplify W by using the fact that it's a rational function of x and y implicitly equated to 0. Thus,

  1. Its denominator doesn't matter and can be discarded. This is the purpose of numer below.
  2. The "content" (common factors of all numerator terms) can be discarded. This is the purpose of op&<2 in 
    op~&<1 @ op&<2 @ evala@AFactors below.  (My &< is defined below.) That @-expression is equivalent to
    e-> op~(1, evala(AFactors(e))[2]).
  3. Exponents on factors of the numerator can be discarded. This is the purpose of 
    op~&<1
restart
:
#abbreviations for things that I commonly use in @-expressions:
`&<`:= curry: `&>`:= rcurry:   
:
WattEq:= r^2 - (b^2 - (a*sin(theta) + sqrt(c^2 - a^2*cos(theta)^2))^2);


#Don't use decimal values for the parameters.
params:= (`=`~ &< [a,b,c])~([[31,11,30]/10, [1,sqrt(2),1], [21,22,6]/10]);


rng:= -b..b:
V:= [x,y]:
abc:= e-> eval(e, params[Params]):
 
Params:= 1: #which parameter values to use

W:= (op~&<1 @ op&<2 @ evala@AFactors @ numer @ evala@Norm @ eval[recurse])(
    abc(WattEq),
    #polar to cartesian: 
    [r= sqrt(x^2+y^2), (cos,sin)(theta)=~ V[]/~r]
):
if nops(W)<>1 then WARNING("multiple factors") fi:
W:= mul(W);


plots:-implicitplot(abc([W, V[]=~ rng])[], scaling= constrained);

For the purely functional/algebraic plot, we continue with the same W, and solve it for x and y separately. This'll lead to much duplication of the plotted points, but that doesn't matter; the gaps left by each will be filled by the other. The solutions for x are plotted in parametric form (e.g., plot([f(y), y, y= a..b])) so that the axes are switched. The solutions are substituted into a procedure template that does hardware-float evaluation. An attempt is made to distinguish these two cases:

  1. If a value has an imaginary part of very small magnitude, it's assumed that that is due solely to rounding error and Re(e) is returned;
  2. If there's an imaginary part of larger magnitude, it's assumed that the value is truly non-real and undefined is returned.

The plot command doesn't care if some points are included that have undefined as one of their coordinates. It simply skips them. However, if a function has an undefined coordinate for all of its points, plot will issue a warning. To avoid these warnings, I remove those functions that always produce non-real values. There are 12 functions to start (as can be seen from the degree of W in and y), and this filtering removes 4 of them.

So, here's the code that does all that: 

Digits:= 15: eps:= 1e-13: Ev:= evalhf: #All 3 may need to be adjusted.

#some testing values to decide real-or-complex issue; 99 is arbitrary:
T:= abc(b)/99 *~ {$1..99}:
 
(Px,Py):= (remove &< (w-> w~(T)={undefined}))~(
    (((w-> subs~&<(_P=~ w)) @ subs&<(V=~ z) @ [solve])~(W, V))( 
        z-> local e:= Ev(_P); `if`(abs(Im(e)) < eps, Re(e), undefined)
    )
)[]:
plot([`[]`~(Px, z->z, abc(rng))[], Py[]], abc(rng), scaling= constrained);

the

The above takes 0.516 seconds on my computer, and, because using evalhf makes the filtering and plotting nearly instantaneous, the vast majority of that time is spent in solve.

There's a chance (version-dependent I guess) that the solve will return RootOf expressions. In that case, add the explicit option. The polynomial is bicubic in both x and y, so it definitely can be explicitly solved.

There's a chance that some spurious horizontal lines will appear on the plot. These can be corrected by adjusting Digitseps, or both. I haven't investigated the origins of these spurious lines. There's a remote chance, I guess, that you'll need to change evalhf to evalf to get higher precision, although I haven't yet seen a case where that's actually needed.

First, here are some things that you're doing wrong in your code snippet:

  1. By row, I assume that you mean Row from LinearAlgebra. Note that its name is capitalized.
  2. When the left operand of := is a function call (such as your row(out, k)), it is not evaluated (i.e., the function is not called) before the assignment[*1]. Rather, these assignments are stored in something called the remember table of row. (At the moment, it's not useful to explain remember tables further; but I'll be happy to when it becomes actually useful.)
  3. To do an assignment to a left operand that needs to be evaluated first, use assign, as in assign(f(x)= y).
  4. I assume that out is a matrix, but it would need to be already created, and have the correct dimensions for your idea to have any chance of working. (And there are several other reasons why this won't work, but they're too complicated to explain right now. Anyway, there are several other easy ways to do what you want.)

Second, just for the sake of showing you some briefer syntax, here's a slight variation on Joe's Answer:

#shortcut "hand-typed" matrix constructor syntax:
M:= <1, 2, "A"; 2, 3, "B"; 3, 4, "A">;

M[[seq](`if`(M[i,3]="A", i, NULL), i= [rtable_dims](M)[1])];

Third, there's a matrix-based superstructure called a DataFrame that's perfect for fancy row-selection operations such as I suspect you have in mind. A DataFrame is just an Object whose underlying structure is an ordinary Matrix with many OOP-style fancy indexing methods overloaded onto the regular square-bracket indexing operator.

In the code below, is the matrix from above, and columns is a list of names for the columns. There's nothing special about the names that I chose ([1, 2, key]).

DM:= DataFrame(M, columns= [1, 2, key]):
DM[DM[key]=~ "A"];


Footnote [*1]There's only one exception to this that I'm aware of: When function cat is the left operand of :=, then it IS evaluated before the assignment is made, as in

cat(x, 3):= 4;

which assigns 4 to x3.

 

Like this:

W:= `[]`~(P, Q);

An operator is called a binary infix operator if it can take 2 arguments and can be placed between (rather than in front of) them. So, most familar arithmetic operators (such as +, -, *, /, =, <, ^, and others) are binary infix operators. The way that you wrote your Question suggests that you want to use oper in infix form. This is easy in Maple, and it's one of my favorite features of Maple's syntax. The operator's name should begin with &, followed by letters or a wide variety of punctuation marks.

AnyOp:= proc(`&op`, x, y) x &op y end proc:
AnyOp(`*`, 7, 3), AnyOp(`+`, 7, 3); 

                             
21, 10

When the operator is not being used in infix form, it must be in back quotes (notice `&op` immediately after proc). This is true for all infix operators, not just those beginning with (notice `*` and `+` above). On the other hand, when they are used in infix form, they must NOT have quotes.

In the examples above, the operation that is passed to `&op` need not be infix as long as it can accept two arguments (see iquo and irem below):

AnyOp:= (`&?`, x, y)-> x &? y:
AnyOp~([`*`, `+`, `-`, `/`, `mod`, iquo, irem, `^`], 7, 3);

                  [21, 10, 4, 7/3, 1, 2, 1, 343]

Since mod is one of Maple's infix operators, even though it's just alphabetic characters, it needs the quotes when it's passed.

 

The plot can be made perfectly (no gaps, no jaggedness) and easily (no fancy options or extra points or precision needed) by these steps:

  1. Express your Watt equation as an expression implicitly equated to 0.
  2. Convert the expression from polar coordinates to Cartesian (x,y). (There's a command for that, but rather than looking up its syntax, I just used the usual well-known substitutions.)
  3. Rationalize the resulting algebraic function with evala@Norm. (The resulting rational function can be simplified significantly, which might be useful if you're going to do further algebra on it; but, it plots with no trouble without simplification.)
  4. Put in your numeric parameters.
  5. Plot with implicitplot.
restart:
WattEq:= r^2 - (b^2 - (a*sin(theta) + sqrt(c^2 - a^2*cos(theta)^2))^2);
params:= [a= 3.1, b= 1.1, c= 3.]:  rng:= -1.2..1.2:
W:= (evala @ Norm @ eval[recurse])(
    WattEq,
    #polar to cartesian: 
    [r= sqrt(x^2+y^2), (cos,sin)(theta)=~ (x,y)/~r]
);
  

plots:-implicitplot(eval(W, params), (x,y)=~ rng, scaling= constrained);

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