Carl Love

Carl Love

26488 Reputation

25 Badges

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

MaplePrimes Activity


These are answers submitted by Carl Love

Like this:

DrawSpaceCurve:= (
    fnc::And(list(algebraic), 3 &under nops),
    Evals::list(name= complexcons), 
    xr::(name=range(realcons))
)->
    plots:-spacecurve(eval(fnc, Evals), xr, _rest)
:
DrawSpaceCurve([x, y, x^2*y], [y=2], x= -5..5, color= blue, axes= frame);

 

Verify:= proc(n::And(odd, positive, Not(1)))
local s;
    for s from 0 to isqrt(iquo(n,2))+1 do
        if isprime(n-2*s^2) then return true fi
    od;
    false
end proc
:
for n from 3 by 2 to 9000 do 
    if Verify(n) then next 
    else printf("Failure for n=%d", n); break
    fi
od
:
Failure for n=5777

To see all failures, simply remove break.

The above can be made a little bit more efficient, at the expense of a bit of clarity. I'll post a followup. Note that my code doesn't use any sequences or lists; however, that approach isn't bad either. Though the one thing that I don't like about it is that it requires foreknowledge of the upper limit (the 9000 in this case).

Note that there are far more primes than squares. Thus, in Verify, it's more efficient to iterate over the squares and check primality than to iterate over the primes and check squareness.

The contourplot command works for parameterized surfaces. You just need to make the function a 3-list instead of a 3-vector:

Ennep:= [u - u^3/3 + u*v^2, v - v^3/3 + v*u^2, u^2 - v^2]:
plots:-contourplot(Ennep, u= -2 .. 2, v= -2 .. 2);

The list form is more convenient for most purposes, including your original plot3d.

In your worksheet, you've made two definitions of (or assignments to) A:

  • (A1): A:= n -> w*int(V*cos(w*n*t), t = 0 .. pw)/Pi:
  • (A2): A:= Vector(Ns, n -> F(n)):

I think that the best solution to your problem woud be to use some variable other than A for the second definition.[*1] That completely corrects the problem. However, it's also very important that you understand the reason why your original code led to the erroneous result: It's a subtle order-of-evaluation issue, and many such issues occur in Maple code (and also other symbolic-computation languages). In between definitions (A1) and (A2), you've made two definitions that directly refer to A:

  • (F1): f(t):= a/2 + sum(A(n)*cos(w*n*t), n= 1..100) + sum(B(n)*sin(w*n*t), n= 1..100);
    It's important to view the output of the above line to see that Maple has interpretted it as if you had entered
    f := t -> a/2 + sum(A(n)*cos(w*n*t), n= 1..100) + sum(B(n)*sin(w*n*t), n = 1..100);
    I recommend that when you define a function, you use the form f:= t-> ... rather than f(t):= ... because Maple will not always interpret the latter as the former, and whether this reinterpretation happens is crucial for the rest of this discussion.
  • (F2): F:= n -> sqrt(A(n)^2 + B(n)^2);

The fact that you've made two definitions ((F1) and (F2)) with direct references to A is immaterial here; what I'm about to say would still apply if you had made one or any other number. The crux of the matter is which definition of A, (A1) or (A2), is being used. It would be natural for you to think that (A1) is used because (F1) and (F2) are defined before (A2) is defined. But this isn't what happens. When a function is defined with the arrow -> (or with proc (see below)), its symbols (such as AB, and w---and even the library symbols sin, cossqrt, and sum; and even the kernel symbols +, /, *, and ..!) use whatever definition is in effect at the time that the function (f or F in this case) is used rather than at the time that it was defined.  Thus, your usage F(5) that returns 0.7... is using (A2). This interpretation of the symbol is called its evaluation, and I'll assume that you can now understand why this overall crucial issue is called the order of evaluation. There are several ways to change the order of evaluation; however, since it's not necessary for you to do so in order to correct your problem, and since I've just taught you a huge paradigm shift into symbolic computation, I won't overload this Answer`those concepts.[*2]

There's another issue that doesn't currently cause a problem, but it could easily become a problem later. The definition of given at (F1) relies on n not having any definition at the time that f is used. The solution to this is to make local to f. This is a concept that applies to almost all computer languages, whether or not they're symbolic. Change the definition (F1) to

f:= proc(t)
local n;
    a/2 + sum(A(n)*cos(w*n*t), n= 1..100) + sum(B(n)*sin(w*n*t), n = 1..100)
end proc;

If you use 1D input (aka Maple Input or Maple Notation), the above can be abbreviated to

f:= t-> local n; a/2 + sum(A(n)*cos(w*n*t), n= 1..100) + sum(B(n)*sin(w*n*t), n = 1..100);

