Carl Love

Carl Love

27599 Reputation

25 Badges

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

MaplePrimes Activity


These are answers submitted by Carl Love

This code does essentially the same thing as @dharr 's but does it with the index-free and functional syntax that I prefer. The entire computation---the construction of the Array---is done by a single line of code. The remaining two lines are simply to put that Array into a neatly displayed Table.

restart
:

Given: Two  leagues (lists of teams). It isn't necessary that the teams or the leagues have the same number of members, althoughh that

happens to be true in this case.

new:= [
    ["Pablo G", "Lucy Z", "Matt C"], ["John S", "Rosalie B", "Morgan B"],
    ["Daniel N", "Cohen M", "Steve S"], ["Thomas S", "Peter M", "Tony W"],
    ["Susanne A", "Stephen J", "Jorah S"], ["Sapan B", "Rishab S", "Jesse B"],
    ["Dean B", "Jesse W", "Hak Y"]
]:
past:= [
    ["Jesse W", "Daniel N", "Lucy Z"], ["John S", "Jesse B", "Rosalie B"],
    ["Jorah S", "Peter M", "Thomas S"], ["Stephen J", "Pablo G", "Matt C"],
    ["Rishab S", "Cohen M", "Steve S"], ["Tony W", "Dean B", "Morgan B"],
    ["Susanne A", "Hak Y", "Sapan B"]
]:

Goal: Find all intersections of new teams with past teams.

#Array of the intersections of the pairs of the Cartesian product of the 2 leagues.
#The raw Array is more suitable than the displayed Table for programmatic access:
new_Xinter_past:= ArrayTools:-GeneralOuterProduct({op}~(new), op@`intersect`, {op}~(past));

"[[["Lucy Z",,,"Matt C","Pablo G",,,],[,"John S","Rosalie B",,,,"Morgan B",],["Daniel N",,,,"Cohen M","Steve S",,],[,,"Peter M","Thomas S",,,"Tony W",],[,,"Jorah S","Stephen J",,,"Susanne A"],[,"Jesse B",,,"Rishab S",,"Sapan B"],["Jesse W",,,,,"Dean B","Hak Y"]]]"

#Display as a Table. The middle operator (()-> ...) of the @-construct is
#for the elegant display of Table entries, without extraneous punctuation:
(Tabulate @ (()-> if nargs=0 then else [args, if nargs=1 then "" else fi] fi)~ @ DataFrame)(
    new_Xinter_past, rows= ["new "||($nops(new))], columns= ["past "||($nops(past))]
):

Download PairwiseIntersection.mw

The Table didn't transfer to MaplePrimes, but you can see it in the worksheet. It's identical to @dharr's both in content and displayed form.

You ( @Christopher2222 ) are a Moderator, and thus you can read the deleted content. If you do this, you will see that the vast majority of Moderator-deleted (as opposed to author-deleted) content is spam.

You are wrong about lost knowledge: It is the failure to delete inappropriate content that will cause the greater loss of knowledge, because good writers will stop using MaplePrimes, and it will become irrelevant. Do you think periodicals shouldn't have editors and peer review?

The defined coordinate systems are stored in an unexported table of an  undocumented package, as given in the title. Here's a procedure to access it:

Coords:= proc(C::{name,function}, V::list(name))
uses PC= Plot:-CoordinateSystems;  #undocumented package
local
    opq:= kernelopts(opaquemodules= false),    
    r:= PC:-CoordinateSystemsTable[String(`if`(C::name, C, op(0,C)))](
        V[],`if`(C::name, [][], op(C))
    )
;
    kernelopts(opaquemodules= opq);
    r
end proc
:
Coords(spherical_physics, [r, theta, phi]);
      [r*cos(phi)*sin(theta), r*sin(theta)*sin(phi), r*cos(theta)]

Since the package is undocumented, don't rely too much on the syntax remaining exactly the same as the way I just accessed it in the procedure.

