Carl Love

Carl Love

28015 Reputation

25 Badges

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

MaplePrimes Activity


These are answers submitted by Carl Love

One way:

ST:= table((parse@lhs=rhs)~(S)):
MS:= [mu, sigma]:
MS =~ index~(ST[2], MS); rhs~(%)[];

 

In cases where the known function f(x) has certain relatively simple forms (such as polynomial), a usable solution for unknown g(x) can be obtained like this:

restart:
f:= x-> 2*x+5:
h:= f*g:
h2:= (D@@2)(h)(x);
                 h2 := 4 D(g)(x) + (2 x + 5) @@(D, 2)(g)(x)

#We need to "hide" the symbolic derivatives because they 
#will confuse the solve command:
S:= {D(g)(x)= g1, (D@@2)(g)(x)= g2}: 
h2i:= subs(S, h2);
                         h2i := 4 g1 + (2 x + 5) g2

sol:= solve({h2i < 0, g1 > 0, g2 < 0}, x, parametric);

          sol :=  / [[  4 g1 + 5 g2    ]]                         
                  | [[- ----------- < x]]      And(g2 < 0, 0 < g1)
                 <  [[     2 g2        ]]                         
                  |                                               
                  \          []                     otherwise     

#Thus, the x-interval(s) of concavity for h are where this inequality
#is true:
Concavity_condition:= op(2, sol)[][];
                                           4 g1 + 5 g2    
                  Concavity_condition := - ----------- < x
                                              2 g2        

#Since we know g2 < 0, that can be simplified to
Concavity_condition2:= %f(x)*g2 + 4*g1 < 0
                 Concavity_condition2 := %f(x) g2 + 4 g1 < 0

#In this case, I did that simplification by hand. 
#Note that I used % to make f(x) inert.

 

