Algebraically solving system of nonlinear equations with parameters

asked 2022-11-21 13:11:58 +0100

Alesekk
13 ●3

updated 2023-07-07 08:48:48 +0100

FrédéricC

5251 ●3 ●44 ●116

Hello,

I am kind of new to the sage and I would like to solve system of two nonlinear equations that look like this:

a,b,c,d = var('a,b,c,d')  
solve([c==a*(a/(a+2*b))+2*b*(a/(a+b)),d==a*(b/(a+2*b))+b*(b/(a+b))],a,b)

I have an idea how the result looks like based on Wolfram alpha, but I would like to verify it using sage. The problem is, that sage only gives me anwer a=0 and b=0, which is clearly not the solution, because I would have divide by zero error. Do you have any ideas what to use or how to solve it? Thanks for all replies.

edit retag flag offensive close merge delete

add a comment

2 Answers

Sort by » oldest newest most voted

answered 2022-11-26 10:49:14 +0100

vdelecroix
7645 ●20 ●85 ●168 http://www.labri.fr/pe...

Your system is an algebraic system that could be solved with variable elimination. The first step consists in building the ideal of solutions.

sage: R.<a,b,c,d> = QQ[]
sage: F1 = c - (a*(a/(a+2*b))+2*b*(a/(a+b)))
sage: F2 = d - (a*(b/(a+2*b))+b*(b/(a+b)))
sage: I = R.ideal([F1.numerator(), F2.numerator()]).radical()
sage: I
Ideal (a*b + 2*b^2 - b*c - 2*b*d,
         a^2 - 4*b^2 - a*c + 2*b*c - 2*a*d + 4*b*d,
         2*b^3 - 2*b^2*c + b*c^2 - 4*b^2*d - a*c*d + 3*b*c*d - 2*a*d^2 + 2*b*d^2)
of Multivariate Polynomial Ring in a, b, c, d over Rational Field

Now, you want to know the expression of a in terms of b,d. To do so, just eliminate b

sage: Ib = I.elimination_ideal([b])
sage: Ib
Ideal (a^4 - a^3*c + a^2*c^2 - a*c^3 - 2*a^3*d + 6*a^2*c*d - 4*a*c^2*d + 8*a^2*d^2 - 4*a*c*d^2) of Multivariate Polynomial Ring in a, b, c, d over Rational Field
sage: Ib.gens()[0]
a^4 - a^3*c + a^2*c^2 - a*c^3 - 2*a^3*d + 6*a^2*c*d - 4*a*c^2*d + 8*a^2*d^2 - 4*a*c*d^2
sage: Ib.gens()[0] // a
a^3 - a^2*c + a*c^2 - c^3 - 2*a^2*d + 6*a*c*d - 4*c^2*d + 8*a*d^2 - 4*c*d^2