By using a Newton's method iteration with pure-integer builtin (i.e., GMP) arithmetic, I can find roots of integers significantly faster than with any method that uses evalf. While most of Maple's integer-arithmetic commands whose names begin with i are builtin, that isn't true of iroot: It's written in Maple, and it uses evalf.

My Newton procedure is the second procedure listed in the table methods in the worksheet below. The Newton iteration for finding the cube root r of is simply

r-> (2*r^3 + n)/(3*r^2)

restart
:

#Make garbage collection single threaded, because multi-threaded makes
#"real-time" comparisons difficult:
kernelopts(gcmaxthreads= 1)
:

#
# The procedures to be tested and compared
# ----------------------------------------
methods:= table([
    #iroot with floor correction:
    Iroot= (n-> local r:= iroot(n,3); `if`(r^3 > n, r-1, r)),

    #Integerized Newton's method:
    Newton= proc(n)
    local r1:= integermul2exp(1, trunc(ilog2(1+n)/3 + 1/2)), r, r2;
        do r1:= iquo(2*(r2:= (r:= r1)^2)*r + n, 3*r2) until r-r1 in {-1,0,1};
        min(r,r1)
    end proc,

    #Two of @dharr's methods:
    root@evalf= (n-> floor(root(evalf(n),3))),
    evalf@root= (n-> floor(evalf(root(n,3)))),

    #iroot's own method, minimalized:
    evalf@`^`= (n-> trunc(evalf(n^(1/3))))
]):
meth_names:= [indices](methods, nolist);

[`@`(evalf, `^`), `@`(root, evalf), Iroot, Newton, `@`(evalf, root)]

# Generate random test cases:
#
nWords:= 3:  iterations:= 2^16:
R:= rand(2^((nWords-1)*64)..2^(nWords*64)-1):
#Arbitrary key for consistent test data on different runs and platforms.
#Change the key to generate different random test data.
randomize(23):
testdata:= table(
    [meth_names[], accuracy]=~ ['['R'()$iterations]' $ 1+numelems(methods)]
):

#Set sufficient floating-point precision (see showstat(iroot), line 39)
#(only needed for methods that directly use evalf):
d:= length(2^(nWords*64)-1):  Digits:= iquo(d,3)+5+length(d);

26

# The tests:
#
#Accuracy test: If all methods agree, then all will be listed in a single set:
entries(ListTools:-Classify(j-> methods[j]~(testdata[accuracy]), meth_names), nolist);

{Iroot, Newton, `@`(evalf, `^`), `@`(evalf, root), `@`(root, evalf)}

#Efficiency test:
for j in meth_names do
    f:= methods[j];
    printf("%a:\n", j);
    CodeTools:-Usage(seq[reduce= ()](f(n), n= testdata[j]));
    printf("\n")
od:

evalf@`^`:
memory used=1.59GiB, alloc change=-16.00MiB, cpu time=4.27s, real time=7.89s, gc time=750.00ms

root@evalf:
memory used=1.79GiB, alloc change=0 bytes, cpu time=3.69s, real time=9.42s, gc time=578.12ms

Iroot:
memory used=485.51MiB, alloc change=0 bytes, cpu time=953.00ms, real time=1.96s, gc time=218.75ms

Newton:
memory used=211.56MiB, alloc change=0 bytes, cpu time=750.00ms, real time=1.42s, gc time=203.12ms

evalf@root:
memory used=3.62GiB, alloc change=0 bytes, cpu time=13.97s, real time=25.77s, gc time=1.73s
 

Newton's method is the winner.

 

The Newton-method procedure for any particular root be generated by this procedure:

NewtonIroot:= (k::And(posint, Not(1)))-> subs(
    _k= k,
    proc(n::nonnegint)
    local r1:= integermul2exp(1, trunc(ilog2(1+n)/_k + 1/2)), r, r2;
        do r1:= iquo((_k-1)*(r2:= (r:= r1)^(_k-1))*r + n, _k*r2)
        until r-r1 in {-1,0,1};
        min(r,r1)
    end proc
):

