Advanced Search

dharr

Dr. David Harrington

6506 Reputation

21 Badges

20 years, 47 days
University of Victoria
Professor or university staff
Victoria, British Columbia, Canada

Social Networks and Content at Maplesoft.com

Maple Application Center
Contact Dr. David Harrington
I am a retired professor of chemistry at the University of Victoria, BC, Canada. My research areas are electrochemistry and surface science. I have been a user of Maple since about 1990.

MaplePrimes Activity


These are answers submitted by dharr

Using RootOf...

dharr 6506
October 08 2024
0 0

As the others have pointed out the sum involves roots of a degree 4 polynomial. If you put some values in, then this will give you directly what you want. (I used some random values so didn't get realistic output.)

Download RootOf.mw

alternative method...

dharr 6506
October 07 2024
2 1

Not sure what happened. Here's a way to get around that. Probably simpler to combine all constants into A and B at the beginning.

Download odes.mw

convert to set...

dharr 6506
October 05 2024
1 0

Something like this? (The matrices don't show on Mapleprimes; download the worksheet and run to see them.)

> 

restart
> 

with(LinearAlgebra)
> 

A := `<,>`(`<|>`(2*x, z), `<|>`(-2*x-2, x-2)); B := `<,>`(`<|>`(4*y, 2*y), `<|>`(w, 3))

Matrix(2, 2, {(1, 1) = 2*x, (1, 2) = z, (2, 1) = -2*x-2, (2, 2) = x-2})

Matrix(%id = 36893490406780077108)

> 

Meqns := A-B-IdentityMatrix(2)

Matrix(%id = 36893490406780082884)

By default solve assumes =0

> 

eqns := convert(Meqns, set); solve(eqns)

{x-6, z-2*y, -2*x-2-w, 2*x-4*y-1}

{w = -14, x = 6, y = 11/4, z = 11/2}

 

 

Download eqns.mw

(For the more usual A.x = b with A a matrix of constants, x a vector of variables and b a vector of constants, use LinearSolve)

fsolve...

dharr 6506
October 03 2024
1 4

For future reference, please upload your worksheet with the green up-arrow in the Mapleprimes editor.
Your nested piecewise expression can be changed to a regular piecewise by simplify(..., piecewise).
You need to solve for y not y_.
Since you want a numerical solution, use fsolve, not solve.
Then it solves and gives the solution 0. You are asking when your piecewise expression 2*AA_=A_, but one of the conditions on each side is that it is 0 for y<=0, and so y=0 is a correct solution.

fsolve.mw

can't be done in that form...

dharr 6506
October 03 2024
1 2
> 

restart;
> 

f:=(deltab*(-2*ub^2*mu*M^2*deltab*((-2+deltab)*(3/4+(phi0-1/2)*kappa)*(kappa-1/2)*M^4-(5*((alpha^2-6/5)*phi0*(kappa-1/2)*deltab+((-9/5*phi0-1/5)*alpha^2+13/5*phi0)*kappa+(11/10*phi0+3/10)*alpha^2-13/10*phi0))*(3/2+phi0-kappa)*M^2-2*alpha^2*phi0*(3/2+phi0-kappa)^2)*(3/2+phi0-kappa)*sqrt(1/(mu*ub^2))+sqrt(2)*(((3*(((mu*ub^2+2/3)*kappa^2+(-(2/3*mu)*ub^2-1)*kappa+(1/12)*mu*ub^2)*deltab^2-4*ub^2*mu*(kappa-1/6)*(kappa-1/2)*deltab+4*ub^2*mu*(kappa-1/6)*(kappa-1/2)))*(kappa-1/2)*M^6+(10*((alpha^2-6/5)*((mu*ub^2+1)*kappa-(1/2)*mu*ub^2-3/2)*deltab^2-(14/5*(ub^2))*mu*((alpha^2-12/7)*kappa-13/14*(alpha^2)+1)*deltab+(8/5*((alpha^2-3)*kappa-2*alpha^2+2))*ub^2*mu))*(kappa-3/2)*(kappa-1/2)*M^4-8*ub^2*mu*(alpha^2*(kappa-1/2)*deltab-3*alpha^2*kappa-(1/2)*kappa+1/4)*(kappa-3/2)^2*M^2-8*ub^2*mu*alpha^2*(kappa-3/2)^3)*(-phi0)^(3/2)+(2*kappa^2*deltab^2*(kappa-1/2)*M^6+(10*((alpha^2-6/5)*(2*kappa^2+(mu*ub^2-2)*kappa-(1/2)*mu*ub^2-3/2)*deltab^2-(14/5*(ub^2))*((alpha^2-12/7)*kappa-6/7*(alpha^2)+15/14)*mu*deltab+(8/5*(ub^2))*((alpha^2-3)*kappa-7/4*(alpha^2)+9/4)*mu))*(kappa-1/2)*M^4-(16*(-(1/8)*alpha^2*(kappa-3/2)*(kappa-1/2)*deltab^2+mu*ub^2*alpha^2*(kappa-1/2)*deltab-3*ub^2*mu*((alpha^2+1/6)*kappa-(1/4)*alpha^2-1/12)))*(kappa-3/2)*M^2-24*ub^2*mu*alpha^2*(kappa-3/2)^2)*(-phi0)^(5/2)+(20*kappa*(alpha^2-6/5)*deltab^2*(kappa-1/2)*M^4+(4*alpha^2*(kappa-3/2)*(kappa-1/2)*deltab^2-8*mu*ub^2*alpha^2*(kappa-1/2)*deltab+24*ub^2*((alpha^2+1/6)*kappa-(1/3)*alpha^2-1/12)*mu)*M^2-24*ub^2*mu*alpha^2*(kappa-3/2))*(-phi0)^(7/2)+(2*(deltab^2*(kappa-1/2)*M^2-4*mu*ub^2))*alpha^2*(-phi0)^(9/2)+((((mu*ub^2+1/2)*kappa-(1/2)*mu*ub^2-3/4)*deltab^2-4*ub^2*mu*(kappa-1/2)*deltab+4*ub^2*mu*(kappa-1/2))*(kappa-1/2)*M^4+2*ub^2*mu*(alpha^2+1)*(-2+deltab)*(kappa-3/2)*(kappa-1/2)*M^2+4*ub^2*mu*alpha^2*(kappa-3/2)^2)*M^2*sqrt(-phi0)*(kappa-3/2)))*sqrt(-phi0/(mu*ub^2))+(1/2*((((mu*phi0*ub^2-2*(phi0-1/2)^2)*kappa^2+(3/2-3*phi0-mu*phi0*ub^2)*kappa-9/8+(1/4)*mu*phi0*ub^2)*deltab^2-4*ub^2*mu*phi0*(kappa-1/2)^2*deltab+4*ub^2*mu*phi0*(kappa-1/2)^2)*M^4-(2*(-(10*(alpha^2-6/5))*(3/4+(phi0-1/2)*kappa)*deltab^2+ub^2*mu*(alpha^2+1)*(kappa-1/2)*deltab-2*ub^2*mu*(alpha^2+1)*(kappa-1/2)))*phi0*(3/2+phi0-kappa)*M^2+4*alpha^2*(mu*ub^2-(1/2)*deltab^2*phi0)*phi0*(3/2+phi0-kappa)^2))*M^2*deltab*sqrt(2)*(3/2+phi0-kappa)*sqrt(1/(mu*ub^2))+(-(((mu*ub^2+4)*kappa^2+(-mu*ub^2-7)*kappa+(1/4)*mu*ub^2+3/2)*deltab^2-4*ub^2*mu*(kappa-1/2)^2*deltab+4*ub^2*mu*(kappa-1/2)^2)*(-2+deltab)*(kappa-1/2)*M^6-(2*((5*(alpha^2-6/5))*(kappa-3/2)*(kappa-1/2)*deltab^3+((ub^2*(alpha^2+2)*mu-8*alpha^2+14)*kappa^2+(-ub^2*(alpha^2+2)*mu-53/2+16*alpha^2)*kappa+(1/4)*ub^2*(alpha^2+2)*mu-6*alpha^2+33/4)*deltab^2-4*ub^2*mu*(alpha^2+2)*(kappa-1/2)^2*deltab+4*ub^2*mu*(alpha^2+2)*(kappa-1/2)^2))*(kappa-3/2)*M^4-(8*(-(1/2)*alpha^2*(kappa-3/2)*deltab^2+ub^2*(alpha^2+1/2)*mu*(kappa-1/2)*deltab-2*ub^2*(alpha^2+1/2)*mu*(kappa-1/2)))*(kappa-3/2)^2*M^2-8*ub^2*mu*alpha^2*(kappa-3/2)^3)*(-phi0)^(3/2)+(-6*kappa*(-2+deltab)*deltab^2*(kappa-1/3)*(kappa-1/2)*M^6+(-40*kappa*(alpha^2-6/5)*(kappa-3/2)*(kappa-1/2)*deltab^3+((60*alpha^2-88)*kappa^3+(-2*ub^2*(alpha^2+2)*mu-134*alpha^2+180)*kappa^2+(2*ub^2*(alpha^2+2)*mu+73*alpha^2-77)*kappa-(1/2)*ub^2*(alpha^2+2)*mu-21/2*(alpha^2)+15/2)*deltab^2+8*ub^2*mu*(alpha^2+2)*(kappa-1/2)^2*deltab-8*ub^2*mu*(alpha^2+2)*(kappa-1/2)^2)*M^4-(16*(kappa-3/2))*((1/8)*alpha^2*(kappa-3/2)*(kappa-1/2)*deltab^3-(5/4*(alpha^2))*(kappa+1/10)*(kappa-3/2)*deltab^2+ub^2*(alpha^2+1/2)*mu*(kappa-1/2)*deltab-2*ub^2*(alpha^2+1/2)*mu*(kappa-1/2))*M^2-24*alpha^2*(((1/6)*kappa-1/4)*deltab^2+mu*ub^2)*(kappa-3/2)^2)*(-phi0)^(5/2)+(-40*deltab^2*((alpha^2-6/5)*(kappa-1/4)*(kappa-1/2)*deltab+(-3/2*(alpha^2)+11/5)*kappa^2+(3/2*(alpha^2)-19/10)*kappa-3/8*(alpha^2)+17/40)*M^4+(-4*alpha^2*(kappa-3/2)*(kappa-1/2)*deltab^3+40*alpha^2*(kappa-1/5)*(kappa-3/2)*deltab^2-8*ub^2*(alpha^2+1/2)*mu*(kappa-1/2)*deltab+16*ub^2*(alpha^2+1/2)*mu*(kappa-1/2))*M^2-24*alpha^2*(kappa-3/2)*(((1/2)*kappa-3/4)*deltab^2+mu*ub^2))*(-phi0)^(7/2)-2*alpha^2*(((kappa-1/2)*deltab-10*kappa+3)*deltab^2*M^2+(-9+6*kappa)*deltab^2+4*mu*ub^2)*(-phi0)^(9/2)-(1/2)*deltab^2*(8*alpha^2*(-phi0)^(11/2)+((-2+deltab)*(kappa-1/2)*M^2-3+2*kappa)*M^4*sqrt(-phi0)*(kappa-3/2)^2))*sqrt(-phi0/(mu*ub^2))*sqrt(-phi0^3)*Omega^2/(sqrt(-phi0)*sqrt(-phi0^3/(mu^3*ub^6))*ub^4*mu^2*(M^2*sqrt(-phi0)*deltab*sqrt(2)*(kappa-1/2)*sqrt(-phi0/(mu*ub^2))-(1/2)*M^2*deltab*sqrt(2)*(3/2+phi0-kappa)*sqrt(1/(mu*ub^2))-2*(-phi0)^(3/2)+(((-kappa+1/2)*deltab+2*kappa-1)*M^2+3-2*kappa)*sqrt(-phi0))^3*M^2):

Required form

> 

f1:=(M^2-M1^2)/(M^2*(M^2-M0^2));

(M^2-M1^2)/(M^2*(M^2-M0^2))

> 

limit(f1*M^2, M=infinity);

1

Find numerator and denominator high and low degrees. Could be (M^6-const)/(M^2*(M^6-const)) perhaps (but no, see below), but can't be of the expected form

> 

f:=normal(f):
nf:=collect(numer(f),M):
degree(nf,M),ldegree(nf,M);
df:=collect(denom(f),M):
degree(df,M),ldegree(df,M);

6, 0

8, 2

Nonzero powers of M in the numerator and denominator for all even powers. We could factor the denominator into M^2(M^6 + ...) but that's about it.

> 

seq([i,evalb(coeff(nf, M, i) <>0)], i = 6..0,-1);
seq([i,evalb(coeff(df, M, i) <>0)], i = 8..0,-1);

[6, true], [5, false], [4, true], [3, false], [2, true], [1, false], [0, true]

[8, true], [7, false], [6, true], [5, false], [4, true], [3, false], [2, true], [1, false], [0, false]

> 

NULL

 

NULL

Download f1.mw

isolve with help...

dharr 6506
September 30 2024
4 0
> 

restart

Equation recast in form for ThueSolve

> 

eq := -N*x*y+x^2+y^2 = N

-N*x*y+x^2+y^2 = N

Try for N=9;

> 

eqN := eval(eq, N = 9)

x^2-9*x*y+y^2 = 9

Finds no solutions - for the quadratic case just passes to isolve, which is not up to the task.

> 

NumberTheory:-ThueSolve(eqN)

[]

> 

isolve(eqN)

Solve just solves the quadratic for arbirary y. But it is clear that if the discriminant were a square and y were even, then x would be integer.

> 

ans := [solve(eqN)]

[{x = (9/2)*y+(1/2)*(77*y^2+36)^(1/2), y = y}, {x = (9/2)*y-(1/2)*(77*y^2+36)^(1/2), y = y}]

Want discriminant a square

> 

deq2 := discrim((lhs-rhs)(eqN), x) = z^2

77*y^2+36 = z^2

This time we get 12 solutions, each in terms of an arbitrary integer variable _Z1

> 

sols := [isolve(deq2)]; nops(%)

12

Look at the first solution, for say _Z1=1

> 

test := simplify(eval(sols[1], _Z1 = 1))

{y = -240, z = -2106}

Find the corresponding x value

> 

test2 := simplify(eval(ans[1], test))

{-240 = -240, x = -27}

And check it is a solution for the original equation

> 

eval(eqN, `union`(test, test2))

9 = 9

From the form of the equation, changing the signs of x and y gives a positive solution.

> 

soln := `~`[`=`]({x, y}, eval({-x, -y}, `union`(test, test2))); eval(eqN, %)

{x = 27, y = 240}

9 = 9

Check that the original form of the equation returns N

> 

eval((x^2+y^2)/(x*y+1), soln)

9

For N=49

> 

eqN := eval(eq, N = 49); ans := `~`[`=`]([x, x], [solve(eqN, x)]); deq2 := discrim((lhs-rhs)(eqN), x) = z^2; sols := [isolve(deq2)]; nops(%); test := simplify(eval(sols[1], _Z1 = 1)); test2 := {simplify(eval(ans[1], test))}; soln := `~`[`=`]({x, y}, eval({-x, -y}, `union`(test, test2))); eval(eqN, soln)

x^2-49*x*y+y^2 = 49

[x = (49/2)*y+(1/2)*(2397*y^2+196)^(1/2), x = (49/2)*y-(1/2)*(2397*y^2+196)^(1/2)]

2397*y^2+196 = z^2

12

{y = -16800, z = -822514}

{x = -343}

{x = 343, y = 16800}

49 = 49

For N=729

> 

eqN := eval(eq, N = 729); ans := `~`[`=`]([x, x], [solve(eqN, x)]); deq2 := discrim((lhs-rhs)(eqN), x) = z^2; sols := [isolve(deq2)]; nops(%); test := simplify(eval(sols[1], _Z1 = 1)); test2 := {simplify(eval(ans[1], test))}; soln := `~`[`=`]({x, y}, eval({-x, -y}, `union`(test, test2))); eval(eqN, soln)

x^2-729*x*y+y^2 = 729

[x = (729/2)*y+(1/2)*(531437*y^2+2916)^(1/2), x = (729/2)*y-(1/2)*(531437*y^2+2916)^(1/2)]

531437*y^2+2916 = z^2

12

{y = -14348880, z = -10460294154}

{x = -19683}

{x = 19683, y = 14348880}

729 = 729

> 

NULL
> 

NULL

Download Thu.mw

Strange square brackets...

dharr 6506
September 23 2024
1 2

For Var__phi you have Omega^2[  ... ] (The ] is at the end of the expression). What do the square brackets mean?

Is this supposed to be Omega^2*( ... ) or some unspecified function Omega( ... )^2 or something else.

intersection point plot...

dharr 6506
September 23 2024
1 0

Your second case can be done by

display(p2, pointplot(eval([x, y], sol[]), symbol = circle, symbolsize = 12, color = black))

Linear_System.mw

.mpl and .mw...

dharr 6506
September 23 2024
1 0

You are confused about the worksheet (*.mw) in which Maple runs, and the external text file of Maple commands that are read in to the worksheet (*.txt). Download the following two files to the same directory. Open the EKHAD.mw file by double-clicking (NOT the one you made) and follow the instructions (it has a bit more explanation). Hope this is clear.

EKHAD.mw

EKHAD.txt

(normally I would call the last one EKHAD.mpl, but the Mapleprimes site won't let me upload that, so EKHAD.txt will do.)

 

centered cube with scaling...

dharr 6506
September 22 2024
1 4

If you are willing to center the cube, then you can use scaling to do this. I used geom3d to easily get the vertices, but you could use cuboid. You can reorder lbls to get the order you need.

> 

restart
> 

with(geom3d); _local(D)
> 

cube(poly, point(cen, [0, 0, 0]), side = 2); vs := vertices(poly); lbls := [A, B, C, D, A[1], B[1], C[1], D[1]]

poly

[[1, 1, 1], [1, 1, -1], [1, -1, 1], [1, -1, -1], [-1, 1, 1], [-1, 1, -1], [-1, -1, 1], [-1, -1, -1]]

[A, B, C, D, A[1], B[1], C[1], D[1]]

Scale coords by 110% and combine with labels.

> 

tplot := zip(proc (coords, lbl) options operator, arrow; [(1.1*coords)[], lbl] end proc, vs, lbls)

[[1.1, 1.1, 1.1, A], [1.1, 1.1, -1.1, B], [1.1, -1.1, 1.1, C], [1.1, -1.1, -1.1, D], [-1.1, 1.1, 1.1, A[1]], [-1.1, 1.1, -1.1, B[1]], [-1.1, -1.1, 1.1, C[1]], [-1.1, -1.1, -1.1, D[1]]]

> 

pcube := draw(poly); ptxt := plots:-textplot3d(tplot, 'font' = ["times", "bold", 20]); plots:-display(pcube, ptxt)

> 

NULL

Download cube.mw

Equivalent...

dharr 6506
September 20 2024
4 1

The standard way of numerically solving second or higher order diiferential equations or systems is to convert them to an equivalent system of first order equations. For example

{diff(y(x), x, x) = -y(x), y(0) = 1, D(y)(0) = 0};

becomes

{diff(w(x), x) = -y(x), w(x) = diff(y(x), x), y(0) = 1, w(0) = 0};

They both have the solution y(x) = cos(x), and the second also has the "intermediate" solution w(x) = -sin(x). This way, the numerical solvers only need to be able to solve first order odes. In fact convertsys is mainly intended as an internal routine used by dsolve; normally you would just pass the higher order odes to dsolve. 

solve,identity...

dharr 6506
September 17 2024
0 2

Your E2 is too complicated, so I'm not sure exactly what is going on. But in principle, I think you want to do something like this:

> 

restart;
> 

eq:=add(a[i]/r^i,i=1..4)=add(rand(1..10)()*x^i/r^i,i=1..4)

a[1]/r+a[2]/r^2+a[3]/r^3+a[4]/r^4 = 7*x/r+10*x^2/r^2+6*x^3/r^3+2*x^4/r^4

> 

solve(identity(eq,r),{seq(a[i],i=1..4)});

{a[1] = 7*x, a[2] = 10*x^2, a[3] = 6*x^3, a[4] = 2*x^4}

> 

NULL

Download solve.mw

partial answer...

dharr 6506
September 17 2024
1 0

There must be at least one typo in the notes. Probably you need to follow the derivation up to there closely to see how to resolve this.

> 

restart
> 

with(LinearAlgebra)

Assume a general form for the matrix. a is all terms that have opposite signs in B[1,1[ and B[2,2]; b is all terms that are the same sign in these entries.

> 

B := Matrix(2, 2, [a+b-Omega, c, c, -a+b+Omega])

Matrix(%id = 36893490796494478140)

> 

solve(Determinant(B), Omega)

a+(b^2-c^2)^(1/2), a-(b^2-c^2)^(1/2)

So it is clear which terms have to have opposite signs -the sign of 4*alpha^2*Delta in A[2,2] is incorrect and is changed (in red). The signs of 4*alpha^2*beta*Delta/Lambdamust the same so one of them must be changed (two possibilities - change in magenta).

> 

A := Matrix([[-2*alpha^2*gamma+beta^2*Delta-4*alpha^2*Delta-4*alpha^2*beta*Delta/Lambda+2*beta*gamma*Lambda+2*beta*Delta*Lambda+gamma*Delta^2+Delta*Lambda^2-Omega, -2*alpha^2*gamma-4*alpha^2*beta*Delta/Lambda], [-2*alpha^2*gamma-4*alpha^2*beta*Delta/Lambda, -2*alpha^2*gamma-beta^2*Delta+4*alpha^2*Delta-4*alpha^2*beta*Delta/Lambda-2*beta*gamma*Lambda+2*beta*Delta*Lambda+gamma*Delta^2-Delta*Lambda^2+Omega]])

Matrix(%id = 36893490796494480436)

> 

s1, s2 := solve(Determinant(A), Omega)

We have the required terms before the sqrt

> 

s1 := expand(s1)

Delta*Lambda^2-4*alpha^2*Delta+beta^2*Delta+2*beta*gamma*Lambda+(Delta^4*Lambda^2*gamma^2+4*Delta^3*Lambda^3*beta*gamma-8*Delta^3*Lambda*alpha^2*beta*gamma+4*Delta^2*Lambda^4*beta^2-16*Delta^2*Lambda^2*alpha^2*beta^2-4*Delta^2*Lambda^2*alpha^2*gamma^2-8*Delta*Lambda^3*alpha^2*beta*gamma)^(1/2)/Lambda

However, the term with the sqrt is not correct. Squaring it must be b^2-c^2 = (b+c)*(b-c). Not clear how to get this result.

> 

simplify(((s1-s2)*(1/2))^2)

Delta*(Delta^2*Lambda*gamma+2*Delta*Lambda^2*beta-8*Delta*alpha^2*beta-4*Lambda*alpha^2*gamma)*(Delta*gamma+2*Lambda*beta)/Lambda

> 

NULL

Download DR_1.mw

 

DocumentTools:-Tabulate...

dharr 6506
September 15 2024
1 0

Here's one way. (alignment = center is actually the default). There are other options you can play with. For me on Windows the export as .pdf renders everything as you see it. (The actual Table is not rendered on Mapleprimes.)

See ?DocumentTools,Tabulate

> 

restart
> 

p1 := plot(x^2, x = -1 .. 1); p2 := plot(x^3, x = -1 .. 1)
> 

DocumentTools:-Tabulate([[p1], [p2]], alignment = center, width = 50)
> 

NULL

Download plots.mw

progress?...

dharr 6506
September 14 2024
1 4

Maple can find the condtion for 3 real roots but it has the inequality the wrong way round. SCR submitted. But then additional conditions for them to be positive seem difficult.

> 

restart
> 

with(SolveTools)
> 

f := p*x^2+x^3+q*x+r

p*x^2+x^3+q*x+r

Returns a piecewise function

> 

sols := Polynomial(f, x, domain = real); nops(%)

4

For the first condition there are 3 real roots, so this condition is what we want.

> 

cond1 := op(1, sols); sol1 := op(2, sols); nops(%)

0 <= (1/27)*q^3-(1/108)*q^2*p^2-(1/6)*q*p*r+(1/4)*r^2+(1/27)*r*p^3

3

> 

simplify(`~`[`*`](cond1, 108))

0 <= 27*r^2+(4*p^3-18*p*q)*r-q^2*p^2+4*q^3

Condition appears to be incorrect - following parameters satisfy the condition but don't give 3 real roots

> 

params := {p = 1, q = 1, r = 1}; eval(cond1, params); solve(eval(f, params))

{p = 1, q = 1, r = 1}

0 <= 4/27

-1, I, -I

(For the second condition, there is only 1 real root, so this is not what we want.)

> 

cond2 := op(3, sols); sol2 := op(4, sols); nops(%)

(1/27)*q^3-(1/108)*q^2*p^2-(1/6)*q*p*r+(1/4)*r^2+(1/27)*r*p^3 < 0

1

Assuming the condition was really that the expression was less than zero, we now we want additional conditions that make the roots positive...

> 

sol1[1]

(1/6)*(36*q*p-108*r-8*p^3+12*(12*p^3*r-3*p^2*q^2-54*p*q*r+12*q^3+81*r^2)^(1/2))^(1/3)-6*((1/3)*q-(1/9)*p^2)/(36*q*p-108*r-8*p^3+12*(12*p^3*r-3*p^2*q^2-54*p*q*r+12*q^3+81*r^2)^(1/2))^(1/3)-(1/3)*p

> 

NULL

Download cubic.mw

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