Carl Love

Carl Love

28055 Reputation

25 Badges

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

MaplePrimes Activity


These are answers submitted by Carl Love

The ".mtx" (aka MatrixMarket) file format is a matrix format, not necessarily a graph format. Since the entire mathematical structure of a graph is encoded in its weight matrix (or adjacency matrix if you don't care about weights), this doesn't present any serious difficulties; you just need the tiny additional step of creating a Matrix when going from Graph to ".mtx" file or vice versa.

The example file that you linked is a sparse representation likely suitable for the vast majority of graphs. It's not necessary that you understand the representaion in order to use this. Only nonzeros are represented. The 1st column is the row numbers, the 2nd is the column numbers, and the 3rd is the weights. As you said, the 3rd column could be all 1s, making it an adjacency matrix.

To go from a Maple Graph G to a ".mtx' file, use

ExportMatrix("my_graph.mtx", GraphTheory:-AdjacencyMatrix(G));

To go in the other direction, use

G:= GraphTheory:-Graph(ImportMatrix("my_graph.mtx"));

The GraphTheory:-Graph command doesn't care whether you pass it an adjacency matrix or a weight matrix; either case will be handled correctly without you needing to specify which one it is.

You have three errors:

  1. The syntax of pdsolve(..., numeric) requires you to change {ICs}, BCs to 
    {ICs, BCs[]}
  2. You have the expression exp(-1000*(x-l)^2) several times. What is x? I guess it should be y.
  3. You have 8 BCs; it wants 6. I removed the 1st and 3rd (arbitrarily).

Making these changes, the pdsolve command itself runs without error. You may still encounter numeric errors when you try to get output. Indeed, I think that it's likely that you will.

I'd use grid= [33, 129]. So, that's 33 r values and 129 theta values. It takes 50 seconds on my computer:

CodeTools:-Usage(
    plot3d(
        [A,B,C], r= 1/5..4/5, theta= -Pi..Pi, view= [(-8..8)$3],
        shading= zhue, grid= [33, 129], style= surface, axes= none
    )
);

You have Vectors named err1 and err2 to hold your errors. The command LinearAlgebra:-Norm is akin to an "absolute value" for vectors---it's a way to assess their magnitude as a nonnegtive real number. So, use 

if LinearAlgebra:-Norm(err1, 2) < tol then ...

The 2 indicates the 2-norm (aka Euclidean norm). Other common norms are 1 (sum of abs of elements) and infinity (max of abs). See ?LinearAlgebra,Norm.

Like this:

local D:= Matrix(A, shape= diagonal);
(L,U):= map2(Matrix, A-D, shape=~ [triangular[lower], triangular[upper]])[];

By setting

Digits:= 50;

and adding option complex to your fsolve command (and no other changes), I get this solution in under 2 seconds:

{A = -2.7553365135418814642586082436429575890825402826031,
B = -0.70285804987973303586180028708027467941012949957141}

Note that this solution is real and didn't require initial specification of intervals.

Using the 3rd-party package DirectSearch (available for free download from the Maple Applications Center), and using its command SolveEquations with option AllSolutions, I get 2 solutions: the same one that I show above, and the one reported by Rouben. I'm not confident that this program, or any numeric method, can guarantee the number of solutions to a system of more than 1 non-polynomial equation, even when all univariate restrictions are meromorphic (as is the case here).

@Rouben Rostamian  The procedure can be shortened like this:

tangent_plane:= P-> [P, (sign*primpart@sort@`.`)((<1,3,5>, <x,y,z>)-~'<P>') = 0]:
tangent_plane~(L);

@Laurenso: For your purposes here, the single command sign is equivalent to your signum@lcoeff, but is faster and more robust. It is designed for this type of operation (exact work on polynomials); signum is not.

Here is my version, which uses lists of integers rather than arrays. It's a very short code. The run time is just slightly longer than dharr's.

restart
:
nextone:= L-> local i:= 2, j, n:= nops(L);
    ["", (do for i from (j:= i) + 1 to n while L[i]=L[j] do od; i-j, L[j] until i > n)]
:
LookAndSay:= (n::posint, start::nonnegint)-> local L:= ["", start];
    <start, (to n-1 do parse(cat((L:= nextone(L))[])) od)>
:
CodeTools:-Usage(LookAndSay(50,1));
memory used=309.93MiB, alloc change=27.45MiB,
cpu time=3.52s, real time=3.52s, gc time=281.25ms

Regarding the timing of your Maple code (the StringTools-based code in your Question): You can't blame the compiled StringTools code. The vast majority of the time is spent doing garbage collection at the Maple kernel level, not in the compiled code. This is due to the large number of temporary substrings that you create.

You just needed some minor syntax corrections:

R:= plots:-implicitplot(
    [x = -3/4, x = -1/3, x = 1/3, x = 3/4, y = -3/4, y = -1/3, y = 1/3, y = 3/4], 
    x = -3 .. 3, y = -3 .. 3, color = red, scaling = constrained, gridrefine = 3
);
z:= x + y*I:
plottools:-transform(unapply(evalc([(Re, Im)(1/z)]), x, y))(R);

The circles don't show as complete circles because the lines are not complete (infinitely long) lines; they're line segments.

To get a single solution, you just need to make 3 trivial changes:

  1. Correct beta_1 = 0.2 to beta_1:= 0.2;
  2. Correct all pi to Pi;
  3. Change solve to fsolve.

If you need all solutions, or all real solutions, it's a little trickier.

The example code that you show

f[N1_]=sum[1/(n^3*sin^2[n]), {n,1,N1}];
DiscretePlot[f[x], {x,0,400},PlotRange->All];

is Mathematica, not Maple.

For Maple, you can do a continuous plot of it as

plot(eval(sum(exp(-lambda)*lambda^n/n!, n= 0..N1), lambda= 15.4), N1= 0..50);

If you insist on a discrete plot, simply replace N1 with trunc(N1) inside the sum:

plot(eval(sum(exp(-lambda)*lambda^n/n!, n= 0..trunc(N1)), lambda= 15.4), N1= 0..50);

You need to set your working directory to something that your OS gives you permission to write a file to. This can be done with the currentdir command. For example, on my computer, this works:

currentdir("/users/carlj/desktop");
ExportMatrix(matlabData, A, target= MATLAB, format= rectangular, mode= ascii);

Let's suppose that:

  1. is any of the arguments of exp that occur in your expressions.
  2. You're able to write a procedure TransA that takes any such and returns what it should be replaced with.
  3. Expr is any expression or list, set, etc., of expressions containing some exp functions that you want to change.

Then, this command will do what you want:

subsindets(Expr, specfunc(exp), exp@TransA@op)

Your procedure TransA must return something for any A that it's passed. If it simply returns an A unchanged, that's fine.

If A itself contains exp subexpressions, there's no need for TransA to transform their arguments, because the subsindets will always process them in order from most deeply nested to least deeply nested.

There's nothing special about exp here. The above process will work for any function of a single argument. With slight modifications, it'll also work for multi-argument functions.

Since the probability is very high that the desired condition will be met anyway, I think that the best strategy is to check that the condition is met, and, if it's not, regenerate the matrix. Like this (using 1-D input):

do
    A:= LinearAlgebra:-RandomMatrix((10,10), 'generator'= -10..10)
until andseq(rtable_scanblock(A, [i, ..], 'NonZeros') > 2, i= 1..10);

You can also generate the rows satisfying the condition one at a time and catenate them, like this (1-D input) (I've added the density= 0.3 to give the problem a more-realistic level of difficulty):

A:= <(
    to 10 do
        do
            R:= LinearAlgebra:-RandomVector(
                10, 'generator'= -10..10, 'density'= 0.3
            )^+
        until rtable_scanblock(R, [], 'NonZeros') > 2
    od
)>;

 

You wrote:

  • The bases that exit for one selected period (your proc) are pairwise conected to another base in the set.  I can see this in the plots.

    So  it's possible to run only the half of the loop with some extra calculations for pairfinding. I think this may reduce the time to 60% of running the fulll loop.

Yes, both of those observations are essentially correct, but actually the time reduction is usually much greater. The "orbits" of a given base with respect to a given period are not just pairs but have size totient(period). So, for example, for period = 11, the bases come in groups of size totient(11) = 10. Knowing any one base in that orbit, the other 9 are trivial to generate. Here is a program for it:

restart:
Periods:= module()
export
    #prime factorization in convenient list form (see ?ifactors):
    Ifactors:= (n::posint)-> (thisproc(n):= ifactors(n)[2]),

    #prefactored form of the totient, making it more efficient to compute the 
    #Carmichael number at the same time:
    TotSeq:= (n::posint)-> 
    local p; 
        (thisproc(n):= [seq]((p[1]-1)*p[1]^(p[2]-1), p= Ifactors(n))),

    #Euler's totient function:
    Totient:= (n::posint)-> (thisproc(n):= mul(TotSeq(n))),

    #Carmichael's lambda function:
    Carmichael:= (n::And(posint, Not(1)))-> 
    local T:= TotSeq(n);
        (thisproc(n):= ilcm(`if`(n::even and T[1]>2, [T[1]/2, T[2..][]], T)[])),

    #list of prime divisors of n:
    PrimeDivisors:= (n::posint)-> (thisproc(n):= op~(1, Ifactors(n))),

    #Carmichael number reduced by its prime divisors:
    PreCarmichael:= (n::And(posint, Not(1)))->
    local L:= Carmichael(n);
        (thisproc(n):= iquo~(L, PrimeDivisors(L))),

    (* Find any element of order Carmichael(n). The optional table Ex contains
    elements to NOT check. The optional 3rd argument St is the element 
    to start with. *) 
    Generator:= proc(
        n::And(posint, Not({1,2})), Ex::table:= table('sparse'), St::posint:= 2
    )
    local E:= PreCarmichael(n), a, e;
        for a from St to n-1 do
            if Ex[a]=0 and igcd(a,n)=1 and andseq(a&^e mod n <> 1, e= E) then 
                return a 
            fi
        od;
        FAIL
    end proc,

    #Find the elements of order ord. The optional 3rd argument NG is the
    #maximum number of generators to use.
    Bases:= proc(
        n::And(posint, Not({1,2})),
        ord::And(posint, Not(1)),
        NG::posint:= infinity
    )
    local 
        Ex:= table('sparse'), L:= Carmichael(n), LP:= iquo(L, ord), 
        E:= select(igcd=1, [$1..ord-1], ord), EG:= select(igcd=1, [$1..L-1], L), 
        G, e, g:= 1, r, ng:= 0, h, TT:= Totient(Totient(n))
    ;
        if irem(L, ord) <> 0 then return {} fi;
        E:= E[2..] -~ E[..-2];
        EG:= EG[2..] -~ EG[..-2];
        {
            to NG while ng <> TT and (g:= Generator(n, Ex, g+1)) <> FAIL do
                G:= e-> (thisproc(e):= g&^e mod n);
                Ex[(h:= g)]:= 1; ng++;           
                for e in EG do 
                    if Ex[(h:= h*G(e) mod n)]=0 then Ex[h]:= 1; ng++ fi 
                od;
                r:= g&^LP mod n;
                G:= e-> (thisproc(e):= r&^e mod n); 
                ((h:= r), seq((h:= h*G(e) mod n), e= E))
            od
        }
    end proc  
;
end module
:
den:= 11842585:  #same number that you used
NumberTheory:-Divisors(Periods:-Carmichael(den));
{1, 2, 3, 4, 6, 7, 11, 12, 13, 14, 21, 22, 26, 28, 33, 39, 42, 
  44, 52, 66, 77, 78, 84, 91, 132, 143, 154, 156, 182, 231, 273, 
  286, 308, 364, 429, 462, 546, 572, 858, 924, 1001, 1092, 1716, 
  2002, 3003, 4004, 6006, 12012}

bases11:= CodeTools:-Usage(Periods:-Bases(den, 11));
memory used=287.81MiB, alloc change=0 bytes, 
cpu time=4.28s, real time=3.79s, gc time=703.12ms

bases11 := {207496, 299716, 338141, 484156, 514896, 561006, 
  568691, 737761, 853036, 1029791, 1098956, 1237286, 1367931, 
  1452466, 1544686, 1590796, 1598481, 1752181, 1767551, 1844401, 
  2105691, 2128746, 2159486, 2267076, 2282446, 2359296, 2482256, 
  2620586, 2781971, 2797341, 3135481, 3143166, 3189276, 3312236, 
  3389086, 3573526, 3658061, 3673431, 4026941, 4188326, 4418876, 
  4457301, 4541836, 4603316, 4634056, 4733961, 4933771, 5202746, 
  5248856, 5448666, 5487091, 5663846, 5717641, 5733011, 5871341, 
  5963561, 6086521, 6148001, 6224851, 6247906, 6401606, 6601416, 
  6662896, 6747431, 6793541, 6901131, 7031776, 7208531, 7254641, 
  7277696, 7308436, 7431396, 7546671, 7631206, 7723426, 7777221, 
  8207581, 8292116, 8307486, 8338226, 8445816, 8660996, 8722476, 
  8799326, 8807011, 8822381, 8960711, 8976081, 9052931, 9091356, 
  9175891, 9268111, 9314221, 9321906, 9337276, 9368016, 9490976, 
  9606251, 9690786, 9752266, 9783006, 9852171, 9882911, 10082721, 
  10205681, 10397806, 10505396, 10636041, 10782056, 10812796, 
  10858906, 10912701, 11020291, 11035661, 11296951, 11373801, 
  11550556, 11627406, 11665831, 11811846}

 

First 21 22 23 24 25 26 27 Last Page 23 of 395