showstat((Iroot5:= NewtonIroot(5)));


proc(n::nonnegint)
local r1, r, r2;
   1   r1 := integermul2exp(1,trunc(1/5*ilog2(1+n)+1/2));
   2   do
   3       r1 := iquo(4*(r2 := (r := r1)^4)*r+n,5*r2)
       until r-r1 in {-1, 0, 1};
   4   min(r,r1)
end proc
 

r:= rand();

441235504072

n:= (r+1)^5-1:  Iroot5(n);

441235504072

 

Download Iroot.mw

The votes awarded after deletion do apply to the author's reputation. It has always been true (at least for the 10 years or so that I've been using MaplePrimes) that votes awarded to an item which is subsequently deleted continue to apply to the author's reputation.

I totally agree with acer that items Deleted As Spam should not be retained. There are now several Posts where the majority of Replies have been deleted as spam yet remain as clutter.

The syntax of Maple's hypergeom always has exactly three arguments. The first two arguments are arbitrary-length lists of algebraic expressions, and the third is a single algebraic expression rather than a list. If it's a 2F1 function, then the first list has 2 elements, and the second has 1.

The biggest and ugliest (IMO) differences between Maple and Mathematica syntax are that Maple's lists are enclosed in square brackets ([...]) rather than curly braces ({...}) and that Maple's function arguments are always in parentheses ((...)) rather than square brackets. (I consider the Mathematica to be uglier.)

In another Answer, @janhardo shows the correct syntax of a simple Maple hypergeom function (of class 2F1): 

hypergeom([1,1], [2], -z)

Your example could be entered as

-r^2/3*((1-beta)/`γ`)^n*
    hypergeom(
        [n, n/(1-beta)], [(beta-1-n)/(beta-1)], -r^(3/n*(1-beta))*chi^2/`γ`
    )

I used `γ` instead of gamma because gamma is by default a pre-defined constant (similar to Pi). If you mean the constant, then just use gamma. If it's intended to be a variable or parameter, there are several workarounds, one of which is `γ`.

My original Answer was wrong because I made M^*.M a symmetric matrix when actually it's hermitian. I thank dharr for pointing out my error. Here's a corrected Answer. The incorrect one is at the bottom of this Answer for the sake of continuity of reading.

I will use for your because the default sans-serif font on MaplePrimes essentially makes l unreadable. The trick to doing your problem is expressing L as Lx + I*Ly where Lx and Ly are assumed real, and simplifying conjugated expressions with simplify(..., conjugate) and evalc.

In Maple, M^%H (or M^*) is the HermitianTranspose of M, and A.B is the matrix product of matrices A and B. Spelling out MatrixMatrixMultiply or HermitianTranspose is unnecessary.

Because I use evalc (which assumes variables are real), the assume(a::complex) is needed.

The final eigenvectors below can be further simplified, but it requires subtlety.

restart:

interface(showassumed= 0):

assume(Lx::real, Ly::real, a::complex, k>0, Lx^2+Ly^2 < 1);

L:= Lx+I*Ly:

B:= simplify(k*(1-L*conjugate(L)));

-k*(Lx^2+Ly^2-1)

M:= <L-B, -a*B; 0, L+B>;

Matrix(2, 2, {(1, 1) = Lx+I*Ly+k*(Lx^2+Ly^2-1), (1, 2) = a*k*(Lx^2+Ly^2-1), (2, 1) = 0, (2, 2) = Lx+I*Ly-k*(Lx^2+Ly^2-1)})

MS:= simplify(evalc~(M^%H.M), conjugate);

