Maple Questions and Posts

These are Posts and Questions associated with the product, Maple

given

ode:=2*x^(1/2)*diff(y(x),x)-y(x) = -sin(x^(1/2))-cos(x^(1/2)); 
ic:=y(infinity) = y__0; 
sol:=dsolve([ode,ic]);

It gives  

This solution satisfies the ode itself. Now cos(sqrt(x)) when x=infinity is  -1..+1

But IC says y(infinity)=y0  so odetest do not verify the IC and gives this

odetest(sol,[ode,ic]);

I think dsolve should not have returned a solution at all. 

What do the experts here think of this result?

Maple 2026.1 on windows 10

I'm trying to rewrite my initialization routines a bit, and want to find out if there are any other procedures with a certain name that should be called.
However there seems to be a problem with the scope. I'm trying to explain that with a little program.

The question is - why is InitSpecific not found by the InitCommon procedure here?

Module1 := module()

 

``

Module2 := module()

Module1:-InitCommon()

"false"

(1)

Module2:-InitSpecific()

"InitSpecific"

(2)

NULL

Download InitCommon.mw

The uploaded worksheet references two youtube videos.

The first one displays the animation of a simple device rotating about an axis tilted at a small angle from the device's principal axis having an intermediate moment of inertia.

The animation and accompanying verbal description demonstrate the Dzhanibekov effect.
The second video contains the first video's narrator's equations which produce the values used in creating the animation.

The uploaded worksheet contains my failed attempt to reproduce these values.

Please suggest the Maple 2020 compatible statements which correctly produce these values.

Dzhanibekov_effect.mw

I am new to evala (and the math behind).

On ?evala,Sqrfree at the first bullet point I was wondering if there isn't a "u" missing here

Can someone confirm?

I asked an LLM to provide an expansion of the MacDonald function of arbitrary order (a modified Bessel function of the second kind with purely imaginary order and positive argument), K(I*y,r), as a weighted sum of MacDonald functions of integer order. It came back with

         K(I*y,z)=2*sinh(Pi*y)/Pi* [K(0,r)/2*y+sum( (-1)^n*y*BesselK(0,r)/(y^2+n^2),n=1..infinity)]

(see below for more readable text)

I evaluated the LHS and RHS using Maple 2026 for various choices of y and r and found numerical agreement using both "sum" and "Sum".  I was very pleased until I realized that the RHS isn't a convergent series!

Can anyone explain to me how Maple pulls this off! 

