Carl Love

Carl Love

21698 Reputation

24 Badges

8 years, 353 days
Natick, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity

These are answers submitted by Carl Love

Your ODE can be easily solved symbolically by the methods of first-year calculus: just simple integration. The solution is an elementary function. All those exp functions are just obfuscation because they don't depend on y. All that matters is

ode:= diff(theta(y), y$2) = -(A*sinh(M*y)+B*cosh(M*y))^2;
int(rhs(ode), y);
S:= int(%+C1, y)+C2;
Cs:= solve({eval(S, y=-sigma)=0, eval(S, y=sigma)=1}, {C1,C2});
combine(simplify(eval(S, Cs)));

Then substitute the appropriate expressions for A and B.

If you use dsolve(..., numeric) for this (or any) BVP, it will use a finite difference method automatically.


Another way to avoid "deep" evaluation is to do it in a procedure:

a:= 1:
save a, "A.m":
NewNames:= proc(com::string)
local OldNames:= {anames('user')};
    parse(com, statement);
    {anames('user')} minus OldNames
end proc
b:= 2:
NewNames("read `A.m`;");


You have unbalanced parentheses in ODE11ODE12, and ODE13.

I see this as a bug (or flaw) in the design of DataFrame, but given that that design has been implemented and that that implementation has been released for several years, it wouldn't be easy to fix it and maintain backward compatability.

@vv wrote: Probably such numerical labels should not be accepted.

I disagree. Here is how I would've designed it: Allow anything to be labels, but if labels are used then those labels rather than row and column numbers should be the only way to index. While it would be difficult to change DataFrame itself in this manner (because of the backward compatibility), it would be very easy to write a wrapper for it that implemented this idea. Since a DataFrame is itself just a wrapper for a matrix, it'd still be trivial for a programmer to index the underlying structure by rows and columns.

I worked on large databases with significant user interfaces (point-of-sale systems, inventory, etc.) for several years in the late '80s and early '90s. I learned the following the hard way: There is no good reason that the end user should ever see, let alone refer to, the actual record numbers (or any other field that is used as the primary key). Indeed, if you allow them to see those keys, several problems will gradually arise:

  1. Gradually (and perhaps at first subconsciously) they begin to ascribe significance or meaning to certain numbers or keys or ranges thereof. That's just the way the human mind works.
  2. Once that happens, occasions will arise where they will want to change some of those keys. Real-life example: A customer's status changes from retail to contractor. The user wants to change the customer's number in whatever way indicates to them that that's a contractor, even though I've explained to them that it'd be much less work for me (and cost for them) to just add a 1-bit field for contractor status.
  3. Doing that can be a major project because those numbers (or primary keys) need to be linked to many other files in the overall database.

The way around this is to use a secondary key, which will seem to the end user to be the primary key, yet they can do however they please with it. But the secondary key only exists in that one file (and its index file). All linkages to other files are via the primary key. In the case of DataFrame, the labels are those secondary keys.

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:

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:

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)                 
                               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 (x - y(x))     
                             - ------------- = 1
                               (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})})->
            C, And(name, indexed)^posint,
            e-> `*`(
                'convert(op([1,0], e), `local`)[(op@op)(1,e)]' 
                $ op(2,e)
                         [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


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:

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();

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. 

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