Matrix(2, 2, {(1, 1) = (Lx^2+Ly^2-1)^2*k^2+2*Lx*(Lx^2+Ly^2-1)*k+Lx^2+Ly^2, (1, 2) = a*k*(Lx^2+Ly^2-1)*(Lx-I*Ly+Lx^2*k+Ly^2*k-k), (2, 1) = k*(Lx^2+Ly^2-1)*conjugate(a)*(Lx+I*Ly+k*(Lx^2+Ly^2-1)), (2, 2) = (Lx^2+Ly^2-1)^2*k^2*abs(a)^2+Ly^2+(Lx-k*(Lx^2+Ly^2-1))^2})

EV:= LinearAlgebra:-Eigenvectors(MS, output= vectors);

Matrix(2, 2, {(1, 1) = -2*a*((-Lx^2-Ly^2+1)*k+I*Ly-Lx)/(-sqrt((Lx^2+Ly^2-1)^2*k^2*abs(a)^4+(4*(Lx^2+Ly^2-1)^2*k^2+4*Lx^2+4*Ly^2)*abs(a)^2+16*Lx^2)+k*(Lx^2+Ly^2-1)*abs(a)^2-4*Lx), (1, 2) = -2*a*((-Lx^2-Ly^2+1)*k+I*Ly-Lx)/(sqrt((Lx^2+Ly^2-1)^2*k^2*abs(a)^4+(4*(Lx^2+Ly^2-1)^2*k^2+4*Lx^2+4*Ly^2)*abs(a)^2+16*Lx^2)+k*(Lx^2+Ly^2-1)*abs(a)^2-4*Lx), (2, 1) = 1, (2, 2) = 1})

eval(map(simplify, EV, {L=`&ell;`, Lx-I*Ly= CL}), CL= conjugate(`&ell;`));

Matrix(2, 2, {(1, 1) = -2*a*(`&ell;`*conjugate(`&ell;`)*k+conjugate(`&ell;`)-k)/(k*(-`&ell;`*conjugate(`&ell;`)+1)*abs(a)^2+sqrt(k^2*(`&ell;`*conjugate(`&ell;`)-1)^2*abs(a)^4+(4*k^2*(`&ell;`*conjugate(`&ell;`)-1)^2+4*`&ell;`*conjugate(`&ell;`))*abs(a)^2+4*(`&ell;`+conjugate(`&ell;`))^2)+2*conjugate(`&ell;`)+2*`&ell;`), (1, 2) = 2*a*(`&ell;`*conjugate(`&ell;`)*k+conjugate(`&ell;`)-k)/(k*(`&ell;`*conjugate(`&ell;`)-1)*abs(a)^2+sqrt(k^2*(`&ell;`*conjugate(`&ell;`)-1)^2*abs(a)^4+(4*k^2*(`&ell;`*conjugate(`&ell;`)-1)^2+4*`&ell;`*conjugate(`&ell;`))*abs(a)^2+4*(`&ell;`+conjugate(`&ell;`))^2)-2*conjugate(`&ell;`)-2*`&ell;`), (2, 1) = 1, (2, 2) = 1})

 

Download SymbEigenvectors.mw

 

Original incorrect worksheet:

restart:

L:= Lx+I*Ly:

B:= simplify(k*(1-L*conjugate(L))) assuming real;

-k*(Lx^2+Ly^2-1)

M:= <L-B, -a*B; 0, L+B>;

Matrix(2, 2, {(1, 1) = Lx+I*Ly+k*(Lx^2+Ly^2-1), (1, 2) = a*k*(Lx^2+Ly^2-1), (2, 1) = 0, (2, 2) = Lx+I*Ly-k*(Lx^2+Ly^2-1)})

MS:= Matrix(evalc~(M^%H.M), shape= symmetric) assuming a::complex;

