dharr

Dr. David Harrington

4449 Reputation

20 Badges

18 years, 281 days
University of Victoria
Professor or university staff
Victoria, British Columbia, Canada

Social Networks and Content at Maplesoft.com

Maple Application Center
I am a professor of chemistry at the University of Victoria, BC, Canada, where my research areas are electrochemistry and surface science. I have been a user of Maple since about 1990.

MaplePrimes Activity


These are replies submitted by dharr

@Christian Wolinski Definitely unexpected to get a simpler RootOf here - documentation for convert,radical says "The conversion can fail if Maple cannot find radical expressions for the roots or if the correct radical expression cannot be selected. If the conversion fails, the RootOf remains unchanged." [my bold].

@sursumCorda If you know which RootOf you want the answer in terms of, then Algfield will help.

restart;

a:=[RootOf(_Z^3 - 3*_Z^2 - 10*_Z - 1, index = 3), -4/5*RootOf(_Z^3 - 3*_Z^2 - 10*_Z - 1, index = 3)^2 + 19/5*RootOf(_Z^3 - 3*_Z^2 - 10*_Z - 1, index = 3) + 3/5, 1];
b:=[RootOf(_Z^3 - 3*_Z^2 - 10*_Z - 1, index = 3), RootOf(_Z^3 + 10*_Z^2 + 3*_Z - 1, index = 3), 1];
c:=[RootOf(_Z^3 - 4*_Z^2 + _Z + 1, index = 1)/RootOf(_Z^3 - 4*_Z^2 + _Z + 1, index = 3), RootOf(_Z^3 - 4*_Z^2 + _Z + 1, index = 2)/RootOf(_Z^3 - 4*_Z^2 + _Z + 1, index = 3), 1];

[RootOf(_Z^3-3*_Z^2-10*_Z-1, index = 3), -(4/5)*RootOf(_Z^3-3*_Z^2-10*_Z-1, index = 3)^2+(19/5)*RootOf(_Z^3-3*_Z^2-10*_Z-1, index = 3)+3/5, 1]

[RootOf(_Z^3-3*_Z^2-10*_Z-1, index = 3), RootOf(_Z^3+10*_Z^2+3*_Z-1, index = 3), 1]

[RootOf(_Z^3-4*_Z^2+_Z+1, index = 1)/RootOf(_Z^3-4*_Z^2+_Z+1, index = 3), RootOf(_Z^3-4*_Z^2+_Z+1, index = 2)/RootOf(_Z^3-4*_Z^2+_Z+1, index = 3), 1]

evalf(a)evalf(b)evalf(c)

[-1.924983964, -9.679389670, 1.]

[-1.924983964, -9.679389672, 1.]

[-1.924983964, -9.679389671, 1.]

From a[2] to b[2]

qb := evala(Algfield(a[2], {RootOf(_Z^3+10*_Z^2+3*_Z-1, index = 3)})); evala(eval(a[2], qb[1]))

RootOf(_Z^3+10*_Z^2+3*_Z-1, index = 3)

From a to c

qc := evala(Algfield(a, {RootOf(_Z^3-4*_Z^2+_Z+1, index = 1), RootOf(_Z^3-4*_Z^2+_Z+1, index = 2), RootOf(_Z^3-4*_Z^2+_Z+1, index = 3)})); evala(eval(a, qc[1]))

[RootOf(_Z^3-4*_Z^2+_Z+1, index = 1)^2-2*RootOf(_Z^3-4*_Z^2+_Z+1, index = 1)-1, 3*RootOf(_Z^3-4*_Z^2+_Z+1, index = 1)^2-10*RootOf(_Z^3-4*_Z^2+_Z+1, index = 1)-4, 1]

I think Maple will always prefer a polynomial over a rational function. Notice here that since they can all be written in terms of one index value, so that is simpler. The splitting field for this polynomial has only one algebraic number (by default index=1)