The last polynomial is the (implicit) expression of a in terms of c and d. Since this is a cubic polynomial in a, you could use the Cardano formula (see https://en.wikipedia.org/wiki/Cubic_e...). This is what you obtain when using solve on a cubic equation

sage: poly = Ib.gens()[0] // a
sage: poly
a^3 - a^2*c + a*c^2 - c^3 - 2*a^2*d + 6*a*c*d - 4*c^2*d + 8*a*d^2 - 4*c*d^2
sage: solve([SR(poly) == 0], SR.var('a'))
[a == -1/3*(1/2)^(2/3)*((c + 2*d)^2 - 3*c^2 - 18*c*d - 24*d^2)*(-I*sqrt(3) + 1)/(2*(c + 2*d)^3 + 27*c^3 + 108*c^2*d + 108*c*d^2 + 36*sqrt(1/3)*sqrt((c^3 + 2*c*d^2 + 14*d^3)*(c + 2*d))*(c + 2*d) - 9*(c^2 + 6*c*d + 8*d^2)*(c + 2*d))^(1/3) - 1/6*(1/2)^(1/3)*(2*(c + 2*d)^3 + 27*c^3 + 108*c^2*d + 108*c*d^2 + 36*sqrt(1/3)*sqrt((c^3 + 2*c*d^2 + 14*d^3)*(c + 2*d))*(c + 2*d) - 9*(c^2 + 6*c*d + 8*d^2)*(c + 2*d))^(1/3)*(I*sqrt(3) + 1) + 1/3*c + 2/3*d,
 a == -1/3*(1/2)^(2/3)*((c + 2*d)^2 - 3*c^2 - 18*c*d - 24*d^2)*(I*sqrt(3) + 1)/(2*(c + 2*d)^3 + 27*c^3 + 108*c^2*d + 108*c*d^2 + 36*sqrt(1/3)*sqrt((c^3 + 2*c*d^2 + 14*d^3)*(c + 2*d))*(c + 2*d) - 9*(c^2 + 6*c*d + 8*d^2)*(c + 2*d))^(1/3) - 1/6*(1/2)^(1/3)*(2*(c + 2*d)^3 + 27*c^3 + 108*c^2*d + 108*c*d^2 + 36*sqrt(1/3)*sqrt((c^3 + 2*c*d^2 + 14*d^3)*(c + 2*d))*(c + 2*d) - 9*(c^2 + 6*c*d + 8*d^2)*(c + 2*d))^(1/3)*(-I*sqrt(3) + 1) + 1/3*c + 2/3*d,
 a == 2/3*(1/2)^(2/3)*((c + 2*d)^2 - 3*c^2 - 18*c*d - 24*d^2)/(2*(c + 2*d)^3 + 27*c^3 + 108*c^2*d + 108*c*d^2 + 36*sqrt(1/3)*sqrt((c^3 + 2*c*d^2 + 14*d^3)*(c + 2*d))*(c + 2*d) - 9*(c^2 + 6*c*d + 8*d^2)*(c + 2*d))^(1/3) + 1/3*c + 2/3*d + 1/3*(1/2)^(1/3)*(2*(c + 2*d)^3 + 27*c^3 + 108*c^2*d + 108*c*d^2 + 36*sqrt(1/3)*sqrt((c^3 + 2*c*d^2 + 14*d^3)*(c + 2*d))*(c + 2*d) - 9*(c^2 + 6*c*d + 8*d^2)*(c + 2*d))^(1/3)]

edit flag offensive delete link

Comments

@vdelecroix :

Bravo, bravo, bravissimo ! But ther's a fly in the ointment :

The same process can be used to eliminate a and find a cubic whose three roots are candidate solutions for b. The problem is that Sage is unable to check any of the nine candidate pairs of solutions : the result of substituting them in the original rational fraction never satisfies is_zero(). Using numerical approximations, one can select three of these pairs as "numerically reasonable", but not prove them correct.

OTOH, the three Mathematica-generated solutions do check (a bit slowly...).

Is this worth opening a ticket ?

Emmanuel Charpentier ( 2022-11-27 18:52:35 +0100 )edit

It is true that we would prefer pairs (A(c,d), B(c,d)) of algebraic functions for this instance. As you said there is a subtle pairing of the roots of these two cubic polynomials. Note that checking numerically is enough: if you have three roots a1, a2, a3 for a and b1, b2, b3 for b you can use numerics to prove that (a1, b2) and (a1, b3) are wrong (ie does not satisfy the system for a particular choice of c, d).

That being said I would not buy the sage symbolic answer as it involves a lot of square roots and third roots that are not carefully determined.

vdelecroix ( 2022-11-30 13:46:25 +0100 )edit

add a comment

answered 2022-11-24 11:09:16 +0100

Emmanuel Charpentier
7774 ●6 ●50 ●153

Partial (non-)answer, to be completed if divine inspiration strikes ;-)...

Solving this system in SR doesn't seem to be possible.

sage: var("a, b, c, d")
(a, b, c, d)
sage: eq1 = c==a*(a/(a+2*b))+2*b*(a/(a+b)) ; eq2 = d==a*(b/(a+2*b))+b*(b/(a+b))

Solving the first equation for b is trivial, leading to two solutions