Matrix(2, 2, {(1, 1) = (Lx+k*(Lx^2+Ly^2-1))^2+Ly^2, (1, 2) = a*((Lx+k*(Lx^2+Ly^2-1))*k*(Lx^2+Ly^2-1)-I*Ly*k*(Lx^2+Ly^2-1)), (2, 1) = a*((Lx+k*(Lx^2+Ly^2-1))*k*(Lx^2+Ly^2-1)-I*Ly*k*(Lx^2+Ly^2-1)), (2, 2) = conjugate(a*k*(Lx^2+Ly^2-1))*a*k*(Lx^2+Ly^2-1)+(Lx-k*(Lx^2+Ly^2-1))^2+Ly^2})

EV:= LinearAlgebra:-Eigenvectors(MS, output= vectors) assuming (Lx,Ly)::~real, abs(L)^2 < 1, k>0;

Matrix(2, 2, {(1, 1) = -2*a*((-Lx^2-Ly^2+1)*k+I*Ly-Lx)/(-sqrt(k^2*(Lx^2+Ly^2-1)^2*abs(a)^4-8*k*Lx*(Lx^2+Ly^2-1)*abs(a)^2+4*a^2*k^2*(Lx^2+Ly^2-1)^2-8*a^2*(I*Ly-Lx)*(Lx^2+Ly^2-1)*k+(4*a^2+16)*Lx^2-(8*I)*a^2*Lx*Ly-4*Ly^2*a^2)+k*(Lx^2+Ly^2-1)*abs(a)^2-4*Lx), (1, 2) = -2*a*((-Lx^2-Ly^2+1)*k+I*Ly-Lx)/(sqrt(k^2*(Lx^2+Ly^2-1)^2*abs(a)^4-8*k*Lx*(Lx^2+Ly^2-1)*abs(a)^2+4*a^2*k^2*(Lx^2+Ly^2-1)^2-8*a^2*(I*Ly-Lx)*(Lx^2+Ly^2-1)*k+(4*a^2+16)*Lx^2-(8*I)*a^2*Lx*Ly-4*Ly^2*a^2)+k*(Lx^2+Ly^2-1)*abs(a)^2-4*Lx), (2, 1) = 1, (2, 2) = 1})

eval(map(simplify, EV, {L=`&ell;`, Lx-I*Ly= CL}), CL= conjugate(`&ell;`));

Matrix(2, 2, {(1, 1) = -2*a*(k*`&ell;`*conjugate(`&ell;`)+conjugate(`&ell;`)-k)/(k*(-`&ell;`*conjugate(`&ell;`)+1)*abs(a)^2+sqrt(k^2*(`&ell;`*conjugate(`&ell;`)-1)^2*abs(a)^4-4*k*(`&ell;`+conjugate(`&ell;`))*(`&ell;`*conjugate(`&ell;`)-1)*abs(a)^2+4*a^2*(`&ell;`*conjugate(`&ell;`)-1)^2*k^2+8*a^2*conjugate(`&ell;`)*(`&ell;`*conjugate(`&ell;`)-1)*k+4*`&ell;`^2+8*`&ell;`*conjugate(`&ell;`)+4*conjugate(`&ell;`)^2*(a^2+1))+2*conjugate(`&ell;`)+2*`&ell;`), (1, 2) = 2*a*(k*`&ell;`*conjugate(`&ell;`)+conjugate(`&ell;`)-k)/(k*(`&ell;`*conjugate(`&ell;`)-1)*abs(a)^2+sqrt(k^2*(`&ell;`*conjugate(`&ell;`)-1)^2*abs(a)^4-4*k*(`&ell;`+conjugate(`&ell;`))*(`&ell;`*conjugate(`&ell;`)-1)*abs(a)^2+4*a^2*(`&ell;`*conjugate(`&ell;`)-1)^2*k^2+8*a^2*conjugate(`&ell;`)*(`&ell;`*conjugate(`&ell;`)-1)*k+4*`&ell;`^2+8*`&ell;`*conjugate(`&ell;`)+4*conjugate(`&ell;`)^2*(a^2+1))-2*conjugate(`&ell;`)-2*`&ell;`), (2, 1) = 1, (2, 2) = 1})

