@Rouben Rostamian  Doesn't the flexibility of the rope have something to do with it? I see no parameter for it.

It seems odd to me that the y-coordinate of the free end is constant. Is that correct?

@vv The command plots:-inequal can handle intersections, unions, and combinations thereof; although it's true that by  using the list/set syntax common to most plotting commands you'll just get an intersection. For a combination of intersections and unions, see the last example on the ?inequal help page.

@vv Yes, certainly it's hopeless for expressions containing user-defined numeric functions that don't control error.

@vv The part that I didn't understand (but now do) was why this doesn't work:

Digits:= 10: #Or just accept the default.
evalf[2](evalf[Digits](p2-p1));

The reason is that Digits = 2 inside the outer evalf.  But if Digits is copied to another variable, that variable won't be changed by the outer evalf.

I'm working on a method to ameliorate the catasthropic cancellation without requiring the user to figure out the number of digits required to do that.

@Carl Love Aha! I figured out how your sigfig works. Each nested level of evalf temporarily changes the value of the global environment variable Digits. I've been trying to figure that out for years, so thanks for the inspiration. With that in mind, here's a reasonable version of sigfig:

sigfig:= proc(e::uneval, d1::posint, d2::posint:= :-Digits)
    evalf[d1](evalf[d2](eval(e)))
end proc:

Usage:

​​​​​​sigfig(p2 - p1, 2);

or if perhaps more precision is required to get the internal computation right:

sigfig(p2 - p1, 2, 2*Digits);

@acer Thank you for teaching me about that difference. Yes, I was mixing up multithreading and hyperthreading.

Do you see any way that Maple could be at fault?

@Glowing Your sigfig trick is very interesting. I can't figure out why it works, but it does work! In particular, I don't understand why using a new variable, Digits0, is any different than using Digits itself.

@Glowing Any ratio greater than 1 shows that there is some hyperthreading. If your O/S is not scheduling optimally or if there's contention over memory caches, etc., it's not Maple's fault. Your ratio is about 4. A factor-of-4 improvement can hardly be characterized by "the benefit of hyperthreading really disappears." Was your computer performing any significant tasks unrelated to Maple at the same time?

@Glowing I don't see any problem in your test run. The hyperthreading performance looks great.

@shkarah I see that someone deleted your Question from today, probably because they thought that it was too similar to this one. The difference that I saw was that in the new Question, you attached code that showed the derivation of Eqs 17-19 from Eq 16. I want to let you know that your derivation was totally correct! The only problem was applying dsolve to Eqs 17-19 in a loop with assign. By simply replacing assign with an equivalent usage of eval, I got dsolve to produce the correct solutions. Here is my worksheet. The only change that I made to your code is that I replaced your last line with my last 4 lines.
 

Download HomotopyPerturbation.mw

@Stretto wrote:

• You say that gcd and mod are for polynomials and that for integers, the most common case, we have to prefix i.

No, not quite. Integer constants are themselves polynomials, so gcd and mod work fine for them also. So, the most-general commands are gcd and mod and they work for the vast majority of cases, both polynomial and integer. Your trouble comes from passing x, which is explicitly a polynomial, to a command which needs to handle (and distinguish) both polynomials and integers. It is only in this case, where x has no assigned value, that you have to clarify the situation by using irem, so that Maple will know that x represents an integer whose value has yet to be assigned.

Now, are you going to say that there's some math book, in number theory or otherwise, that says that x is not a polynomial? So, you're trying to force Maple to view something that is manifestly a polynomial as if it were an integer whose value hasn't been assigned. That's the unusual case, and it's handled by irem (or igcd).

• mod can be generalized to real w.l.o.g.

I think that extending the definition to reals does incur a loss of generality (l.o.g.) because when mod is given a fraction, it finds the multiplicative inverse w.r.t. the modulus. So, for example, 7/2 mod 3 = 2 not 1/2. 

• gcd and mod should be for integer and pgcd and pmod should then be used for polynomials.