sage: Solb1 = solve(eq1, b, solution_dict=True)
sage: Solb1
[{b: -1/4*(3*a^2 - 3*a*c + sqrt(-7*a^2 + 6*a*c + c^2)*a)/(2*a - c)},
 {b: -1/4*(3*a^2 - 3*a*c - sqrt(-7*a^2 + 6*a*c + c^2)*a)/(2*a - c)}]

Substituting this solution in the second equation leads to non-linear equations, e. g. :

sage: eq3 = eq2.subs(Solb1[0])
sage: eq3
d == -1/2*(3*a^2 - 3*a*c + sqrt(-7*a^2 + 6*a*c + c^2)*a)*a/((2*a - c)*(2*a - (3*a^2 - 3*a*c + sqrt(-7*a^2 + 6*a*c + c^2)*a)/(2*a - c))) + 1/4*(3*a^2 - 3*a*c + sqrt(-7*a^2 + 6*a*c + c^2)*a)^2/((4*a - (3*a^2 - 3*a*c + sqrt(-7*a^2 + 6*a*c + c^2)*a)/(2*a - c))*(2*a - c)^2)

This equation cannot be solved by Maxima's solver : calling it without extra arguments fails :

sage: Solam = solve(eq3, a, algorithm="maxima", solution_dict=True) # Apparent success...
sage: len(Solam)
4

This solution isn't a solution :

sage: [u[a].variables() for u in Solam]
[(), (a, c, d), (a, c, d), (a, c, d)]

The three last "solutions" are just rewritings of eq3 ; the last is a triviality :

sage: Solam[0]
{a: 0}

leading to `{b:0} when substituting in the first equation, unacceptable as noted before...

Calling

sage: Solamp = solve(eq3, a, algorithm="maxima", solution_dict=True, to_poly_solve=True)

didnt't return after > 30 minutes

Similarly, algorithm="sympy" didn't return after > 30 minutes. Trying giac fails (Giac's solver seems limited to polynomials...). Frica's solver gives an answer seeminlgly using intermediary computations, which do not seem to be capturable "simply" by Sage (to be explored further...).

Another possibility is to use Sage's apparatus for polynomials :

sage: R1.<c0, c1>=QQbar[]
sage: R2.<c2, c3>=FractionField(R1)[]
sage: J=R2.ideal([R2(u.subs(dict(zip((a, b, c, d), (c0,c1, c2, c3))))) for u in (eq1, eq2)])

This ideal, which is, in some sense, the "best" solution to your problem, can be defined in Sage

sage: J
Ideal (c2 + (-c0^3 + (-3)*c0^2*c1 + (-4)*c0*c1^2)/(c0^2 + 3*c0*c1 + 2*c1^2), c3 + (-c0^2*c1 + (-2)*c0*c1^2 + (-2)*c1^3)/(c0^2 + 3*c0*c1 + 2*c1^2)) of Multivariate Polynomial Ring in c2, c3 over Fraction Field of Multivariate Polynomial Ring in c0, c1 over Algebraic Field

One can even be sure that this solution is a finite set :

sage: J.dimension()
verbose 0 (3998: multi_polynomial_ideal.py, groebner_basis) Warning: falling back to very slow toy implementation.
verbose 0 (1082: multi_polynomial_ideal.py, dimension) Warning: falling back to very slow toy implementation.
0

Alas :

sage: J.variety()

[ Snip... ]

NotImplementedError: root finding for this polynomial not implemented

OTOH, I'm far from an expert in this domain, and might have gone in completely wrong direction...

My tentative concusion isn't that Sage cannot solve this system, but rather that I can't solve this system with Sage but by using its mathematica interface to the Wolfram engine...

edit flag offensive delete link

add a comment

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Algebraically solving system of nonlinear equations with parameters

2 Answers

Comments

Your Answer

Question Tools

Stats

Related questions

Algebraically solving system of nonlinear equations with parameters edit

2 Answers

Comments

Your Answer

Question Tools

Stats

Related questions

Algebraically solving system of nonlinear equations with parameters