Carl Love

Carl Love

28040 Reputation

25 Badges

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

MaplePrimes Activity


These are replies submitted by Carl Love

@vv Thanks. Now I understand option threadsafe. It's just a way to flag a procedure as threadsafe in case the kernel incorrectly suspects otherwise; it overrides that suspicion.

Google has a much superior search algorithm, which is applicable here because Maple help pages are Google indexed (since they're somewhat duplicated as regular web pages). 

It burns me up that MaplePrimes Answers are not properly searchable by Google! The best hit that you'll get is the Answers page of the user who wrote the Answer. Since I have about 4000 Answers, someone Google searching for one of my Answers only knows that they have a positive hit; they'd need to wade through all the Answers to find it.

I got your email about this problem, where you attached a worksheet, which I re-attach here.

No Maple code that I have ever written (or ever will write) is intended to be used by Maple's 2D Input. Considering the length of the code, I refuse to look for the error in the 2D Input. If you redo it in the 1D Input form (code in a reddish monospaced font) that I originally posted it in, then I'll look at it. 


 

restart

Digits := trunc(evalhf(Digits));

15

(1)

ODEs := [diff(f(eta), `$`(eta, 3))+A^2+f(eta)*(diff(f(eta), `$`(eta, 2)))-(diff(f(eta), eta))^2+beta*((diff(g(eta), eta))^2-g(eta)*(diff(g(eta), `$`(eta, 2)))-1), lambda*(diff(g(eta), `$`(eta, 3)))+(diff(g(eta), `$`(eta, 2)))*f(eta)-g(eta)*(diff(f(eta), `$`(eta, 2)))]:

`<,>`(ODEs[])

Vector(2, {(1) = diff(f(eta), eta, eta, eta)+A^2+f(eta)*(diff(f(eta), eta, eta))-(diff(f(eta), eta))^2+beta*((diff(g(eta), eta))^2-g(eta)*(diff(g(eta), eta, eta))-1), (2) = lambda*(diff(g(eta), eta, eta, eta))+(diff(g(eta), eta, eta))*f(eta)-g(eta)*(diff(f(eta), eta, eta))})

(2)

LB, UB := 0, 1:

BCs := [`~`[`=`](([D(f), f, g, (D@@2)(g)])(LB), [1+B1*((D@@2)(f))(0), 0, 0, 0])[], `~`[`=`](([D(f), D(g)])(UB), [A, 0])[]];

[(D(f))(0) = 1+B1*((D@@2)(f))(0), f(0) = 0, g(0) = 0, ((D@@2)(g))(0) = 0, (D(f))(1) = A, (D(g))(1) = 0]

(3)

Params := Record(A = .9, B1 = .5, beta = .5, lambda = .5):

NBVs := [-((D@@2)(f))(1) = `C*__f`]:

Cf := `C*__f`:

Solve := module () local nbvs_rhs, Sol, ModuleApply, AccumData, ModuleLoad; export SavedData, Pos, Init;  nbvs_rhs := `~`[rhs](:-NBVs); ModuleApply := subs(_Sys = {:-BCs[], :-NBVs[], :-ODEs[]}, proc ({ A::realcons := Params:-A, B1::realcons := Params:-B1, beta::realcons := Params:-beta, lambda::realcons := Params:-lambda }) Sol := dsolve(_Sys, _rest, numeric); AccumData(Sol, {_options}); Sol end proc); AccumData := proc (Sol::{Matrix, procedure, list({name, function} = procedure)}, params::(set(name = realcons))) local n, nbvs; if Sol::Matrix then nbvs := seq(n = Sol[2, 1][1, Pos(n)], n = nbvs_rhs) else nbvs := `~`[`=`](nbvs_rhs, eval(nbvs_rhs, Sol(:-LB)))[] end if; SavedData[params] := Record[packed](params[], nbvs) end proc; ModuleLoad := eval(Init); Init := proc () Pos := proc (n::name) local p; option remember; member(n, Sol[1, 1], 'p'); p end proc; SavedData := table(); return  end proc; ModuleLoad() end module:

colseq := [red, green, blue, brown]:

Pc := B1 = .5, A = .1, beta = .5:

Ps := [B1 = .5, A = .1, beta = .5]:

Pv := [lambda = [.5, 1, 1.5, 2]]:

for i to nops(Ps) do plots:-display([seq(plots:-odeplot(Solve(lhs(Pv[i]) = rhs(Pv[i])[j], Ps[i][], Pc), [eta, theta(eta)], 'color' = colseq[j], 'legend' = [lhs(Pv[i]) = rhs(Pv[i])[j]]), j = 1 .. nops(rhs(Pv[i])))], 'axes' = 'boxed', 'gridlines' = false, 'labelfont' = ['TIMES', 'BOLDOBLIQUE', 16], 'caption' = nprintf(cat(`$`("\n%a = %4.2f, ", nops(Ps[i])-1), "%a = %4.2f\n\n"), `~`[lhs, rhs](Ps[i])[]), 'captionfont' = ['TIMES', 16]) end do

Error, (in dsolve/numeric/process_input) invalid argument: (B1 = .5)[]

 

``


 

Download 23.mw

@acer I'm sure that I was the original author of that code, but it looks like it's been put through the wringer a few times, extracting all comments and formatting, and incorrectly modified a bit, and converted to 2D Input and then back to plaintext.

@emendes 

Re: Homorbital: As I said, that one is my own definition. Perhaps the concept is defined elsewhere with a different name,

Re: counterparts: Don't think in terms of what is removed. Nothing is removed in this program that uses an Iterator. Rather, only items that are desired for the final output are included in the first place. That's why this program is memory conservative. The more conditions that you add, the more time it will use, but the less memory.

Re: arguments: Yes, I'll add those things as arguments.

You say that parms can change. How exactly? Will there always be two indices? Is the first index always $1..nIs?

Can you think of names for the 6 conditions? It would be useful to have key names to replace the 1-6 used as table indices.

Can you describe how the tables of sets for conditions 4-6 are created? It would be useful to create them programmatically from some higher concept if possible.

Here is a worksheet with all-new and faster code, including passing the arguments. I do not pass nIs because it can be derived from parms.

Symmetries.mw

Please test scenarios that gave errors before. I might have fixed some of them. I can't test because I don't have Linux or a machine with 36 processors.

@nm Maple does have a limited ability to translate some Mathematica. See ?MmaTranslator.

@crashoverwide Oops, I have any intermittent GUI bug that causes minus signs to show up as K, even on graph axes. I just automatically translate them in my mind; unfortunately that's how common this bug is. So I just realized that I posted such a graph above. So, here's the graph with regular minus signs.

Does anyone know anything about this bug. It happens on both my computers. Plus signs show as C, and there are many other weird transliterations.

@dingtianlidi Here's a bonus problem, since you seem interested. Using the same type of graphical analysis, we can convert the same integration region to polar coordinates. For this, consider a ray emanating from the origin. Through which curve does it enter the region? Through which curve does it exit the region?
 

f1:= sqrt(4*x-x^2): f2:= 2*sqrt(x):

plots:-display(
    plots:-shadebetween(
        piecewise(x<2, f1, x), f2, x= 0..4,
        color= magenta
    ),
    plots:-shadebetween(f1, x, x= 2..4, color= cyan),
    plot(x, x= 0..4, color= black, linestyle= dash, thickness= 2),
    gridlines= false
);

Polar:= f-> solve(eval(f, [x,y]=~ r*~[cos,sin](theta)), r):

g1:= Polar(y=f1);

0, 4*cos(theta)

g2:= Polar(y=f2);

0, 4*cos(theta)/sin(theta)^2

g3:= Polar(x=4);

4/cos(theta)

Int(r, [r,theta]=~ [g1[2]..g3, arctan(0,4)..arctan(4,4)]) +   #cyan
Int(r, [r,theta]=~ [g1[2]..g2[2], arctan(4,4)..arctan(4,0)]); #magenta

Int(r, [r = 4*cos(theta) .. 4/cos(theta), theta = 0 .. (1/4)*Pi])+Int(r, [r = 4*cos(theta) .. 4*cos(theta)/sin(theta)^2, theta = (1/4)*Pi .. (1/2)*Pi])

value(%);

32/3-2*Pi


 

Download Polar_conversion.mw

@Steve Roper It's safer to generalize that if-statement to 

if not [args]::list(numeric) then return 'procname'(args) fi

Thus it'll work correctly regardless of the procedure's parameter names.

You want to "simplify" the given irrational algebraic constant? It's irrational, so it doesn't get much simpler than it already is. It's equal to 5*2^(-59/6), but that's hardly any simpler. What the point? If you're solving a practical physics problem about dimensioned quantities (such as temperature), then exact answers (other than 0, infinity, etc.) are irrelevant. Just use evalf.

@vv Yes, it's an impressive program. I guess that the input inequalities must be converted to polynomial? If that's the case, you may be able to avoid most solve failures by using solve(..., parametric, real).

@JAMET Now verifying the identity is trivial by hand or with Maple:

is((z1+z2)^2/(z1*z2) = (z1/z2 + z2/z1 + 2));
             
true

 

@JAMET Let z1 and z2 both equal 1. Then the left side is 2 and the right side is 4. They're not equal.

@JAMET I think that you transcribed something wrong: (z1+z2)^2/(z1+z2) = (z1+z2), obviously.

@emendes Here is new code that processes all 6 conditions. You'll see that I added a mechanism to make it easy to exclude conditions (3 is excluded here) and to specify their order (I used 1,2,5,4,6).

I also included to kludge to workaround the error that occurs immediately after restart.

OrbitPartition:= module()
option
    `Conceptual author: emendes`,
    `Maple code author: Carl Love <carl.j.love@gmail.com> 2020-May-7`
;
uses It= Iterator;
local 
    #Problem setup:
    #==============
    nIs:= 3, #number of distinct first indices 
    Is:= {$1..nIs}, #set of first indices
    Pairs1:= combinat:-choose(Is,2),

    #fundamental set of lists (indexed variables without a stem):
    i, j, parms:= {seq(seq([i,j], j= 0..9), i= Is)}, 

    permutations_by_index:= [
        table([2= 3, 3= 2]), #permutations of 1st index
        table([2= 3, 3= 2, 5= 6, 6= 5, 7= 9, 9= 7]) #... of 2nd index
    ],

    #Exclusion conditions:
    #---------------------
    #Build a table of sets of the 2nd index grouped by their 1st index:
    ClassifyI:= proc(S)
    option cache;
        (op~)~(2, ListTools:-Classify(x-> x[1], S))
    end proc,

    op2:= proc(S) option cache; op~(2,S) end proc,
    
    TS_5t:= table([0=(), 1=1, 2=2, 3=3, 4=1, 5=(1,2), 6=(1,3), 7=2, 8=(2,3), 9=3]),
     TS_5:= proc(j) option remember; TS_5t[j] end proc,
     TS_6:= table([{1,2}={0,1,2}, {1,3}={0,1,3}, {2,3}={0,2,3}]),

    Condition:= table([
        1= proc(S) option cache; op~(1,S) = Is end proc,
        2= (S-> max(op2(S)) >= 4),
        3= (S-> nops({entries(ClassifyI(S))}) = nIs),
        4= (S->
            Condition[1](S) and
            not ormap(k-> ClassifyI(S)[k] subset [{0,1,4}, {0,2,7}, {0,3,9}][k], Is)
           ),
           5= (S-> TS_5~({op2(S)[]}) = Is),
           6= (S->
               Condition[1](S) and
               not ormap(p-> `union`(map2(index, ClassifyI(S), p)[]) subset TS_6[p], Pairs1)
              )
       ]),
    Conditions:= [1,2,5,4,6], #The order is for efficiency.
            
    #Decide whether a permutation's orbit is allowed to be represented:
    AllowOrbit:= S-> andmap(k-> Condition[k](S), Conditions),
    
    #=================
    #End of setup code

    #Set permutations to identity for unmentioned indices:
    T:= proc(T::table, k) option remember; `if`(assigned(T[k]), T[k], k) end proc,

    #Combine individual indices' permutations into a permutation function for lists:
    conds:= S-> map(x-> T~(permutations_by_index, x), S),

    #For any tuple of lists, compute its orbit under `conds`. Return a
    #representative of the orbit.
    Orbit:= proc(S)
    local r:= S, R, ao;
        do 
            if AllowOrbit(r) then R[r]:= () else return fi 
        until assigned(R[(r:= conds(r))]);        
        {indices(R, 'nolist')}[1] #Return lexicographic min as representative.
    end proc,

    Combos, #:=It:-Combination(nops(parms), tuple_size)

    DoChunk:= (rn::[posint,posint])->  #starting rank and number of iterates
        local 
            Combo_:= Object(Combos, 'rank'= rn[1]),
            has:= ModuleIterator(Combo_)[1], C:= output(Combo_)
        ;       
        (to rn[2] while has() do Orbit(parms[[seq(C+~1)]]) od),
        
    ModuleApply:= proc(tuple_size::posint, {sequential::truefalse:= false})
    local it;
        Combos:= It:-Combination(nops(parms), tuple_size);
        for it to 2 do #kludge to workaround error on 1st run after restart
            try
                return
                    `if`(sequential, map, Threads:-Map['tasksize'= 1])(
                        DoChunk,            
                        {It:-SplitRanks(Number(Combos))[]}
                    )
            catch:
                if it=1 then
                    printf(
                        "Error trapped, "
                        "so efficiency measures are invalid for this run.\n"
                    )
                else
                    error
                fi
            end try
        od
    end proc 
;
end module:

Please ask for definitions of any terms that you're not familiar with, especially those that are italicized.

From a mathematical perspective, tables and functions operating on finite sets are the same thing. This is probably obvious to you. It's only other nonmathematical factors that determine whether they're coded as tables or procedures. And a procedure with option remember or option cache is a hybrid of a table and a procedure.

Note that the tables in permutations_by_index, after being adjusted in T, are not merely functions of {$1..3} and {$0..9}; they are permutations (i.e., bijections with the same domain and range). I assume that this will always be the case---that you'd never have a similiar problem where these tables were not permutations. This observation is key to the rest of this article. Also assume for the rest of this article that the permutations operate on finite sets.

Theorem 1: The Cartesian product of permutations is also a permutaion, operating on the Cartesian product of the sets that the factor permutations operate on.
Proof: It's obvious; or you could likely find a proof in most combinatorics textbooks.

In our code, conds is precisely the Cartesian product of permutations_by_index, as adjusted by T.

Theorem 2: Let P be permutation on A. Let kA be the set of all k-combinations of A. Then the extension of to kA is also a permutation. In particular, applied to any k-combination is a k-combination.
Proof: I leave it for the reader.

Definition 1: Let P be a permutation and S a member of its domain. The orbit of S is the set {S, P(S), P(P(S)), ...}, which is necessarily finite. (This is a standard definition.)

Theorem 3: The orbits of a permutation defined on form a partition of A.
Proof: It's fairly obvious, and you could likely find it in numerous books.

Understanding all of the above is fundamental to understanding our program.

Definition 2: Let be a boolean predicate defined on a set A and let P be a permutation of A. We'll call C homorbital if it being true for one member of an orbit implies that it's true for all members. We call it necessarily homorbital if it's homorbital for any permutation. (These are my own definitions.) 

Now we're ready to address our specific conditions. Let A and B be sets with permutations PA and PB, and let P be the natural permutation defined on the set SS of k-combinations of the Cartesian product AxB. Let C1 be the condition defined on SS by C1(S) is true iff the set of first elements of S is all of A. (This is our condition 1). Then C1 is necessarily homorbital.
Proof: PA(A) = A for any permutation PA, by the definition of permutation.

With the same setup, let C2 be the condition that the maximum of the second elements of S is greater than or equal to 4. Then C2 is homorbital for our specific permutations.
Proof: Our permutation of the second position doesn't exchange any element greater than or equal to 4 with one less than 4.

Our condition 3 is necessarily homorbital. I won't bother proving this because you're no longer using it.

Conditions 4 and 5 are not homorbital, but you said that their being false necessarily disqualifies their counterparts, which amounts to the same thing as being homorbital for the purpose of our program.

What about condition 6? Does it also disqualify the counterpart(s)?

 

First 198 199 200 201 202 203 204 Last Page 200 of 709