``

Below, I have put your data in exactly the same order as you gave it; however, I've grouped it with square brackets and commas. Please note where the brackets are; the pattern should be obvious. Line breaks, spacing, and indentation are irrelevant; they're only used to enhance readability.

Pts:= [
    [
        [80, 1, 1.80074340195803263], [80, 2, 1.79160585795803287], [80, 3, 1.7805446209580329], 
        [80, 4, 1.76755968995803269], [80, 5, 1.75265106495803280], [80, 6, 1.73581874695803272],
        [80, 7, 1.71706273595803272]
    ],
    [                        
        [85, 1, 1.71097991748770582], [85, 2, 1.70184237348770577], [85, 3, 1.69078113648770581], 
        [85, 4, 1.67779620548770559], [85, 5, 1.66288758048770571], [85, 6, 1.64605526248770562],
        [85, 7, 1.62729925148770591]
    ],                          
    [
       [95, 1, 1.53550477309879818], [95, 2, 1.52636722910871089], [95, 3, 1.51530599205435690], 
       [95, 4, 1.50232106105435669], [95, 5, 1.48741243605435680], [95, 6, 1.47058011805435672], 
       [95, 7, 1.45182410705435672]
    ],
    [                    
       [100, 1, 1.44826461660029512], [100, 2, 1.43912707260029623], [100, 3, 1.42806583560029569],               [100, 4, 1.41508090460029576], [100, 5, 1.40017227960029588], [100, 6, 1.38333996160029579],       
       [100, 7, 1.36458395060029579]
    ]
]: 
plots:-display(PLOT3D(MESH(Pts)), labels= [x,y,z]);

Note that PLOT3D and MESH are all uppercase letters. Some help for them can be found at ?plot,structure; however, the best way to learn about them is to deconstruct the results of plot3d (lowercase) commands with op.

An equivalent command, producing exactly the same plot, is

plots:-surfdata(Pts, labels= [x,y,z])

Either way, the trickier part is gettiing the data arranged correctly in a list of lists of lists (as shown above) or a 3D-Array (as shown by acer). Once that's done, displaying it is easy.

 

The trick is to insert a fake symbolic exponent and substitute 1 for it after the extraction, like this:

Powers:= (e, T::type)-> local One;
    eval(
        indets(combine(evalindets(e, T, `^`, One), power), T^anything),
        One= 1
    )
:
e:=_C1^2+_C1+_C1*_C2^3+a+b:
Powers(e, suffixed(_C, nonnegint));
                       /        2     3\ 
                      { _C1, _C1 , _C2  }
                       \               / 