If you can make that table into a 3-column Matrix (I'll call it A), then do

plots:-pointplot3d(A, connect)

Or if you can make three VectorXYZ, then do

plots:-pointplot3d(X, Y, Z, connect)

Or if you can make the table into a list of 3-element sublists, then do

plots:-pointplot3d(L, connect)

Like this:

restart:
interface(prettyprint= 1):

#Your original input:
ode:= diff(y(x),x) = 2*(2*y(x)-x)/(x+y(x)):
ic:= y(0)=2:

#Transform to inverse function:
ode1:= 1/diff(x(y),y) = subs(y(x)= y, x= x(y), rhs(ode));
ic1:= (op@lhs=x@rhs)(ic);

#dsolve(implicit), then transform back:
sol:= subs(x(y)=x, y= y(x), dsolve({ode1,ic1}, 'implicit'));

                                 1       2 (2 y - x(y))
                      ode1 := -------- = --------------
                               d            x(y) + y   
                              --- x(y)                 
                               dy
                      
                               ic1 := 0 = x(2)

               /2 x - y(x)\       /x - y(x)\                              
   sol := -3 ln|----------| + 2 ln|--------| - ln(y(x)) + I Pi + ln(2) = 0
               \   y(x)   /       \  y(x)  /                              

#A little bit of algebra shows that that's equivalent to your
#"book solution":
exp(combine(lhs(sol), symbolic)) = exp(rhs(sol));
                                           2    
                               2 (x - y(x))     
                             - ------------- = 1
                                           3    
                               (2 x - y(x))     

 

Here's a procedure for it. It counts all the indices in one pass:

Count:= (C::{name, name^posint, `*`({name, name^posint})})->
    (Statistics:-Tally@op~@[op]@indets)(
        subsindets(
            C, And(name, indexed)^posint,
            e-> `*`(
                'convert(op([1,0], e), `local`)[(op@op)(1,e)]' 
                $ op(2,e)
            )
        ),
        indexed
    )
:
C:=x*u[]*u[1,2,2]^3*u[1,1,2]:
Count(C);
                         [1 = 5, 2 = 7]

 

As nm and acer have explained, the difference that you're seeing is not due to a difference between map and ~ but due to the difference between matrix/vector multiplication in the old linalg package vs the modern A.B. However, there are two differences between map and that average users should be aware of:

  1. can only map over lists, sets, rtables (which includes Vectors, Matrixes, and Arrays), and tables (I'll call these the neutral containers for the rest of this post); whereas map can map over almost any composite structure (functions, polynomials, etc.). Function arguments that aren't one of those types get treated as singletons. Example 1a: Let L be a list of pairs with each pair being a list; then op~(2, L) returns the same as map2(op, 2, L). Example 1b: The list of prime factors of any integer n is returned by op~(1, ifactors(n)[2]). In both these examples, the first argument of opor 2, is being ignored by (and still seen by op).
  2. map only maps over one object; whereas ~ maps simultaneously over all neutral containers that appear as arguments as long as their sizes (and sometimes shapes) conforms. Example 2: f~([1,2], 3, {a, b}).

 

Any plots of the same number of dimensions (either 2 or 3) can be combined into a single plot with the command plots:-display. The plots don't need to have anything in common other than the number of dimensions. So, you could do this:

p1:= plots:-arrow(...);
p2:= plots:-arrow(
...);
p3:= plots:-arrow(
...);
plots:-display(p1, p2, p3, 
...);

where the ... in the arrow commands represents exactly the arguments in your posted worksheet, and the ... in the display command represents possible other options that would apply to the plot as a whole (such as scaling= constrained). The above could also be done as

plots:-display(
    plots:-arrow(
...),
    plots:-arrow(
...),
    plots:-arrow(
...),
    ...
);

I didn't see any vector field defined in your worksheet; however, if there had been, you could use one the several commands for plotting vector fields (such as plots:-fieldplot3d) and include that in the display command.

All the information is in the paper, but the portion that you showed only gives this obviously periodic function:

x:= t-> exp((cos(t)-cos(T))/4);

They say T=Pi/2t>T, and gamma=0.4 (although your screenshot doesn't show where gamma is used). We can plot x like this:

T:= Pi/2:
plot(x(t), t= 1..24, view= [1..24, -0.2..1.8]);

All that you need is 

[entries](op~(1, ListTools:-Classify(O-> O:-sol, SOL)), nolist);

Even behind the scenes, Classify only does a single pass through the list.

Also, you can eliminate the problematic building of a list one element at a time yet still use a list like this:

restart;
ode_solution_type:= module()
    option object;
    export sol:= (), stuff;
end module:

SOL:= [
    for n to 5 do
        tmp:= Object(ode_solution_type);
        tmp:-sol:= y(x)=`if`(n=3, 1, 0);
        tmp:-stuff:= rand();
        tmp
    od
]:

As I told you a day or two ago:

Theorem: If the prime factorization of n is n = p1^e1*p2^e2*...*pb^eb, then its number of proper divisors is 
(e1+1)*(e2+1)*...*(eb+1) - 1.

Note that this number depends only on the exponents, not the prime factors.

This can be encoded in Maple 13 as a one-line procedure:

NPropDivs:= (n::posint)-> `*`((op~(2, ifactors(n)[2])+~1)[]) - 1:

Or you can use the command

numtheory[sigma][0](n) - 1;

Or Tom's method can be extended to arbitrary positive integers as

combinat:-numbcomb((`$`@op)~(ifactors(n)[2])) - 1;

Although it's kind of roundabout, it's still orders of magnitude faster than what you're doing. 

As we've discussed before, you are a Moderator. The editing privilege is reserved for Moderators, who are presumed to be responsible enough to use it with discretion and moderation. 

In addition to the situations mentioned by dharr, there are numerous others where editing should be done. All of these are commonly occurring, not things I'm making up:

  1. There are a huge number of blank lines in the post, often so that the entire first screen is blank.
  2. The entirety of the post is typed into the title box and the main body is blank or nearly blank.
  3. The product specifier is wrong or absent or all possible boxes have been rudely checked.
  4. The post has been misclassified as Post, Question, Answer, Reply.

Several high-profile crowd-sourced information sites, such as StackExchange and Wikipedia, use a system where contributors gradually gain trust and editorial abilities. 

The problem that you have displaying n:= 2^(2^26)+1 is only a limitation of the prettyprinted GUI display; it's not a computational limitation of Maple. If you suppress the display, Maple can compute the number in 16 milliseconds. Then it can be converted to a string of ~20 million base-10 digits in 4 seconds. Then you can easily display any substring of those digits. Or you can skip the string conversion steps and get any subsequence of the digits through pure arithmetic:

#The % signs below are only needed to get an accurate timing;
#they aren't needed to do the computation:
N:= CodeTools:-Usage(value(2%^(2%^26)%+1)):
memory used=16.06MiB, alloc change=16.01MiB, 
cpu time=16.00ms, real time=12.00ms, gc time=0ns

#Make base-10 string:
DN:= CodeTools:-Usage(cat("", N)):
memory used=0.99GiB, alloc change=0 bytes, 
cpu time=3.91s, real time=3.93s, gc time=0ns

L:= length(DN);
                         L := 20201782

#The millionth digit and the 100 following digits:
P:= 10^6..10^6+100:
DN[10^6..10^6+100];
"83772626128948024825321458257987526943544879027080675670992321924911766331675065174616704431088906993"

#Same thing done with pure arithmetic:
CodeTools:-Usage(irem(iquo(N, 10^(L-rhs(P))), 10^(rhs(P)-lhs(P)+1)));
memory used=141.49MiB, alloc change=15.26MiB, 
cpu time=469.00ms, real time=479.00ms, gc time=0ns

 83772626128948024825321458257987526943544879027080675670992321924911766331675065174616704431088906993

When working with large objects, it's crucial to suppress the display, because there's no way to interrupt the display phase once the computation phase has ended. This is true even if no output has been displayed yet because the GUI is still figuring out how to format it.

Let and n be positive integers greater than 1, and let be a nonnegative integer. (Without loss of generality, we can further restrict b < a.) Being someone interested in primes, you should prove the following (it's easy):

Proposition: If p = a*n+b is prime, then a and b are relatively prime (i.e., gcd(a,b) = 1).

The statement "All primes greater than 5 are of the form 6*n+1 or 6*n+5" is an easy consequence of this proposition. Another is that 3, 5, 7 is the only possible sequence of overlapping pairs of twin primes, i.e., 5 is the only prime that belongs to two distinct twin-prime pairs.

It's a fairly deep theorem, due to Dirichlet, that if gcd(a,b) = 1, then there are in fact an infinite number of primes of the form a*n+b, and the sum of their reciprocals diverges.

I think that most Maple commands work in Maple Flow. So, since I don't have Flow, I'll show you how to do it in Maple. Let me know if this works for you:

#Generate a random example:
Data:= LinearAlgebra:-RandomMatrix(9,2):

#Compute least-squares fit line:
F:= Statistics:-LinearFit(a+b*x, Data, x);
          F := -27.2507801140355 + 0.668921237072979 x

plot([F, Data], x= (min..max)(Data[..,1]), style= [line, point]);

 

In your PDF followup question, you propose the following solution:

sq:= anything^{-1/2, 1/2}:
subsindets(
    subsindets(expr, sq, e-> _O*e), 
    And(`*`, satisfies(e-> membertype(sq, {op}(e)) and has(e,exp))),
    e-> eval(simplify(sqrt(e^2)), _O= 1)
);

The problem with this is that it leaves some _Os in the result. That's because _O is added to every sqrt but it only gets removed from the sqrts that are paired with exp. The solution is to wrap the whole thing in eval(..., _O= 1)

sq:= anything^{-1/2, 1/2}:
eval(
    subsindets(
        subsindets(expr, sq, e-> _O*e), 
        And(
            `*`, 
            satisfies(
                e-> andmap(membertype, [sq, specop(exp)], {op}(e))
            )       
        ),
        e-> simplify(sqrt(e^2))
    ),
    _O= 1
);

I also corrected your use of has(e, exp), which is incorrect because it finds exps that are inside the sqrts also.

This trick with the _O is specifically to handle the situation that you discussed at length with Edgardo: `*` with only two operands; the _O becomes the implied coefficient 1. Thus, there's no need to treat two operands and more than two operands as separate cases.

And note that anything can only be matched by exactly one operand; it cannot be matched by the absence of an operand or by a sequence of operands. So, I think the above solution is much easier than using any type that begins with `&*`(...and/or is specific to a specific number of operands.

First 61 62 63 64 65 66 67 Last Page 63 of 394