Program for Polynomial Reduction

This file may be copied and/or distributed.

This is one of three inter-related files. To locate the other two files, access Three Surface Intersection Problem.

The embedded data in this program listing corresponds to "Section Two" of the program output file which will be found by following the above link. To run other data cases, edit the a, b, c, d, e data below. Note that this very basic program coding does not provide for a data input file mechanism.

--------------------------

   
/* Three surface program Test Case 2.wpd    November 21, 2011 */
/* Author: David Ziskind, Marlborough, CT */
/* 11/21/11: Changed test on solution_error from 1.e-3 to 1.e-2 .*/

/* begin one time code */

a[1]:-1.03$   a[2]:2$        a[3]:-0.01$
b[1]:-1$      b[2]:0$        b[3]:0$
c[1]:1$       c[2]:0.01$     c[3]:0.04$
d[1]:1.01$    d[2]:0.02$     d[3]:2$
e1  :3.9$     e2  :4.1$      e3  :8.4$

x[1]:x1$  x[2]:x2$  x[3]:x3$

Da : sqrt( sum((x[m]-a[m])^2, m,1,3) ) $
Db : sqrt( sum((x[m]-b[m])^2, m,1,3) ) $
Dc : sqrt( sum((x[m]-c[m])^2, m,1,3) ) $
Dd : sqrt( sum((x[m]-d[m])^2, m,1,3) ) $

Dasqr : sum((x[m]-a[m])^2, m,1,3) $
Dbsqr : sum((x[m]-b[m])^2, m,1,3) $
Dcsqr : sum((x[m]-c[m])^2, m,1,3) $
Ddsqr : sum((x[m]-d[m])^2, m,1,3) $

ratprint : false$
Sbound   : 1.e-9$
case_cntr: 0$

/* Following equations have no part in the solution process but are the
basic equation set (before square operations).  These equations are supplied
to enable checking as to whether the solutions of the polynomial equations do
in fact satisfy the basic equation set. */

eqbasic1 : Db + Dc - e1 $
eqbasic2 : Da + Db - e2 $
eqbasic3 : Dc + Dd - e3 + e2 $ 

"Equations eq1, eq2, eq3 follow:";
eq1 : ratsimp((Dbsqr - e1^2 - Dcsqr)^2 - 4*e1^2*Dcsqr,  x3,x2,x1) ;
eq2 : ratsimp((Dbsqr - e2^2 - Dasqr)^2 - 4*e2^2*Dasqr,  x3,x2,x1) ;
eq3 : ratsimp((Dcsqr - (e3-e2)^2 - Ddsqr)^2 - 4*(e3-e2)^2*Ddsqr,  x3,x2,x1);

x1eqpp : part(eliminate([eq1,eq2,eq3], [x2,x3]), 1)$
x1eq   : ratexpand(x1eqpp)$
disp("x1eq follows", float(x1eq))$
x1roots : realroots(x1eq, Sbound)$
/* ========================== end one time code */

Fx1(i) := block(
x1p    : part(x1roots, i,2),
x1pflt : float(x1p),
x2eqpp : part(eliminate([ev(eq2,x1=x1p), ev(eq3,x1=x1p)], [x3]), 1),
x2eq   : ratexpand(x2eqpp),
disp("x2eq follows", float(x2eq)),
x2roots : realroots(x2eq, Sbound),
nmbr_x2roots : length(x2roots),
display(i, nmbr_x2roots),

Fx2(j) := block(
x2p     : part(x2roots, j,2),
x2pflt  : float(x2p),
x3eqpp  : part(ev(eq3,x1=x1p,x2=x2p), 1),
x3eq    : ratexpand(x3eqpp),
x3roots : realroots(x3eq, Sbound),

Fx3(k) := block(
case_cntr : case_cntr + 1,
x3p : part(x3roots, k,2),
x3pflt : float(x3p),
solution_error :
  float(sqrt(
      ev(eq1,x1=x1p,x2=x2p,x3=x3p)^2
    + ev(eq2,x1=x1p,x2=x2p,x3=x3p)^2
    + ev(eq3,x1=x1p,x2=x2p,x3=x3p)^2
    )),
if solution_error <= 1.e-2
then block(
basic_set_error :
  float(sqrt(
      ev(eqbasic1,x1=x1p,x2=x2p,x3=x3p)^2
    + ev(eqbasic2,x1=x1p,x2=x2p,x3=x3p)^2
    + ev(eqbasic3,x1=x1p,x2=x2p,x3=x3p)^2
    )),
display(i, j, k, x1pflt, x2pflt, x3pflt, solution_error, basic_set_error)
) /* end of if solution_error block */
), /* end of Fx3 block */ 

/* do loop for the level 3 loop */
for kk : 1 step 1 thru length(x3roots) do Fx3(kk)
), /* end Fx2 block */

/* do loop for the level 2 loop */
for jj : 1 step 1 thru length(x2roots) do Fx2(jj)
)$ /* end Fx1 block */

/* ================== begin one time code */
"Show solutions for the polynomial system";
for ii : 1 step 1 thru length(x1roots) do Fx1(ii)$
display(case_cntr, Sbound)$ 

"end of run";
  
Home

Rev. 11/26/11 g