When this is done, the n inside f will not be affected by any change made to any n outside of f. Indeed, it becomes a completely different variable that just happens to have the same name. This need for local does not occur elsewhere in your code because variables that appear on the left side of an arrow -> are automatically local in a sense (technically, they're called parameters, which is a slightly different concept than local).

[*1] There are other situations, however, where redefining a symbol is the best thing to do to maintain the human readability of the code.

[*2] The best way to change the order of evaluation in this case is a command called unapply. After you've digested the rest of material in this Answer, you should read its help page ?unapply.​​​​​​

For any given k, that question (the number of times required) can be answered by 

GroupTheory:-PermOrder(Perm(eq_arrangement(k)));

There's likely an easy relationship between k and the returned order, but I haven't spotted it yet.

The only "edo" I know is a castle in Japan. Did you mean "ode"? If so, the command is dsolve. For a numeric solution, include the keyword numeric in the arguments.

It just needs small changes. The procedure arguments must be numbers, not a list of numbers.

TV := proc(p1,p2)
  sol(parameters=[args]):
  sol(10):
  return sol(eventfired=[1])[];
end proc:

TV(1$2);
ranges := 0.8..1.2, 0.8..1.2:
Optimization:-Maximize(TV, ranges);

 

@one man Thanks to your worksheet, I've finally taken a good look at Draghilev's method. I noticed that you use extremely unnecessarily verbose syntax. I guess that this is because you don't know the easier ways. I've revised your worksheet. I hope that you find this useful. Here are the code and comments. For the prettyprinted output, see the attached worksheet.

restart
:
VC:= VectorCalculus: #package

#Parameterize the curve of intersection of two spheres.
V:= [x,y,z]: #coord vars
(R,C):= ([1, 0.7], [[0, 0, 0], [0.5, 0.25, 0.2]]): #radii & centers
F:= (c-> `+`(((V-~c)^~2)[]))~(C) -~ R^~2; #spheres
             [ 2    2    2      
        F := [x  + y  + z  - 1, 

                   2             2            2       ]
          (x - 0.5)  + (y - 0.25)  + (z - 0.2)  - 0.49]

n:= [$1..nops(V)]:
deqs:= (-1)^~n*~diff(V(s),s) =~
    eval(
        (i-> VC:-Jacobian(F, V[[({n[]} minus {i})[]]], determinant)[2])~
           (n), 
        V=~ V(s) #change coord vars to functions
    )
:
print~(deqs):
                 d                              
              - --- x(s) = -0.8 y(s) + 1.00 z(s)
                 ds                             

                 d                             
                --- y(s) = -0.8 x(s) + 2.0 z(s)
                 ds                            

                 d                              
              - --- z(s) = -1.00 x(s) + 2.0 y(s)
                 ds                             

#We need one point on the curve of intersection for initial conditions. 
#To obtain the point, add an arbitrary constraint. Here, I use x=y. 
#If no solution is obtained, I'd try another constraint. 
ics:= V(0)=~ eval(V, fsolve([F[], x=y]));
       ics := [x(0) = 0.6653163290, y(0) = 0.6653163290, z(0) = -0.3386862336]

sln:= dsolve({deqs[], ics[]}):
print~(evalf[3](sln)): 
     x(s) = 0.612 + 0.368 sin(2.38 s) + 0.0536 cos(2.38 s)
     y(s) = -0.511 sin(2.38 s) + 0.359 cos(2.38 s) + 0.306
     z(s) = -0.583 cos(2.38 s) - 0.281 sin(2.38 s) + 0.245

s_rng:= s= 0. .. 2*Pi/2.38:
plots:-display(
    plottools:-sphere~(C, R, color=~ [green, yellow])[],
    plots:-spacecurve(eval(V(s), sln), s_rng, color= red, thickness= 3),
    transparency= 0.5, axes= normal, scaling= constrained
);


Download worksheet: Draghilev.mw
 

plot(<K[..,1] | y~(K[..,1]) - K[..,3]>, view= [0.35..0.4, DEFAULT]);

As discussed in my Reply to Preben, the following code addresses all the issues raised thus far, and does so in a way that's transparent to the end user.

restart:

unprotect(
    LinearAlgebra:-Eigenvectors, LA_EV_orig, LA_EV_inner, LA_EV_sort
):
LA_EV_orig:= eval(LinearAlgebra:-Eigenvectors):

SetOrder:= (x,y)-> x={x,y}[1]
:
LA_EV_sort:= (R::list)->
#R is [output] from LinearAlgebra:-Eigenvectors. The eigenvalues and 
#corresponding vectors are sorted in "set order".
    if R=[] then return 
    elif R=[1]::Vector then 
        local p:= sort(R[1], SetOrder, ':-output'= ':-permutation');
        if nops(R) > 1 and R[2]::Matrix then
            (R[1][p], R[2][.., p], R[3..][])
        else
            (R[1][p], R[2..][])
        fi
    elif R[1]::list then 
        (sort(R[1], SetOrder@op@curry(op,1)~@`[]`), R[2..][])
    else R[]
    fi
:
LA_EV_inner:= proc(AC::seq(specfunc(_Inert_MATRIX))) 
option cache;
    (LA_EV_sort@[LA_EV_orig])(FromInert~([AC])[], _rest)
end proc
:
LinearAlgebra:-Eigenvectors:= overload([
    proc(
        AC::seq(Matrix(square, algebraic)), 
        O::seq({name, name=anything}):= (),
        $
    )
    option overload;
        LA_EV_inner(ToInert~([AC])[], O)
    end proc,

    LA_EV_orig
]):
protect(LinearAlgebra:-Eigenvectors, LA_EV_orig, LA_EV_inner, LA_EV_sort):
    
#Example:
Sy:= Matrix([[0,-I,0],[I,0,-I],[0,I,0]])/sqrt(2):
for k to 5 do (ev||k, EV||k):= LinearAlgebra:-Eigenvectors(Sy) od:
ev1, EV1;

#Prove that they're all equal:
andmap(LinearAlgebra:-Equal, {ev||(2..5)}, ev1), 
    andmap(LinearAlgebra:-Equal, {EV||(2..5)}, EV1);

                       [-1]  [   -1     1     -1    ]
                       [  ]  [                      ]
                       [ 0]  [   (1/2)         (1/2)]
                       [  ], [I 2       0  -I 2     ]
                       [ 1]  [                      ]
                             [   1      1      1    ]
       
                                 true, true

#Request list-form output:
LinearAlgebra:-Eigenvectors(Sy, output= list);
 

 

The commands (a) evalb(A=B), which compares expressions' identities, and (b) is(A=B), which compares expressions mathematically, possibly under assumptions, are intentionally different, and this is not due solely to efficiency considerations. It wouldn't be "more robust" to change (a) in any way that makes it more like (b). Indeed, it would be far less robust, because (b) must necessarily (see Richardson's Theorem) FAIL is some cases.

The poorly named and mostly redundant command LinearAlgebra:-Equal maps operation (a) over arrays. It's trivial to write a similar command that maps (b) over arrays. There are many possibilities, some shown by the other respondents. My first choice is

Is:= (A::rtable, B::rtable)->
    LinearAlgebra:-Equal(0~(A), 0~(B), _rest) and 
    andmap(is, simplify~(A-~B)=~ 0)
:

This procedure behaves reasonably even if A and B have mismatched sizes or shapes.

Several programming points:

1. The break statement (up until Maple 2021) breaks out of the innermost loop that contains it. Maple 2021 introduced the multi-level break, which is described on pages 1-5 of this document: https://www.maplesoft.com/products/maple/new_features/Maple2021/PDFs/Language.pdf

However, Kitonum has already given the better solution for this particular example: return.

2. The print statement is never an appropriate way for a procedure to return a value. It is only for displaying supplementary information.

3. A procedure should not rely on a particular package having been loaded with with. Nor should a with command ever be used in a procedure. Instead, use a uses clause (example below).

4. The problem of loops nested to an arbitrary depth can usually be solved by using a Cartesian product, like this:

Conremovedges:= proc(
    G::GRAPHLN, 
    Paircrossedges::And(
        {set,list}({set,list}(set)),
        identical({[2,2]}) &under ((nops~@[op])~@{op})
    )
)
uses GT= GraphTheory, It= Iterator;
    for local p in It:-CartesianProduct(Paircrossedges[]) do       
        if 
            (GT:-EdgeConnectivity@GT:-DeleteEdge)(
                G, (local dedges:= {seq}(p)), inplace= false
            ) = 3
        then return dedges
        fi
    od;
    "not found"
end proc
:
Conremovedges(
    GraphTheory:-CompleteGraph(6), 
    [[{1,6},{3,5}],[{2,5},{1,4}],[{2,6},{3,4}]]
);
                    {{1, 6}, {2, 5}, {2, 6}}

​​​​​​

I believe that that is a result of an anti-piracy measure included by Maplesoft in specific version(s) of Maple. It's possible that the measure is triggered by false positives. It has been reported here on MaplePrimes about a dozen times. The solution is to reinstall Maple with correct licensing information. My guesses about all this are based solely on what's been reported here on MaplePrimes, and thus they're highly speculative.

By the way, "Problem with output" is a very poor title for a Question here because it could describe the majority of things that users report.

You need to change self to _self. This is mentioned on page 6 of this PDF: https://www.maplesoft.com/products/maple/new_features/Maple2021/PDFs/Language.pdf

Suppose that your solution is stored in Sol. Then do

NewSol:= eval(Sol, BesselK= ((x,y)-> (2/y)^x/GAMMA(x)));

Note that GAMMA is all capitals.

You might find the command MultiSeries:-asympt useful.

First 56 57 58 59 60 61 62 Last Page 58 of 382