Solver cannot solve fairly simple system of equations

[mikldk]

Steps to reproduce:

In> OldSolve({(x*y)/2-3*a == 0, x^2/4-(3*a)/2 == 0, 45-(3*x+(3*y)/2) == 0}, {x, y, a})
Out> {{x==(6*a)/y,y==(6*a)/Sqrt(6*a),a==a}}
In> Solve({(x*y)/2-3*a == 0, x^2/4-(3*a)/2 == 0, 45-(3*x+(3*y)/2) == 0}, {x, y, a})
Out> {{x==(6*a)/y,y==(6*a)/Sqrt(6*a),a==a}}

Expected: (By manual/numeric/Wolfram Alpha+%3D%3D+0%7D,+%7Bx,+y,+lambda%7D))):

{{x == 0, y == 30, a == 0}, {x == 10, y == 10, a == 50/3}} 

Not sure if the problem occurs because there are two solutions or what causes this?

Background: Finding extrema of the Lagrangian

L := x^2 * (y/4) - a*(3*x + 3*y/2 - 45)
eq1 := D(x) L
eq2 := D(y) L
eq3 := D(a) L
OldSolve({eq1 == 0, eq2 == 0, eq3 == 0}, {x, y, a})

[grzegorzmazur]

Hi,

From the top of my head (I haven't actually checked it), the type of the system is incorrectly recognized and wrong method to find a solution is applied, which results in plainly wrong answer. I'll try to fix this, but it may take some time (I'm quite swamped now and next week I go for short vacations).

Cheers, Grzesiek

[mikldk]

Okay, thanks. I'll also see if I have any chance of trying to dig into it.

[grzegorzmazur]

OK, at least diagnosis is simple: in solve.rep/code.sy Solve() calls Solve'System() which in turn checks if this is a linear system. Linear systems are solved by Gaussian elimination, and everything else by back substitution. In this particular case an approach based on Groebner basis would be the right one, but it's not there (at least yet :)

[mikldk]

That would indeed be a very useful addition! Unfortunately I'm not sure that I can implement that. But I would be very grateful if you can at some point (nothing urgent).

[mikldk]

I have started digging into a simple form of Buchberger's algorithm a bit. (I was thinking to start with that, and then make various optimisations once a simple, robust, but possibly slow version was implemented. To have something to test against.)

I seems like a few features are missing in yacas for working with multivariate polynomials:

• MultiDegree(): Degree() for a varlist returning a list of same length as varlist
• Reduction/multivariate division
• A few helpers
  • LeadingMonomial(f, var) = var^Degree(f, var)
  • LeadingTerm(f, var) = LeadingCoef(f, var) * LeadingMonomial(f, var)

[grzegorzmazur]

Hi,

Sorry for answering so late, but I left computer behind going for vacations :) Actually, what you want is already kind of implemented in yacas, that is the simplest form of the Buchberger's algorithm (see multivar.rep/code.ys:318). Unfortunately it's been virtually untested, and hangs for your system. Fortunately, it's been easy to fix, see 3645d9c346d59cb81d78c89b74b985c272feef16 :)

I'll try to plug it somehow into the solver in the next days.

Cheers, Grzesiek

[mikldk]

No reason to apologise! I hope you had a nice vacation! Thanks for looking into this again. I can call

Groebner({(x*y)/2-3*a, x^2/4-(3*a)/2, 45-(3*x+(3*y)/2)})

and get

{(x*y)/2-3*a,x^2/4+((-3)*a)/2,(-3)*x+((-3)*y)/2+45,9*y*a-90*a,15*y-y^2/2-6*a,(-6)*a^2+100*a,((-15)*y)/2+(-9)*a+225}

I can then pick the ones that are easiest to solve: Solve((-6)*a^2+100*a == 0, a) gives me a = 50/3, and further WithValue(a, 50/3, Solve(9*y*a-90*a == 0, y)) gives y = 10 to finally obtain WithValue({a, y}, {50/3, 10}, Solve((x*y)/2-3*a == 0, x)) which gives x = 10.

So it does indeed seem to work. Thank you very much! I look forward to when the above happens automatically in Solve() :-).

A related question: I cannot call e.g.

MultiDegree((x*y)/2-3*a)
MultiLT((x*y)/2-3*a)

Are these not "exported"/visible for the user?

[grzegorzmazur]

Hi,

There are two issues here. You are right, they're not exported - you can check this looking at the corresponding .def file. But there is one more thing, the rules you've mentioned are defined for "multinomials", that is for a specific sparse representation of polynomials which you obtain by applying MM() to a polynomial.

Cheers, Grzesiek

[mikldk]

Thanks. That explains it.

But sorry, it leads to to two additional questions:

• There seem to be missing something about MultiDegree(), MultiDivide(), MultiGcd(), Groebner() etc. from the documentation. Should I create an issue to not forget it?
• Is there any reason that MultiDegree() and other "helpers" (e.g. MultiLeadingTerm()) are not exported? They seem to be interesting in themselves? Should I create an issue mentioning this?

[grzegorzmazur]

Hi,

Starting from the latter issue - no, there's no actual reason. They are treated as helpers, but actually they deserve a higher status, exactly as you say. So please yes, or, even better, if you could just go ahead and fix it - should probably take not much more effort.

With respect to the documentation, they are mentioned in https://yacas.readthedocs.io/en/latest/book_of_algorithms/multivar.html#implementation-of-multivariate-polynomials but that's obviously not enough. So if you don't mind creating an issue that would be great.

Cheers, Grzesiek

[mikldk]

Fantastic, thanks heaps!