For your personal use, it would be very easy to change gcd as you suggest and to create a pgcd ​​​​​​and pmod. Each of these could be done with 1 short line of code.  Changing mod is a little trickier (because it's builtin and a binary infix operator), but it's also possible.

• [I]f the inputs in to these functions are not defined then how the hell can it make any decisions about anything? No one has explained to me why gcd(a,b) returns 1 when a and b are not defined!

As I said in the other thread, if a and b are not defined, then they're explicitly polynomials, as in a = 1*a+0. Does some math book say that a variable with no assigned value, standing alone, is not a polynomial?

• The goal of any programming language(including maple) is not to make life for the programmer hell by having obscure side conditions that can be difficult to detect.

Personally, I've had no trouble learning Maple, and I did it simply from the help pages, a bit of experimentation, and a bit of reading in these online forums. Overall, it makes more consistent sense to me and has fewer obscure roadblocks than any other language that I know: Basic, Pascal, Algol, Fortran, APL, C, Maxima, Kadol, and dBase/FoxPro.

• Maple is wrong in this case. Ideally it would be fixed.

The things that you want "fixed" (mod and gcd) cannot be changed (for general release) because there's nearly 30 years of Maple code that relies on their current definitions. But, if you were writing a package for publication, or just wanted something for personal use, it'd be fairly easy to change in those cases.

I would guess that anyone who's used Maple moderately to extensively has a few ideas about how they would make some things different if they were building it from the ground up. One thing that I might do differently would be to have a way to explictly mark things as polynomials (akin to the way that matrices are explicitly marked) so that things like just plain x wouldn't be mistaken for one. 

• The fact that mod does work for explicit integers FURTHERS the problem because now mod is overload to work for integers in some cases and not others, but you tell me it only works for polynomials?

Integers are polynomials. It works for explicit rational numbers (including integers), algebraic numbers specified in RootOf form, expressions with variables that evaluate to rationals, and matrices of polynomials. It's not true that mod works "for integers in some cases and not others." It works for integers in all cases. An unassigned x is not an integer. That's obvious.

• What precisely is a polynomial to maple? is 3 not a polynomial? is q a polynomial?

Both 3 and q are polynomials to Maple. Polynomials in Maple are pretty much what you'd expect from the mathematical definition: the algebraic closure of the operations `*` , `+`, and `^` with positive integer exponents over the set of objects of type {constant, name}. Expressed as a recursive type in Maple:

TypeTools:-AddType(
    poly,
    Or({constant, name}, poly^posint, specop(poly, {`*`, `+`}))
);

There is of course already a predefined type polynom. The above is just to make the workings explicit.

It makes no sense for you to compare the way that Maple treats x mod 3 when x has no assigned value with any other language that doesn't even allow variables to not be assigned values (i.e., symbolic variables). But if you're willing to make a comparison only among languages that allow symbolic variables, I'd be interested.

@shkarah Post a fully executed worksheet, beginning with restart, that shows your results. And please post it so that it displays in MaplePrimes. Here is my worksheet. The final solutions are mathematically identical to what's shown in the paper; they just need trivial re-arranging to make that obvious. For example, for f[2](x), you just need to distribute the E^2*a/30 over the other terms, reduce the fractions, and substitute alpha for E and eta for x.
 

Download HomotopyPerturbation.mw

@Glowing With no decimal point, 1000 means exactly 1000. But 1000. means that there's some uncertainty about the fractional part (there'll always be uncertainty if it represents a measurement rather than a count). Controlling that uncertainty is the essence of decimal computation. If there's no uncertainty, you should represent 1014.1 as 1014+1/10.

Adding an exact integer to a decimal number is like adding apples and orange juice---the former is counted and the latter is measured.

@Scot Gould A better command to check whether an object has a particular type is type not whattype. Note that I said has, not is. An object can have many types; whattype only returns the most superficial of these. So, you can do

`if`(type(x, float), ..., ...)

This can be abbreviated to 

`if`(x::float, ..., ...)

The command convert(..., int) is meant to convert its input to an integral, not to an integer. It is intended for higher mathematical functions that are defined as or can be represented as integrals. If there's no reasonable way to represent its input with integrals, it just returns its input unchanged, which is what's happening in your case. If you wrote a procedure named `convert/integer`, then convert(..., integer) would work.

It will help you to know that all (software) floats in Maple are represented as x = M*10^E where M and E are integers that can be accessed via (M,E):= op(x). The mantissa M and the exponent E are always integers. The precision is simply the number of digits of M.

@acer There's a help page for frem, but not for fmod. But a look at the code of `evalf/fmod` shows that all that it does (other than some trivialities concerning 0 and infinity) is

(x,y)-> evalf[2*Digits](x-trunc(x/y)*y)