@dharr Thanks for the workaround.  That's very helpful.

A lesson without an example is not very useful.  Here is a concrete example that illustrates how to rotate an object via quaternions.

Download rotation-quaternion.mw

That's rather odd.  If the point are collinear, then the area would be simply zero.  We don't need conditions for that.  Maple's answer would have made sense before the invention of zero.

@one man Okay, now I see what you are doing.  You specify the paths of the contact point independently on the two surfaces, and then roll one surface against the other while requiring the contact point to stay on both paths.  In that case there is one degree of freedom, and all is well.

As to the Latinization of Dragilev, I see that "Dragilev's Method" is written "Метод Драгилева" in Russian, so the Latinized name is "Dragilev".  There is no reason to insert an "h" there.

Edit: The impresario of Ballets Russes, Sergei Diaghilev (Сергей Дягилев), spelled his Latinized name with a "gh", but that was influenced by the Italian orthography, as he lived in Italy, died in Venice, and is buried there.

@one man As far as I know, Dragilev's method applies to systems with one degree of freedom.  The method may be extended to systems with more degrees of freedom but then the resulting equations will be PDEs, not ODEs.

The rolling torus has two degrees of freedom -- it can roll, and it can also spin about its point of contact without rolling.  I can't tell how your worksheet handles that.  Perhaps it makes an assumption on the type of motion that reduces two degrees of freedom to one.

@Carl Love Yes, you are right and I agree with everything that you have written.  I should have been explicit that the roots occur in pairs, and the equal spacing refers to the spacing of the pairs. When I wrote about the roots of f(x) = x^2 + e, I was thinking of the leftmost and rightmost roots in your diagram.

That looks very nice.  Vote up!

Looking through the worksheet, I don't see anything like equations of dynamics.  What does determine the motion of the little torus?

@one man After the transformation x^3 = s, the roots within each cluster are pretty much equally spaced, so it would be difficult to miss any one of them.

There is, however, the following consideration.  Suppose that we are looking for the roots of the function f(x) = x^2 + e.  We see that there are no roots if e > 0 and two roots if e < 0.  But if the size of e is smaller than what the floating point representation can resolve, deciding whether e is positive or negative is ill-posed, and that can lead to lost roots or false roots.

That same thing can happen in the problem that you have posed.  To detect that, I we would do the calculation with some setting of Digits, and then repeat with a larger value of Digits and check if the number of roots changes.

But you are asking for an a priori guarantee for detecting such exceptional cases.  I don't see a good way of doing that.

@one man The number of roots within each cluster is proportional to the length of the interval between the consecutive zeros of the function cos(0.25*surd(s,3)).  These intervals become longer as s grows, and so the cluster sizes become larger.  The first cluster, located in 1800 < s < 2181, is the smallest one, and as we have seen, it contains 122 roots.  In the same way we find that the next cluster corresponds to 15128 < s < 16646 and contains 482 roots.

As the clusters grow in size the farther away we go from the origin, I don't understand what you mean by the "most distant" cluster.

@nm Yes, I agree.   If Maple cannot determine the arbitrary constant, it should return NULL instead of something involving the constant.

@acer And it turns out that the "simplify" is not necessary:

@KIRAN SAJJAN Maple's pdsolve() is not the right tool for solving that problem.  I wrote earlier that pdsolve() is designed for solving evolution equations (initial value problems) while what you need is software for solving a boundary value problem. Maple's numerical PDE solver is not designed for that.

Additionally, I noted earlier that you need to take the outflow condition—equation (3) in the paper—into account.  That's missing in your worksheet.  That condition is not optional; the problem cannot be solved without it.

I am afraid that there is no easy solution to this problem.  If you really need to solve it, you will need to write your own finite elements (FEM) solver.  That will be a nontrivial task and I wouldn't suggest going in that direction unless you are an expert numerical analyst and have extensive experience in writing FEM code.

@KIRAN SAJJAN What you have listed under BCs is not enough.  You will also need to account for the condition expressed as equation (3) in the paper that you have cited.

That, however, won' t help.  There is no symbolic solution to that system of PDEs; I wish there were, but there isn't. If there were, then there wouldn't have been a point in publishing the paper that you have cited; they could have just published the symbolic solution instead.

The best you can hope for is a numerical solution.  But for that you need to specify the numerical values of the parameters, such as lambda, epsilon, Pr, etc.  Even after you do that, Maple's numerical PDE solver won't be able to help you find a solution because it is limited to solving evolution equations, that is, PDEs that depend on time, while yours isn't.  Furthermore, the condition (3) noted above is not a PDE so it needs a special treatment.  You need a specialized solver to solve the PDE numerically, and that would be a somewhat advanced research project.

@KIRAN SAJJAN The original system of PDEs is nonlinear, so we don't expect a symbolic solution which is Maple's main strength.  You wrote that you tried solving the system in Matlab.  That may be an option, but I don't know enough about Matlab's PDE solver to be of help.  You may want to ask in a Matlab forum.

As to

    In matlab also I have tied but for only 0 to 1 I'm getting
    or -1 to 0 only getting.

I can't tell what that means and it's not clear to me whether you fully understand the problem's setting.  The domain of the PDEs is the three-dimensional annular region bounded by the planes z=±h and a < r < b in cylindrical coordinates.

@KIRAN SAJJAN In the paper that you have cited, the author expands the solution into infinite series and thus reduces the PDEs to a set of ODEs (11) and boundary conditions (12).  The graphs shown are those of the solutions of that system of ODEs.

You wrote that you know how to solve ODEs, so you should be able to continue from there.