(I asked Maplesoft Tech Support but they said it is above their pay grade... I suspect that Maple is using Borel summability to evaluate the RHS but I haven't been able to verify that)

I apologize, but I can't see how to attach a .mw file, so I've cut and pasted the code below

restart;

T := R(xi)*R(xi) + lambda;

u := A[0] + A[1]*R(xi) + B[1]/R(xi);

d[1] := A[1]*T - B[1]*T/R(xi)^2;

d[2] := 2*A[1]*R(xi)*T - 2*B[1]*T/R(xi) + 2*B[1]*(R(xi)^2 + lambda)*T/R(xi)^3;

expand(((-alpha^2*b^2 + a^2)*alpha^2)/(2*beta)*d[2] + (omega + alpha^2*(alpha^2*l^2 + k^2)/2 - a*C[1]/(-alpha^2*b^2 + a^2))*u[0]/(beta - 2*beta*a^2/(-alpha^2*b^2 + a^2)) + u[0]*u[0]*u[0]);

value(%);

simplify(%);

collect(%, R(xi));


      /      6  4    4      2\      3
 A[1] \-alpha  b  + a  alpha / R(xi) 
 ------------------------------------
             /     2  2    2\        
        beta \alpha  b  + a /        

                  /      6  4    4      2\         
      A[1] lambda \-alpha  b  + a  alpha / R(xi)   
    + ------------------------------------------ + 
                     /     2  2    2\              
                beta \alpha  b  + a /              

                         /                               /     
             1           |/                    B[1] \    |     
   --------------------- ||A[0] + A[1] R(xi) + -----|[0] |beta 
        /     2  2    2\ \\                    R(xi)/    \     
   beta \alpha  b  + a /                                       

                                                  2
   /     2  2    2\ /                    B[1] \    
   \alpha  b  + a / |A[0] + A[1] R(xi) + -----|[0] 
                    \                    R(xi)/    

      1  2  2      6   1 /  2  2    2  2\      4
    + - b  l  alpha  + - \-a  l  + b  k / alpha 
      2                2                        

                                                       \\
      /  1  2  2    2      \      2    2               ||
    + |- - a  k  + b  omega| alpha  - a  omega + a C[1]||
      \  2                 /                           //

                  /      6  4    4      2\
      B[1] lambda \-alpha  b  + a  alpha /
    + ------------------------------------
               /     2  2    2\           
          beta \alpha  b  + a / R(xi)     

            6  4       2         4      2       2     
      -alpha  b  lambda  B[1] + a  alpha  lambda  B[1]
    + ------------------------------------------------
                     /     2  2    2\      3          
                beta \alpha  b  + a / R(xi)           

WHen I open many worksheets at same time, say 10. The new UI do not stack them all (i.e. the tab at the top), forcing one to use the small arrow to navigate to each worksheet.

Is there a way to tell the UI to show all tabs (may be double rows and 3 rows as needed) to make it easier to jump from one worksheet to the other?

I do not know if this is new feature in the new ribbon UI or not. 

Here is screen show where I have 10 worksheets open

There is also a pull down menu, but it only shows 8 worksheets and one can have more open but they do not show. So have to scroll down looking for the rest. Even that does not work well. many times when I try to scroll down, the window closes. It will not give me time to move the mouse to the scroll bar to move it before it closes.

Both of these solutions are not good. Having to use the arrow key to look and navigate for a different worksheet is bad UI design.

How to see all tabs for all open worksheet in same UI?  If the tabs do not fit on one row, why not make second row? If two rows do not fit, make 3rd row. This should be an option for the user. But I did not see one so far. But will keep looking.

I find tabs where all worksheet show much better design that this UI design.   

I only use worksheet and not document mode. Windows 10.

To give you idea what I mean, These are examples found on the net of stacked tabs

 

 

Where in Maple, each tab above will have the name of the worksheet open. Font can be small, is OK.

Is it possible to have this in the new UI for open worksheets?

One option I might try to make my worksheets names much shorter. May be then they will fit all in same window.

The below problem has already occured several times to me. In all such instances Maple did not realise that extracting a factor from a square root is the key for further simplification. Doing this by hand is obvious and often easy when extracted factors are positive.  

Did I overlook something? Are there other ways avoid disassembling an expression with the op command?
Should simplify or other commands be improved to adress such problems?

restart

How to transform the left-hand side by commands that it matches the right-hand side

sqrt(x__0+1)*sqrt(-2*beta^2*x__0-2*beta^2+4)*sqrt(-(x__0+1)*(beta^2-1))/((beta^2*x__0+beta^2-2)*(beta^2*x__0+beta^2-x__0-1)) = 2/(sqrt(-beta^2+1)*sqrt(-2*beta^2*x__0-2*beta^2+4))

(x__0+1)^(1/2)*(-2*beta^2*x__0-2*beta^2+4)^(1/2)*(-(x__0+1)*(beta^2-1))^(1/2)/((beta^2*x__0+beta^2-2)*(beta^2*x__0+beta^2-x__0-1)) = 2/((-beta^2+1)^(1/2)*(-2*beta^2*x__0-2*beta^2+4)^(1/2))

(1)

assumptions := 0 < x__0 and x__0 < 1, 0 < beta and beta < 1

0 < x__0 and x__0 < 1, 0 < beta and beta < 1

(2)

`assuming`([simplify(lhs((x__0+1)^(1/2)*(-2*beta^2*x__0-2*beta^2+4)^(1/2)*(-(x__0+1)*(beta^2-1))^(1/2)/((beta^2*x__0+beta^2-2)*(beta^2*x__0+beta^2-x__0-1)) = 2/((-beta^2+1)^(1/2)*(-2*beta^2*x__0-2*beta^2+4)^(1/2))))], [assumptions])

-(4+(-2*x__0-2)*beta^2)^(1/2)/((-beta^2+1)^(1/2)*(-2+(x__0+1)*beta^2))

(3)

I have tried the usual simplify and combine commands to remove the square root from the numerator.
Extracting a factor for -2 from the square root would probably make further simplification possible but there is no simple command to do so.

Factor_ := -2

-2

(4)

old := simplify([op(denom(-(4+(-2*x__0-2)*beta^2)^(1/2)/((-beta^2+1)^(1/2)*(-2+(x__0+1)*beta^2))))])

[(-beta^2+1)^(1/2), -2+(x__0+1)*beta^2]

(5)

new := old; new[1] := old[1]/Factor_; new[2] := old[2]*Factor_

[-(1/2)*(-beta^2+1)^(1/2), 4-2*(x__0+1)*beta^2]

(6)

subs(1/old[1] = 1/new[1], 1/old[2] = 1/new[2], -(4+(-2*x__0-2)*beta^2)^(1/2)/((-beta^2+1)^(1/2)*(-2+(x__0+1)*beta^2)))

2*(4+(-2*x__0-2)*beta^2)^(1/2)/((-beta^2+1)^(1/2)*(4-2*(x__0+1)*beta^2))

(7)

expand(simplify(2*(4+(-2*x__0-2)*beta^2)^(1/2)/((-beta^2+1)^(1/2)*(4-2*(x__0+1)*beta^2))))

2/((-beta^2+1)^(1/2)*(-2*beta^2*x__0-2*beta^2+4)^(1/2))

(8)

2/((-beta^2+1)^(1/2)*(-2*beta^2*x__0-2*beta^2+4)^(1/2)) = rhs((x__0+1)^(1/2)*(-2*beta^2*x__0-2*beta^2+4)^(1/2)*(-(x__0+1)*(beta^2-1))^(1/2)/((beta^2*x__0+beta^2-2)*(beta^2*x__0+beta^2-x__0-1)) = 2/((-beta^2+1)^(1/2)*(-2*beta^2*x__0-2*beta^2+4)^(1/2)))

2/((-beta^2+1)^(1/2)*(-2*beta^2*x__0-2*beta^2+4)^(1/2)) = 2/((-beta^2+1)^(1/2)*(-2*beta^2*x__0-2*beta^2+4)^(1/2))

(9)

is(2/((-beta^2+1)^(1/2)*(-2*beta^2*x__0-2*beta^2+4)^(1/2)) = 2/((-beta^2+1)^(1/2)*(-2*beta^2*x__0-2*beta^2+4)^(1/2)))

true

(10)

Second approach after "discovering" that content works also on square roots

[op(-(4+(-2*x__0-2)*beta^2)^(1/2)/((-beta^2+1)^(1/2)*(-2+(x__0+1)*beta^2)))]

[-1, 1/(-beta^2+1)^(1/2), (4+(-2*x__0-2)*beta^2)^(1/2), 1/(-2+(x__0+1)*beta^2)]

(11)

mul(`~`[`*`](`~`[content]([-1, 1/(-beta^2+1)^(1/2), (4+(-2*x__0-2)*beta^2)^(1/2), 1/(-2+(x__0+1)*beta^2)]), `~`[primpart]([-1, 1/(-beta^2+1)^(1/2), (4+(-2*x__0-2)*beta^2)^(1/2), 1/(-2+(x__0+1)*beta^2)])))

-2^(1/2)*(-beta^2*x__0-beta^2+2)^(1/2)/((-beta^2+1)^(1/2)*(beta^2*x__0+beta^2-2))

(12)

simplify(-2^(1/2)*(-beta^2*x__0-beta^2+2)^(1/2)/((-beta^2+1)^(1/2)*(beta^2*x__0+beta^2-2))) = rhs((x__0+1)^(1/2)*(-2*beta^2*x__0-2*beta^2+4)^(1/2)*(-(x__0+1)*(beta^2-1))^(1/2)/((beta^2*x__0+beta^2-2)*(beta^2*x__0+beta^2-x__0-1)) = 2/((-beta^2+1)^(1/2)*(-2*beta^2*x__0-2*beta^2+4)^(1/2)))

2^(1/2)/((2+(-x__0-1)*beta^2)^(1/2)*(-beta^2+1)^(1/2)) = 2/((-beta^2+1)^(1/2)*(-2*beta^2*x__0-2*beta^2+4)^(1/2))

(13)

is(%)

true

(14)

NULL

Context: The left-hand side in an integrand which was produced by a change of variables in a elliptic integral. Maple simplifies only halfway which makes validation of the result of the variable change difficult.  

NULL

Related functional programming question: Is a onliner `...`(-(4+(-2*x__0-2)*beta^2)^(1/2)/((-beta^2+1)^(1/2)*(-2+(x__0+1)*beta^2)))from the above content-primpart construct possible?NULL

Download Simplify_radical_02.mw

Help me rewrite the code to create visible Bar chart of different colors. I can't figure out why this code is not giving me a visible bar graph

restart; with(Statistics); with(plots); Data := [45, 38, 51, 67, 74, 91]; P := BarChart(Data, tickmarks = [[1 = "Chemical Vector Control", 2 = "Resistant Cultivars", 3 = "Roguing & Sanitation", 4 = "u1+u2", 5 = "u1+u3", 6 = "Integrated"], default], width = .75); T := textplot([[1, 48, "45%"], [2, 41, "38%"], [3, 54, "51%"], [4, 70, "67%"], [5, 77, "74%"], [6, 94, "91%"]], font = ["TIMES", "BOLD", 12]); display([P, T], title = "Figure 20: Comparative Effectiveness of Optimal Control Strategies", labels = ["Control Strategies", "Reduction in Coinfection Burden (%)"], labelfont = ["TIMES", "BOLD", 14], titlefont = ["TIMES", "BOLD", 16], axes = boxed, gridlines = true, view = [.5 .. 6.5, 0 .. 100], size = [1000, 650])

 
 

NULL

Download Bargraph.mw

restart;
with(plottools);
with(plots);
a := 1;
b := 1;
c := 1;
k := 1;
l := 1;
omega := 1;
A[2] = 2;
alpha := 2;
beta := 1;
kappa := 0.5;
C[1] := 1;
lambda := -1;

omega := (-alpha^6*b^4*lambda + 2*alpha^6*b^2*l^2 - 2*a^2*alpha^4*l^2 + 2*alpha^4*b^2*k^2 + a^4*alpha^2*lambda - 2*a^2*alpha^2*k^2 + 4*a*C[1])/(-4*alpha^2*b^2 + 4*a^2);

a[0] := 0;

a[1] := sqrt(-(-alpha^2*b^2 + a^2)/(4*beta))*alpha;

b[1] := sqrt(-(alpha^2*b^2*lambda*sigma - a^2*lambda*sigma)/(4*beta))*alpha;

sigma := A[1]*A[1] - A[2]*A[2];

T := A[1]*sinh(xi*sqrt(-lambda)) + A[2]*cosh(xi*sqrt(-lambda)) + mu/lambda;

t := diff(T, xi);

S := t/T;

R := 1/T;

mu := 0;

A[1] := 0;

y := 0;

xi := k*x^kappa/kappa + l*y^kappa/kappa - omega*t^kappa/kappa;

  Error, recursive assignment

Dear all,

I'm reporting what seems to me as a bug in the SMTLIB library in maple. 

    |\^/|     Maple 2026 (X86 64 LINUX)
._|\|   |/|_. Copyright (c) Maplesoft, a division of Waterloo Maple Inc. 2026
 \  MAPLE  /  All rights reserved. Maple is a trademark of
 <____ ____>  Waterloo Maple Inc.
      |       Type ? for help.
> SMTLIB:-Satisfiable({x^2=2,y^2=2,x<y});
                                     true

> SMTLIB:-Satisfiable({x^2=2,y^2=2,y<x});
                                     false

> SMTLIB:-Satisfiable({x^2=2,a^2=2,a<x});
                                     true

The Satisfiable command do not output the correct decision on two formulas of equivalent realization by switching x<y (output SAT) to y<x (output UNSAT). I suspect this is because some alphabetical order depandance in the variables as for a<y we get SAT again.

I tried to feed Z3 with the code given by ToString on the problematic formula and I get two different outputs :

  • on the Z3 version 4.8.12 from the ubuntu repository (apt install) I also get the wrong UNSAT output;
  • one the Z3 version 4.17.0 build from the official github repository I finally get the correct SAT output.

Thus, I suspect a version problem in SMTLIB that do not take in account the last updates made in SMT solvers (Z3?).

Many thanks for considering my problem!

A little while ago, I created a video, Engaging and Enlightening Students with Maple Visualizations, that showed a sample of Maple visualizations that would be helpful in teaching math. Doing that allowed me to get reacquainted with some of Maple's plotting features that I hadn't used for a while. As a result, I made a second instructional video for my Maple tips series, Animating a Polyhedron in Maple

I chose this topic because I thought it would show several features in Maple that might not be known to all users. I list them below and encourage you to try them out.

  • The plots:-polyhedraplot command allows you to create a 3-D plot of a polyhedron, including one of 138 polyhedra that Maple knows about.

  • The list of named polyhedra available can be obtained by calling the plots:-polyhedra_supported command.

  • The viewpoint option, which allows you to create an animation by varying the viewpoint through a 3D plot, can be used to rotate the polyhedron.

  • Finally, the Export feature allows you to save the plot animation as an animated GIF.

 

restart;
solve({-alpha^4*b^2*lambda*mu*b[1] + a^2*alpha^2*lambda*mu*b[1] + 6*beta*lambda^2*sigma*a[0]*a[1]^2 + 6*beta*mu^2*a[0]*a[1]^2 - 6*beta*lambda*a[0]*b[1]^2 = 0, -alpha^4*b^2*lambda^2*mu*sigma - alpha^4*b^2*mu^3 + a^2*alpha^2*lambda^2*mu*sigma + a^2*alpha^2*mu^3 - 4*beta*lambda^2*sigma*a[0]*b[1] - 4*beta*lambda*mu*b[1]^2 - 4*beta*mu^2*a[0]*b[1] = 0, -alpha^4*b^2*lambda^2*sigma - alpha^4*b^2*mu^2 + a^2*alpha^2*lambda^2*sigma + a^2*alpha^2*mu^2 + beta*lambda^2*sigma*a[1]^2 + beta*mu^2*a[1]^2 - 3*beta*lambda*b[1]^2 = 0, -alpha^4*b^2*lambda^2*sigma - alpha^4*b^2*mu^2 + a^2*alpha^2*lambda^2*sigma + a^2*alpha^2*mu^2 + 3*beta*lambda^2*sigma*a[1]^2 + 3*beta*mu^2*a[1]^2 - beta*lambda*b[1]^2 = 0, -alpha^6*b^4*lambda^2*mu*b[1] + alpha^6*b^2*l^2*lambda^2*sigma*a[0] + alpha^6*b^2*l^2*mu^2*a[0] - a^2*alpha^4*l^2*lambda^2*sigma*a[0] + alpha^4*b^2*k^2*lambda^2*sigma*a[0] - a^2*alpha^4*l^2*mu^2*a[0] + alpha^4*b^2*k^2*mu^2*a[0] + 2*alpha^2*b^2*beta*lambda^2*sigma*a[0]^3 + a^4*alpha^2*lambda^2*mu*b[1] - a^2*alpha^2*k^2*lambda^2*sigma*a[0] - 6*alpha^2*b^2*beta*lambda^2*a[0]*b[1]^2 + 2*alpha^2*b^2*beta*mu^2*a[0]^3 - a^2*alpha^2*k^2*mu^2*a[0] + 2*a^2*beta*lambda^2*sigma*a[0]^3 + 2*alpha^2*b^2*lambda^2*omega*sigma*a[0] - 6*a^2*beta*lambda^2*a[0]*b[1]^2 + 2*a^2*beta*mu^2*a[0]^3 + 2*alpha^2*b^2*mu^2*omega*a[0] - 2*a^2*lambda^2*omega*sigma*a[0] - 2*a^2*mu^2*omega*a[0] + 2*a*lambda^2*sigma*C[1]*a[0] + 2*a*mu^2*C[1]*a[0] = 0, -2*alpha^6*b^4*lambda^3*sigma - 2*alpha^6*b^4*lambda*mu^2 + alpha^6*b^2*l^2*lambda^2*sigma + alpha^6*b^2*l^2*mu^2 - a^2*alpha^4*l^2*lambda^2*sigma + alpha^4*b^2*k^2*lambda^2*sigma + 2*a^4*alpha^2*lambda^3*sigma - a^2*alpha^4*l^2*mu^2 + alpha^4*b^2*k^2*mu^2 + 6*alpha^2*b^2*beta*lambda^2*sigma*a[0]^2 + 2*a^4*alpha^2*lambda*mu^2 - a^2*alpha^2*k^2*lambda^2*sigma - 6*alpha^2*b^2*beta*lambda^2*b[1]^2 + 6*alpha^2*b^2*beta*mu^2*a[0]^2 - a^2*alpha^2*k^2*mu^2 + 6*a^2*beta*lambda^2*sigma*a[0]^2 + 2*alpha^2*b^2*lambda^2*omega*sigma - 6*a^2*beta*lambda^2*b[1]^2 + 6*a^2*beta*mu^2*a[0]^2 + 2*alpha^2*b^2*mu^2*omega - 2*a^2*lambda^2*omega*sigma - 2*a^2*mu^2*omega + 2*a*lambda^2*sigma*C[1] + 2*a*mu^2*C[1] = 0, -alpha^6*b^4*lambda^3*sigma + alpha^6*b^4*lambda*mu^2 + alpha^6*b^2*l^2*lambda^2*sigma + alpha^6*b^2*l^2*mu^2 - a^2*alpha^4*l^2*lambda^2*sigma + alpha^4*b^2*k^2*lambda^2*sigma + a^4*alpha^2*lambda^3*sigma - a^2*alpha^4*l^2*mu^2 + alpha^4*b^2*k^2*mu^2 + 6*alpha^2*b^2*beta*lambda^2*sigma*a[0]^2 - a^4*alpha^2*lambda*mu^2 - a^2*alpha^2*k^2*lambda^2*sigma - 2*alpha^2*b^2*beta*lambda^2*b[1]^2 + 12*alpha^2*b^2*beta*lambda*mu*a[0]*b[1] + 6*alpha^2*b^2*beta*mu^2*a[0]^2 - a^2*alpha^2*k^2*mu^2 + 6*a^2*beta*lambda^2*sigma*a[0]^2 + 2*alpha^2*b^2*lambda^2*omega*sigma - 2*a^2*beta*lambda^2*b[1]^2 + 12*a^2*beta*lambda*mu*a[0]*b[1] + 6*a^2*beta*mu^2*a[0]^2 + 2*alpha^2*b^2*mu^2*omega - 2*a^2*lambda^2*omega*sigma - 2*a^2*mu^2*omega + 2*a*lambda^2*sigma*C[1] + 2*a*mu^2*C[1] = 0}, {omega, a[0], a[1], b[1]});
 

Recently @salim-barzani asked a question about a paper that involved analysing the different types of roots of polynomials. The appendix in that paper gave the example of the roots of x^4+x^2*e[2]+x*e[1]+e[0]using the analysis in Lu et al, "A complete discrimination system for polynomials", Science in China (Ser. E), 39 (1996) 628-646. The analysis uses the discriminant sequence and extensions. Maple provides this through RegularChains:-ParametricSystemTools:-DiscrminantSequence. For example for this polynomial we find there is a real root of multiplicity 2 and a complex conjugate pair when D__2*D__3 < 0 and D__4 = 0 where the D__i are the ith entries in the discriminant sequence [1, -e[2], -2*e[2]^3+8*e[0]*e[2]-9*e[1]^2, 16*e[0]*e[2]^4-4*e[1]^2*e[2]^3-128*e[0]^2*e[2]^2+144*e[0]*e[1]^2*e[2]-27*e[1]^4+256*e[0]^3].

 

The problem with these conditions is that they are in terms of the D__i and not directly in terms of the e__i parameters. One can derive these conditions and then solve them to find the conditions on the parameters, but Maple has various routines in the RegularChains, RootFinding:-Parametric and SolveTools packages that directly find conditions on parameters to find when there are specified numbers of real or complex roots for polynomial systems. So this post is my attempt to use these tools to find the conditions on the parameters of the above polynomial that give various types of roots. One immediate difficulty is that generally these routines count distinct roots irrespective of multiplicity, and so some indirect analysis is required. There are several different types of commands and analyses that could be used, and my choices here are more to do with my learning experience than an optimum analysis.

 

The first conclusion is that it is possible, although RealComprehensiveTriangularize did not work as I expected when asking for zero real roots (see cases (a) and (b)) (bug?). Assuming it had worked, RealComprehensiveTriangularize could cover all the cases here, though that will not be true for higher-degree polynomials with more parameters. There doesn't seem to be an obvious systematic way of doing this analysis, which is a downside. Another downside is the large number of subcase conditions, which look as if they could be combined into fewer subcases. CellDecomposition works well for cases without multiplicity.


Main worksheet [not all is displayed below]:

Download RootAnalysis4.mw

restart

with(RegularChains); with(ParametricSystemTools); with(RootFinding:-Parametric)

Consider the following polynomial in x with the three real parameters e[0], e[1], e[2]. We would like to know the conditions on these parameters that lead to different numbers of real and complex-conjugate pairs of roots of different multiplicities.

p := x^4+x^2*e[2]+x*e[1]+e[0]

x^4+x^2*e[2]+x*e[1]+e[0]

Consider first how many cases there are. We can set this up as a combinatorial problem in the combstruct package.

sys := {C = Atom, R = Atom, realrts = Set(multiplereal), rts = Prod(realrts, complexrts), complexpr = Prod(C, C), complexrts = Set(multiplecomplex), multiplecomplex = Sequence(complexpr, card > 0), multiplereal = Sequence(R, card > 0)}

Draw := proc (q) options operator, arrow; eval(q, {Epsilon = NULL, Prod = `[]`, Set = (proc () options operator, arrow; args end proc), Sequence = `*`}) end proc

For a degree 4 polynomial there are 9 different cases to consider. Here [C, C] means a (non-real) complex-conjugate pair of roots and R means a real root; the exponents indicate the multiplicities.

all := combstruct:-allstructs([rts, sys], size = degree(p, x)); nops(%); `~`[Draw](all)

9

[[[C, C]^2], [[C, C], [C, C]], [R^2, [C, C]], [R^4], [R, R, [C, C]], [R^2, R^2], [R^3, R], [R, R, R^2], [R, R, R, R]]

These are (in order)
(a) A duplicate pair of complex-conjugate roots

(b) Two distinct pairs of complex-conjugate roots

(c) A real root of multiplicity 2 and a pair of complex-conjugate roots

(d) A real root of multiplicity 4

(e) Two distinct real roots of multiplicity 1 and a complex-conjugate pair

(f) Two distinct real roots each of multiplicity 2

(g) A real root of multiplicity 3 and a real root of multiplicity 1

(h) Three distinct real roots, of multiplicities 2, 1, and 1

(i) Four distinct real roots of multiplicity 1

Declare the variables first and parameters last. np is the number of parameters. Use the suggested order.

vp := SuggestVariableOrder([p = 0], [x]); R := PolynomialRing(vp); np := nops(`minus`({vp[]}, {x}))

[x, e[0], e[1], e[2]]

polynomial_ring

3

Define derivative polynomials. p2 = 0 when there is a root of multiplicity 2 or more; p3 = 0 when there is a root of multiplicity 3 or more and p4 = 0when there is a root of multiplicity 4.

p2 := diff(p, x); p3 := diff(p2, x); p4 := diff(p3, x)

4*x^3+2*x*e[2]+e[1]

12*x^2+2*e[2]

24*x

Discriminant is zero if and only if there are repeated roots.

Delta := discrim(p, x)

16*e[0]*e[2]^4-4*e[1]^2*e[2]^3-128*e[0]^2*e[2]^2+144*e[0]*e[1]^2*e[2]-27*e[1]^4+256*e[0]^3

(d) A real root of multiplicity 4

 

This is perhaps the simplest case and can be done using RealComprehensiveTriangularize. Specifying np = 3 means use the last 3 variables in the PolynomialRing as parameters. The argument 1 means we want the cases where there is one distinct real root. We specify that all of p, p2, p3, p4 are zero so that the 1 real root is the common root of these polynomials, i.e., is a root on multiplicity 4. (I find the cadcell output a little easier to use, but it is not critical.)

rct := RealComprehensiveTriangularize({p = 0, p2 = 0, p3 = 0, p4 = 0}, np, R, 1, output = cadcell); Display(rct, R)

PiecewiseTools:-Is, "Wrong kind of parameters in piecewise"

This means that the conditions on the parameters to get a single real root of multiplicity 4 are

conds := Info(rct[2][1][1], R)

[e[2] = 0, e[1] = 0, e[0] = 0]

and that the polynomial to solve to find this root under these conditions is

poly := Info(rct[1][][2], R)

[x = 0]

which we can check:

eval(p, conds); solve(%, x)

x^4

0, 0, 0, 0

(g) A real root of multiplicity 3 and a real root of multiplicity 1

   

(f) Two distinct real roots each of multiplicity 2

   

(i) Four distinct real roots of multiplicity 1

   

(h) Three distinct real roots, of multiplicities 2, 1, and 1

   

(e) Two distinct real roots of multiplicity 1 and a complex-conjugate pair

   

(c) A real root of multiplicity 2 and a pair of complex-conjugate roots

   

(b) Two distinct pairs of complex-conjugate roots

   

(a) A duplicate pair of complex-conjugate roots

 

Here we want multiplicity 2 but no real roots. The discriminant is expected to be zero, but in most cases is not.

rct := RealComprehensiveTriangularize({p = 0, p2 = 0, p3 <> 0, p4 <> 0}, np, R, 0, output = cadcell); nops(rct[2]); Display(rct[2][1 .. 4], R)

40

[[PIECEWISE([e[0] < RootOf(256*_Z^3-128*e[2]^2*_Z^2+(16*e[2]^4+144*e[1]^2*e[2])*_Z-4*e[1]^2*e[2]^3-27*e[1]^4, index = real[1]), ``], [e[1] < -(2/9)*(-6*e[2]^3)^(1/2), ``], [e[2] < 0, ``]), []], [PIECEWISE([RootOf(256*_Z^3-128*e[2]^2*_Z^2+(16*e[2]^4+144*e[1]^2*e[2])*_Z-4*e[1]^2*e[2]^3-27*e[1]^4, index = real[1]) < e[0], ``], [e[1] < -(2/9)*(-6*e[2]^3)^(1/2), ``], [e[2] < 0, ``]), []], [PIECEWISE([e[0] < -(1/12)*e[2]^2, ``], [e[1] = -(2/9)*(-6*e[2]^3)^(1/2), ``], [e[2] < 0, ``]), []], [PIECEWISE([e[0] = -(1/12)*e[2]^2, ``], [e[1] = -(2/9)*(-6*e[2]^3)^(1/2), ``], [e[2] < 0, ``]), []]]

Consider one of the cases with an equality for e[0], suggesting a zero discriminant.
Find a sample point satisfying the conditions, and see what the roots are like

j := 37; cell := rct[2][j][1]; conds := Info(cell, R); pts := Info(SamplePoints(cell, R), R)[1]; `~`[is](eval(conds, pts))

37

cad_cell

[0 < e[2], e[1] = 0, e[0] = (1/4)*e[2]^2]

[e[0] = 1/16, e[1] = 0, e[2] = 1/2]

[true, true, true]

We find two complex roots of multiplicity 2, which are complex conjugates, as expected

eval(p, pts); solve(%, x)

x^4+(1/2)*x^2+1/16

(1/2)*I, -(1/2)*I, (1/2)*I, -(1/2)*I

Check that discriminant is zero.

eval(Delta, pts)

0

But, for example, the first cell does not have the discriminant zero, and does not give a correct result.

j := 1; cell := rct[2][j][1]; conds := Info(cell, R); pts := Info(SamplePoints(cell, R), R)[1]; `~`[is](eval(conds, pts))

1

cad_cell

[e[2] < 0, e[1] < -(2/9)*(-6*e[2]^3)^(1/2), e[0] < RootOf(256*_Z^3-128*e[2]^2*_Z^2+(16*e[2]^4+144*e[1]^2*e[2])*_Z-4*e[1]^2*e[2]^3-27*e[1]^4, index = real[1])]

[e[0] = 1/2, e[1] = -3/2, e[2] = -1/2]

[true, true, true]

We find a complex conjugate pair and two distinct real roots, which is not expected for this case.

eval(p, pts); fsolve(%, x, complex)

x^4-(1/2)*x^2-(3/2)*x+1/2

-.7473459056-.9001675303*I, -.7473459056+.9001675303*I, .3077440660, 1.186947745

The discriminant is indeed nonzero.

eval(Delta, pts)

-3073/16

We did not yet find the conditions for this case. We can ask for the conditions for different numbers of complex roots (complex in this context includes real).

cmplx := ComplexRootClassification([p], np, R)

[[constructible_set, 1], [constructible_set, 2], [constructible_set, 3], [constructible_set, 4]]

We are interested in the conditions for exactly two distinct roots, which is found from the second constructible_set. There are two subcases.

Display(cmplx[2][1], R)

PIECEWISE([12*e[0]+e[2]^2 = 0, ``], [27*e[1]^2+8*e[2]^3 = 0, ``], [e[2] <> 0, ``]), PIECEWISE([4*e[0]-e[2]^2 = 0, ``], [e[1] = 0, ``], [e[2] <> 0, ``])

Consider the second case

case2 := Info(cmplx[2][1], R)[2]

[[4*e[0]-e[2]^2, e[1]], [e[2]]]

solve(case2[1], {e[0], e[1]}); q1 := eval(p, %); rts2 := [solve(%, x)]

{e[0] = (1/4)*e[2]^2, e[1] = 0}

x^4+x^2*e[2]+(1/4)*e[2]^2

[(1/2)*(-2*e[2])^(1/2), -(1/2)*(-2*e[2])^(1/2), (1/2)*(-2*e[2])^(1/2), -(1/2)*(-2*e[2])^(1/2)]

We saw this before when e[2]<0 as case (f) with real roots of multiplicity 2. Now, for e[2]>0 we indeed have two duplicate pairs of complex conjugate roots

`assuming`([simplify(rts2)], [e[2] > 0])

[((1/2)*I)*2^(1/2)*e[2]^(1/2), -((1/2)*I)*2^(1/2)*e[2]^(1/2), ((1/2)*I)*2^(1/2)*e[2]^(1/2), -((1/2)*I)*2^(1/2)*e[2]^(1/2)]

Consider the first case

case1 := Info(cmplx[2][1], R)[1]

[[12*e[0]+e[2]^2, 27*e[1]^2+8*e[2]^3], [e[2]]]

ans1 := {solve(case1[1], {e[0], e[1]}, explicit)}; q11 := eval(p, ans1[1]); rts11 := [solve(%, x, explicit)]

{{e[0] = -(1/12)*e[2]^2, e[1] = -(2/9)*(-6*e[2])^(1/2)*e[2]}, {e[0] = -(1/12)*e[2]^2, e[1] = (2/9)*(-6*e[2])^(1/2)*e[2]}}

x^4+x^2*e[2]-(2/9)*x*(-6*e[2])^(1/2)*e[2]-(1/12)*e[2]^2

[(1/6)*(-6*e[2])^(1/2), (1/6)*(-6*e[2])^(1/2), (1/6)*(-6*e[2])^(1/2), -(1/2)*(-6*e[2])^(1/2)]

The rts11 subcase polynomial q11 must have real coefficients and therefore only applies for e[2]<0. The roots are real with multiplicity 3 and multiplicity 1 and this case is just case (g) above. The rts12 subcase polynomial q12 (below) also requires e[2]<0 and corresponds to case (g), but with signs reversed.

q12 := eval(p, ans1[2]); rts12 := [solve(%, x, explicit)]

x^4+x^2*e[2]+(2/9)*x*(-6*e[2])^(1/2)*e[2]-(1/12)*e[2]^2

[(1/2)*(-6*e[2])^(1/2), -(1/6)*(-6*e[2])^(1/2), -(1/6)*(-6*e[2])^(1/2), -(1/6)*(-6*e[2])^(1/2)]

Therefore the rts2 case is the solution for case (a); the cell37 example was a special case of this. In fact if we have a duplicate pair of complex-conjugate roots, the polynomial must be a pefect square, as we see it is

factor(q1)

(1/4)*(2*x^2+e[2])^2

NULL

Download RootAnalysis5.mw

 

Any explanation why this happens? notice, I did not supply the x and y ranges, let Maple decide.

restart;

interface(version);

`Standard Worksheet Interface, Maple 2026.1, Windows 10, April 28 2026 Build ID 2011354`

plots:-contourplot(y+sin(x),'colorbar'=false,':-contours' = 2,size=[100,100]);

plots:-contourplot(y+sin(x),'colorbar'=false,':-contours' = 1,size=[100,100]);

Error, (in plot/iplot2d) numeric exception: division by zero

 

 

Download bug_in_contourplot.mw

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