Bill Daly on Thu, 27 May 1999 11:59:44 -0400 |
[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]
Problem with symbolic calculation |
I tried the following GP program. The comments should explain the problem. /* I want to solve (symbolically) the problem of defining a conic y^2 + r*x*y + s*x^2 + t*y + u*x + v = 0 passing through a given set of 5 points [x1,y1] through [x5,y5]. This is a straightforward linear algebra problem which PARI should be able to handle. I start by defining a matrix m whose i-th row is [xi*yi,xi^2,yi,xi,1] and a column vector v = [-y1^2,-y2^2,-y3^2,-y4^2,-y5^2]~. */ m = matid(5); v = vectorv(5,j,0); { for (j = 1, 5, sx = eval(Str("x"j)); sy = eval(Str("y"j)); m[j,1] = sx*sy; m[j,2] = sx^2; m[j,3] = sy; m[j,4] = sx; m[j,5] = 1; v[j] = -sy^2; ); } /* At this point, I should be able to solve for [r,s,t,u,v]~ by computing m^-1*v, but when I tried to invert m, PARI failed to give an answer even with a 64 Meg stack. This doesn't look like such a difficult problem, so I tried to solve the system of equations longhand. */ { forstep (j = 5, 2, -1, for (k = 1, j-1, c = m[k,j]/m[j,j]; v[k] -= v[j]*c; m[k,] -= m[j,]*c; ); ); } /* The above loop should work, but PARI gives a divide-by-zero error in calculating c when k = 1 and j = 2. At this point, we have m[1,2] = x1^2 + ((((y4 - y5)/(x4 - x5))*x3^2 + (((-x4^2 + x5^2)/(x4 - x5)) *y3 + ((y5*x4^2 - x5^2*y4)/(x4 - x5))))/(((-y4 + y5)/(x4 - x5))*x3 + (y3 + ((-y5*x4 + x5*y4)/(x4 - x5)))))*x1 + (((-x3^2 + ((x4^2 - x5^2)/ (x4 - x5))*x3 + ((-x5*x4^2 + x5^2*x4)/(x4 - x5)))/(((-y4 + y5)/ (x4 - x5))*x3 + (y3 + ((-y5*x4 + x5*y4)/(x4 - x5)))))*y1 + ((((y5*x4 - x5*y4)/(x4 - x5))*x3^2 + ((-y5*x4^2 + x5^2*y4)/(x4 - x5))*x3 + ((x5*x4^2 - x5^2*x4)/(x4 - x5))*y3)/(((-y4 + y5)/(x4 - x5))*x3 + (y3 + ((-y5*x4 + x5*y4)/(x4 - x5)))))) m[2,2] = ((((y4 - y5)/(x4 - x5))*x2 + (-y2 + ((y5*x4 - x5*y4)/(x4 - x5)))) *x3^2 + (((-y4 + y5)/(x4 - x5))*x2^2 + (((x4^2 - x5^2)/(x4 - x5))*y2 + ((-y5*x4^2 + x5^2*y4)/(x4 - x5))))*x3 + ((x2^2 + ((-x4^2 + x5^2)/ (x4 - x5))*x2 + ((x5*x4^2 - x5^2*x4)/(x4 - x5)))*y3 + (((-y5*x4 + x5*y4)/(x4 - x5))*x2^2 + ((y5*x4^2 - x5^2*y4)/(x4 - x5)) *x2 + ((-x5*x4^2 + x5^2*x4)/(x4 - x5))*y2)))/(((-y4 + y5)/(x4 - x5)) *x3 + (y3 + ((-y5*x4 + x5*y4)/(x4 - x5)))) and m[1,2]/m[2,2] gives a divide-by-zero error. This can't be right. */ Regards, Bill ********************************************************************** This email and any files transmitted with it are confidential and intended solely for the use of the individual or entity to whom they are addressed. This footnote also confirms that this email message has been swept by MIMEsweeper for the presence of computer viruses. **********************************************************************