solutions of simple polynomial equations not exact

[Henning Siebel <henning.siebel at gmx.de>, Tue, 28 Nov 2006 11:06:16 +0000 (UTC)]:
> a: x^2=2*y$; 
> b: y^2=x/2$;

One way is to simply do what you would do on paper:

- eliminate y: eliminate([a,b], [y])

- solve for x: solve(..., x)

- for each solution xx: map(lambda([xx], ...cons), ...)

  - substitute x into a: ev(a,xx)

  - solve for y: solve(..., y)

  - and build a list containing the solution for x and y: cons(xx, ...)

or in one expression:

map(lambda([xx], cons(xx, solve(ev(a,xx), y))),
    solve(eliminate([a,b], [y]), x));
==> [[x = (2^(1/3)*sqrt(3)*%i-2^(1/3))/2,
y = -(2^(2/3)*sqrt(3)*%i+2^(2/3))/4],
[x = -(2^(1/3)*sqrt(3)*%i+2^(1/3))/2,
y = (2^(2/3)*sqrt(3)*%i-2^(2/3))/4],[x = 2^(1/3),y = 2^(2/3)/2],
[x = 0,y = 0]] //

Not very elegant, but it works.  Numerically,

float(%);
==> [[x = 0.5*(2.182247271943443*%i-1.259921049894873),
y = -0.25*(2.749459273997205*%i+1.587401051968199)],
[x = -0.5*(2.182247271943443*%i+1.259921049894873),
y = 0.25*(2.749459273997205*%i-1.587401051968199)],
[x = 1.259921049894873,y = .7937005259840997],[x = 0.0,y = 0.0]] //

this is the same as what you got from solve (or algsys):

solve([a,b], [x,y]);
==> [[x = 1.259921095381759,y = .7937004973291583],
[x = 1.091123635971721*%i-.6299605249474366,
y = -.6873648184993013*%i-.3968502629920498],
[x = -1.091123635971721*%i-.6299605249474366,
y = .6873648184993013*%i-.3968502629920498],[x = 0,y = 0]] //

HTH,

Albert.