PolynomialTools:-Split(_Z^3-4*_Z^2+_Z+1, _Z, 'K'); K

{RootOf(_Z^3-4*_Z^2+_Z+1)}

But you can check that c[1] is the same as the above result.

evala(RootOf(_Z^3-4*_Z^2+_Z+1, index = 1)/RootOf(_Z^3-4*_Z^2+_Z+1, index = 3))

RootOf(_Z^3-4*_Z^2+_Z+1, index = 1)^2-2*RootOf(_Z^3-4*_Z^2+_Z+1, index = 1)-1

NULL

Download Algfield.mw

@Zeineb The Routh-Hurwitz 'array'  treatment is equivalent so will lead to the same answer. If you got the polynomial from a matrix, then working from the matrix (or even the original differential equations) can sometimes be easier. For example a symmetric matrix proves real roots, but that is harder to tell from the polynomial.

@Axel Vogt Creative approach! (gfun is now a standard package within Maple.)

@C_R I branched it off to here

@C_R Interesting. Yes, the index=real[4] RootOf does not behave like a typical RootOf in some respects. Normally evalf on a RootOf would give a numerical value but here it returns unevaluated.  This would normally be a way to work toward an exact solution, since convert(RootOf(...,numeric value),radical) could then be used.

@sursumCorda A combination of implicitplot (plot_real_curves is sometimes too slow) and algcurves:-singularities can track down the roots. 

restart

with(algcurves)

expr := a^3*b^3*(12*(a*b+a+b)-28*(a+b+1)^2)+4*a^2*b^2*(a+b+1)*(a^2+b^2+(a+1)*(b+1))*(13*(a+b+1)^2-8*(a*b+a+b))+a*b*(7*(a+b+1)^8-40*(a*b+a+b)*(a+b+1)^6+32*(a*b+a+b)^2*(a+b+1)^4+28*(a*b+a+b)^3*(a+b+1)^2-20*(a*b+a+b)^4)+(a+b+1)*(a*b+a+b)*((a+b+1)^4-6*(a*b+a+b)*(a+b+1)^2+6*(a*b+a+b)^2)^2:

plots:-implicitplot(expr, a = -5 .. 5, b = -5 .. 5)

Curve is symmetric - easy to find solutions with a=b. There are three non-negative ones.

simplify(eval(expr, b = a)); solve(%); evalf(%); solns := {a = 0, b = 0}, {a = 1, b = 1}, {a = (1/2)^(1/4), b = (1/2)^(1/4)}

8*(a-1)^2*a*(a^4-1/2)^2

0, 1, 1, (1/2)*2^(3/4), ((1/2)*I)*2^(3/4), -(1/2)*2^(3/4), -((1/2)*I)*2^(3/4), (1/2)*2^(3/4), ((1/2)*I)*2^(3/4), -(1/2)*2^(3/4), -((1/2)*I)*2^(3/4)

0., 1., 1., .8408964155, .8408964155*I, -.8408964155, -.8408964155*I, .8408964155, .8408964155*I, -.8408964155, -.8408964155*I

{a = 0, b = 0}, {a = 1, b = 1}, {a = (1/2)*2^(3/4), b = (1/2)*2^(3/4)}

The plot shows that there are two other solutions along each of the a and b axes - these are the index=1 and index=2 ones

simplify(eval(expr, b = 0)); q := solve(%); evalf(%)

(a+1)*a*(a^4-2*a^3-2*a+1)^2

-1, 0, RootOf(_Z^4-2*_Z^3-2*_Z+1, index = 1), RootOf(_Z^4-2*_Z^3-2*_Z+1, index = 2), RootOf(_Z^4-2*_Z^3-2*_Z+1, index = 3), RootOf(_Z^4-2*_Z^3-2*_Z+1, index = 4), RootOf(_Z^4-2*_Z^3-2*_Z+1, index = 1), RootOf(_Z^4-2*_Z^3-2*_Z+1, index = 2), RootOf(_Z^4-2*_Z^3-2*_Z+1, index = 3), RootOf(_Z^4-2*_Z^3-2*_Z+1, index = 4)