By the way, suffixed(_C, ...is a subtype of symbol, so your And(symbol, ...) is redundant.

It's not a caching issue. What's happening is that the elapsed time between your calls to randomize is much shorter than the "ticks" of the clock. There's never a need to call randomize more than once. Indeed, the effect of calling randomize twice in a short time period is to reset the random-number generator to a specific entry in its sequence, so that's about as unrandomized as you can get.

My feeling about randomize is that it should only be called at the top level, never in a procedure.

Also, there's never a need to pass randomize a custom clock value. The simple call randomize() (with no arguments) already uses a clock-based value that's sufficient for any purpose.

Like this:

collect(T, indets(T, function));
                   
f(x) + 3 + (b + 1) g(x)

@Kitonum is correct that shading regions of polar plots is difficult in Maple, and a "polygon" approach is needed. Here's a way to do it that's perhaps more intuitive than his: I extract the polygon vertices directly from the output of polarplot commands. There's no need to look at the numerical vertices; it's all done rather automatically. Some additional benefits of my approach are that I put the plots on Maple's polar axes rather than Cartesian axes, and no explicit coordinate conversions are required.

I also show three distinct ways to get the area.

restart:

interface(prompt= "")
:

(r1, r2):= (-6*cos(theta),  2-2*cos(theta))
:

#Find theta values at the intersections:
theta__1:= rhs(solve({r1=r2, theta >= 0, theta <= Pi})[]);

(2/3)*Pi

I will extract the point-data matrices (rtables) for two of the relevant arcs for two purposes:
    1. To make filled polygon plots;
    2. To get the area of those polygons via the Shoelace formula (simply to verify integration accuracy). The formula is given below, but you don't
        need to understand it to understand the rest of this worksheet.

 

The point data is stored in Cartesian (rather than polar) coordinates; however, you don't need to know that to understand this worksheet.

(Circle__arc, Cardioid__arc):= (op@indets)~(
    plots:-polarplot~([r1, r2], theta=~ [Pi/2..theta__1, theta__1..Pi]),
    rtable
)[]:
#Reconnect stray end to origin to close the polygon:
Cardioid__arc:= <Cardioid__arc, <0|0>>:

The whole plot:

plots:-display(
    plots:-polarplot~(
        [r1, r2], theta=~ [Pi/2..Pi, 0..Pi], color=~ [red, blue],
        thickness= 2,
        legend=~ [r1, r2], coordinateview= [default, 0..Pi]
    ),
    plots:-polygonplot~(
        [Circle__arc, Cardioid__arc], transparency= .5,
        color=~ [pink, blue], legend=~ ["Area 1", "Area 2"]
    )
);

We want the area of the pink region + the area of the blue region, times 2. We can get that directly like this:

2*(Int(r1^2/2, theta= Pi/2..theta__1) + Int(r2^2/2, theta= theta__1..Pi)):
% = value(%);

2*(Int(18*cos(theta)^2, theta = (1/2)*Pi .. (2/3)*Pi))+2*(Int((1/2)*(2-2*cos(theta))^2, theta = (2/3)*Pi .. Pi)) = 5*Pi

An indirect way to get the area, just as easy, is the area of the circle minus the area the part of the circle outside the cardioid:

2*(Int(r1^2/2, theta= Pi/2..Pi) - Int(r1^2/2-r2^2/2, theta= theta__1..Pi)):
% = value(%);

2*(Int(18*cos(theta)^2, theta = (1/2)*Pi .. Pi))-2*(Int(18*cos(theta)^2-(1/2)*(2-2*cos(theta))^2, theta = (2/3)*Pi .. Pi)) = 5*Pi

And a way to approximate the area is the "Shoelace formula", which gives the area of any polygon simply from the Cartesian coordinates of its vertices. I do this because its easy since we already have the polygons.

Shoelace:= (XY::Matrix)->
    add(XY[..,1]*~XY[[2..,1],2]-XY[..,2]*~XY[[2.., 1],1])/2
:
2*Shoelace(<Circle__arc, Cardioid__arc>) = evalf(5*Pi);

HFloat(15.707812810646793) = 15.70796327

 

Download PolarArea.mw

Use sort with this key:

sort(F, key= [f-> subs(0= infinity, degree~(f, [x4,x3,x2,x1])), nops])

I assume that by "degree of nonlinearity in x4" you mean nothing more than simply "degree in x4 (regardless of the presence of other variables)". If that's not what you mean, let me know. Also, this solution assumes that the polynomials are expanded.

Change
yB_:= y -> sol4[1]:
to
yB_:= MakeFunction(sol4[1], y):

You are correct that the y in the solution is not being seen. This is because there is no y explicitly on the right side of the arrow. In order to make a non-explicit variable into a procedure parameter, you need to use MakeFunction instead of ->.

Creating the set of equations can be done as simply as this:

{seq}(Matrix(2, symbol= a) - Matrix(2, symbol= b) =~ Matrix(2, shape= identity));
       
{a[1, 1] - b[1, 1] = 1, a[1, 2] - b[1, 2] = 0, a[2, 1] - b[2, 1] = 0, a[2, 2] - b[2, 2] = 1}

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