-1., 0., .4354205447, 2.296630263, -.3660254038+.9306048591*I, -.3660254038-.9306048591*I, .4354205447, 2.296630263, -.3660254038+.9306048591*I, -.3660254038-.9306048591*I

solns := solns, {a = 0, b = convert(q[3], radical)}, {a = convert(q[3], radical), b = 0}, {a = 0, b = convert(q[4], radical)}, {a = convert(q[4], radical), b = 0}

{a = 0, b = 0}, {a = 1, b = 1}, {a = (1/2)*2^(3/4), b = (1/2)*2^(3/4)}, {a = 0, b = 1/2+(1/2)*3^(1/2)-(1/2)*2^(1/2)*3^(1/4)}, {a = 1/2+(1/2)*3^(1/2)-(1/2)*2^(1/2)*3^(1/4), b = 0}, {a = 0, b = 1/2+(1/2)*3^(1/2)+(1/2)*2^(1/2)*3^(1/4)}, {a = 1/2+(1/2)*3^(1/2)+(1/2)*2^(1/2)*3^(1/4), b = 0}

Any other points are not part of the main curve, so are singularities. Here they are in homogeneous coordinates, but since the last coordinate is 1, we can just take the first two coordinates. We already had the a=1,b=1 solution and the a=(1/2)^(1/4),b=(1/2)^(1/4) one

s := singularities(expr, a, b)

{[[1, 1, 1], 2, 1, 2], [[1, RootOf(_Z^4-2), 1], 2, 1, 2], [[RootOf(_Z^4-2), 1, 1], 2, 1, 2], [[RootOf(2*_Z^4-1), RootOf(2*_Z^4-1), 1], 2, 1, 2]}

simplify(eval(expr, {a = 1, b = 2^(1/4)}))

0

So the complete set is

solns := solns, {a = 1, b = 2^(1/4)}, {a = 2^(1/4), b = 1}

{a = 0, b = 0}, {a = 1, b = 1}, {a = (1/2)*2^(3/4), b = (1/2)*2^(3/4)}, {a = 0, b = 1/2+(1/2)*3^(1/2)-(1/2)*2^(1/2)*3^(1/4)}, {a = 1/2+(1/2)*3^(1/2)-(1/2)*2^(1/2)*3^(1/4), b = 0}, {a = 0, b = 1/2+(1/2)*3^(1/2)+(1/2)*2^(1/2)*3^(1/4)}, {a = 1/2+(1/2)*3^(1/2)+(1/2)*2^(1/2)*3^(1/4), b = 0}, {a = 1, b = 2^(1/4)}, {a = 2^(1/4), b = 1}

NULL

Download algcurves.mw

For the second one, there are only singularities, and I'll let you break them down into all the subcases.

algcurves2.mw

@C_R  Interesting about index = real[4]. I've not seen that before in many years using Maple, and the stuff on the forum isn't clear about it being a bug or not. Perhaps it is a bug in the documentation. There do seem to be just six solutions, from plotting the curve with algcurves.

restart

with(algcurves)

expr := 36*a^3*b^3+8*a^2*b^2*(9*(a+b+1)^2+4*(a*b+a+b))*(a+b+1)+((a-b)^2-2*(a+b)+1)^2*(a+b+1)^5+a*b*((a-b)^2-2*(a+b)+1)*(17*(a+b+1)^2+4*(a*b+a+b))*(a+b+1)^2:

Curve is symmetric - easy to find solutions with a=b

solve(eval(expr, a = b)); evalf(%)

1, 1, 1/2-(1/2)*3^(1/2), 1/2+(1/2)*3^(1/2), 1/2-(1/2)*3^(1/2), 1/2+(1/2)*3^(1/2)

1., 1., -.3660254040, 1.366025404, -.3660254040, 1.366025404

The positive points here points seem to be isolated from the main curve - and there are 6 points in the closed first quadrant.

plot_real_curve(expr, a, b, abserr = 1.*10^(-3), view = [-5 .. 5, -5 .. 5])

 

Download algcurves.mw

@Christian Wolinski I sumitted an SCR.

Use the green up-arrow to upload the worksheet.

@tomleslie @mmcdara Well, @nm asked for "a more elegent way to remove any entry in piecewise which has undefined in it". I didn't really answer that for the general case, but did for the case in question with a Heaviside function. As always, Maple considers the general case and one needs to be careful in special cases. I'm sure @nm will consider that. However, in the context of laplace transforms outputting the unit step function is fine because integrals are involved. In fact taking the Laplace transform of the different outputs gives the same answer:

restart;

with(inttrans):
expr:=exp(-s)/s;
t1:=invlaplace(expr,s,t);
p1:=convert(t1,piecewise);
_EnvUseHeavisideAsUnitStep:=true;
p2:=convert(t1,piecewise);
laplace(t1,t,s);
laplace(p1,t,s);
laplace(p2,t,s);

expr := exp(-s)/s

t1 := Heaviside(t-1)

p1 := piecewise(t < 1, 0, t = 1, undefined, 1 < t, 1)

_EnvUseHeavisideAsUnitStep := true

piecewise(t < 1, 0, 1 <= t, 1)

exp(-s)/s

exp(-s)/s

exp(-s)/s

 

Download Heaviside.mw

For the general case, Maple straddles the two extremes of treating Heaviside as a function or a distribution, and so for other applications I agree the unit step function might not be appropriate.

@Rouben Rostamian  In Maple 2022, ?coords takes me to a very helpful page "Coordinate Systems Supported in Maple" which makes this clear. (This page wasn't so easy to find in the online help.):

"(**): There are two widely used conventions for how to represent spherical coordinates: one is used more in physics and the other more in mathematics. The only difference is a swap of the second and third coordinate: the vector represented by 
                           "u, v, w"

 in one convention is represented by 
                           "u, w, v"

 in the other. The plotting library, the changecoords command, and the VectorCalculus package allow you to use the coordinate system names spherical_math and spherical_physics to refer to these coordinate systems unambiguously. You can also use the name spherical; its behavior is as follows:
– If neither of the Physics and VectorCalculus packages has been loaded (using the with command), then the plotting library and the changecoords command will interpret spherical as spherical_math.
– If either of the packages Physics or VectorCalculus has been loaded, then the plotting library and the changecoords command will interpret spherical as spherical_physics.
The VectorCalculus package will always interpret spherical as spherical_physics.
 

@Ahmed111 Well, Maple is often a bit inscrutable about the way it does things. Here's my guess. The first time, Maple always found the lambdas together, e.g., (lambda__1-conjugate(lambda__2)) and so kept them together. But the second time, you have things like phi2 or conjugate(psi2), which now have polynomials in lambda__2 or its conjugate. Now I might look at phi2 and conjugate(psi2) and try to combine them (they are very similar after all, but not conjugates), which might give a shorter expression if split into real and imaginary parts.

If you don't want Maple to combine them, just call them something like lambda__1 and clambda__1 until the end of the calculation and then substitute back.

But at some point, you have to decide what is acceptably simple, and Maple might not agree (or at least without contortions).

@arashghgood You needed list1[i][2] instead of just list1[i], and numeric rather than numerical.

ode.mw

@arashghgood I'm confused about what you want. You talked about residues - do you want to do a contour integration? if so of what function - Q, deq? - which contour? upper half plane? That can probably be done.

How exactly was the plot in c made? I think by numerical solution of the same differential equation? if you want to exactly replicate that then it would be helpful to know the method, stepsize, hardware precision, exact initial conditions.

Perhaps a step back to tell us where is problem came from and what you want out of it would be helpful.

4 5 6 7 8 9 10 Last Page 6 of 42