

Subject:
Find roots of the polynomial: a*x^6+b*x^5+c*x^4+d*x^3+e*x^2+f
Category: Science > Math Asked by: johnswanstonega List Price: $100.00 
Posted:
31 Dec 2004 22:30 PST
Expires: 30 Jan 2005 22:30 PST Question ID: 449910 
An algebraic solution is preferred as I'm going to have to solve this problem hundreds of thousands of times with different coefficients. Here are the coefficients: a = 4*n^4*(p^2+q^2) b = 4*n^2*p*(m^2*p^2+n^2*p^2+m^2*q^2+n^2*q^2+n^2*r^2) c = m^4*p^4n^4*p^42*m^2*n^2*p^4m^4*q^4n^4*q^4+ 2*m^2*n^2*q^42*m^4*p^2*q^22*n^4*p^2*q^2n^4*r^4 2*n^4*p^2*r^22*m^2*n^2*p^2*r^2+2*n^4*q^2*r^2+ 2*m^2*n^2*q^2*r^2 d = 4*m^2*n^2*p*r^2*(p^2+2*q^2) e = 2*m^2*p^2*r^2*(m^2*p^2+n^2*p^2+m^2*q^2+n^2*q^2+n^2*r^2) f = m^4*p^4*r^4 where m,n,p,q,r are arbitrary real number constants subject to the following constraints: n>0 m>0 r>0 r<p<r r<q<r p^2+q^2<r^2 Each time I solve this p and q will be different. r,n,m will change only every hundred thousand solutions or so. This is the reason I'm looking for an exact solution. Numerical methods will be too slow for the number of times I need to solve this. Alternately, each n,m,p,q,r I choose will reduce the coefficients to specific numbers. Then in my symbolic math software, I can use a "factor" command to solve each particular polynomial. I don't know what's going on behind the scenes, or if it is fast enough to do a hundred thousand of them in say 3 to 5 seconds on a modern PC, but if it is, then a description and implementation of the ?factor? algorithm in C or C++ would be acceptable as a solution to this problem. For those of you who are curious I will describe the origin of this problem: The problem I'm working on involves performing refraction distortion correction of the path of a laser beam as it passes through a human eye. The air/eye boundary is modeled in my problem as a hemisphere. In actuality it will be an arc of some radius and some extent which will be determined by the data. see: http://answers.google.com/answers/threadview?id=449830 The device can sample data at 256 different points (give or take) in time for each beam shot into the eye (using optical coherence tomography). This process is repeated 512 times (give or take), each time with the laser beam translated by a few 10's of micrometers across the eye. So I have this grid of data say 256 points deep and 512 points across. My job is create an image of the interior of the eye from this data. To do this I place a rectangular grid across the eye and then for each point <p,q> in the grid I ask what is the <x,y> value on the surface of the eye that would cause the laser beam to reflect into <p,q>. If I get lucky this will correspond to a single column in my data grid. Now using the distance between <p,q> and <x,y> I can find which data sample that corresponds to that <p,q> in the eye. If I'm unlucky (which I will be most of the time) I can find the pair of data columns that bracket the <x,y> point and interpolate the number I'm looking for from there. Ok, now here?s the analysis I?ve done to date: Let?s start with a hemisphere of radius r centered at the origin and extending from ?r to r on the xaxis and from 0 to r on the yaxis. Let?s assume that everything outside the hemisphere has an index of refraction n and everything inside has an index of refraction m. Let?s pick any point <x,y> on the hemisphere and any point <p,q> inside the hemisphere. The time it takes light to travel from the point <x,r> to <x,y> is: n/c * (ry) The time it takes light to travel from <x,y> to <p,q> is: m/c * sqrt((xp)^2 + (yq)^2)) And the total time is the sum of these two terms. Knowing that x^2+y^2=r^2 let?s us substitute and eliminate x in the prior equation. Now applying Fermat?s Principle, which is dt/dy = 0 let?s us differentiate the substituted equation of time with respect to y and set it?s value to 0. A bit of algebra leads us to the following equation: m*(p*y/xq) = n*sqrt(r^22*p*x+p^22*q*y+q^2) Squaring both sides and a bit of algebra leads to: [(m^2*p^2)/n^2]* y^2/x^2 ? [2*(m^2*p*q)/n^2] * y/x + [2*p] * x + [2*q] * y + [(m^2*q^2)/n^2]  p^2  q^2  r^2 = 0 We still know that: x^2 + y^2 = r^2 So we have 2 equation and 2 unknowns. Substituting y = sqrt(r^2 ? x^2) and squaring both sides again and quite a bit of algebra leads to the original 6th degree polynomial. For those you who are adventurous: Continuing with the original problem description ? Next I have to repeat this process for each boundary within the eye. The problem is that all boundaries other than the first need to be represented by splines (which are piecewise polynomials with the property of being differential in the 2nd degree at the piecewise boundaries). My thoughts are that I could solve for <x,y> on each polynomial similarly to what I?ve describe above for the hemisphere and exclude the solution if it is not in the "piece" of the spline that polynomial defines. Any thoughts? Suggestions? Ok. Thanks for getting this far. I hope to hear from you soon. J.  
 
 
 
 
 
 
 
 
 
 
 
 
 
 


Subject:
Re: Find roots of the polynomial: a*x^6+b*x^5+c*x^4+d*x^3+e*x^2+f
Answered By: mathtalkga on 30 Jan 2005 22:10 PST 
Hi, John: As you've outlined, a vertical beam of light strikes a point on an upper semicircle: x² + y² = r², y > 0 where the index of refraction changes from n above the semicircle to m below. According to Snell's Law, the light will deflect from a straight path toward the normal line when n < m and away from it when n > m. Given your application of correcting images of the eye, as light passes from air (index of refraction ~ 1) to the cornea, lens, and aqueous media beyond (index of refraction > 1), we will primarily treat the case where n < m. We will refer to ratio k = n/m rather than to the separate indices. By choosing units of length the radius of the semicircle can be assumed r = 1. Applying Fermat's principle you obtained for point (p,q) within the upper semicircle: py  qx  = k SQRT((px)² + (qy)²) x which we wish to solve for (x,y). Polynomial equation and extra roots  In deriving a polynomial from this and the equation x² + y² = 1, a number of "artifact" roots are introduced, ie. values x satisfying that polynomial formulation but not the original one. These extra roots are added by steps that multiply both sides by nonconstant polynomials or that square both sides of an equation. For example the equation x = 1 becomes x² = 1 when both sides are squared and acquires a "false root" x = 1 in the process. Nonetheless the desired root x is among those which satisfy: ax^6 + bx^5 + cx^4 + dx^3 + ex^2 + f = 0 where: a = 4(p² + q²)(k²)² b = 4((k² + 1)(p² + q² + 1)  1)k²p c = (k²  1)(q²)² + 2(((k²)² + 1)*(p²  1)  (k²  1))q² + ((k² + 1)p²)² + 2(k² + 1)k²p² + (k²)² d = 4k²p(p² + 2q²) e = 2((k² + 1)(p² + q² + 1)  1)p² f = (p²)² Note that while both odd and even powers of p appear among these coefficients, only even powers of q are present. This explains some of the artifact roots. Where a solution x for (p,q) in the upper semicircle exists, another solution x for mirrorimage point (p,q) in the lower semicircle may also exist. This accounts for an extra root which exists for the sample data: p = 0.5 q = 0.1 k = (1.0/1.5) = 2/3 While the refracted path through the point on the unit circle: x = 0.68024533659318476 y = 0.73298450327631404 passes through the point (p,q), this artifact solution: x = 0.75808380947988874 y = 0.65215714195618513 passes through (p,q). Similarly k only appears in even powers in the polynomial, and thus the other two real "solutions" actually satisfy for (p,q) or (p,q): py  qx  = k SQRT((px)² + (qy)²) x It is easy to see that these last artifact roots were introduced by squaring both sides of the Fermat's principle formulation. Impossibility of solution by radicals  Earlier in an RFC I reported the result of using MAXIMA and PARI/GP to determine if the particular sixthdegree polynomials produced in this application can be solved by radicals, ie. by algebraic formulas that require roots and fields operations acting on the coefficients of the polynomial. The result was that it is impossible to do so for the particular case p = 0.5, q = 0.1, k = 2/3 which you had used, and therefore impossible in general for these sixthdegree polynomials. However we did note that for the case q = 0, the resulting polynomial for x is only cubic and thus can be solved by radicals. Snell's Law formulation  A third formulation of the problem is possible using Snell's Law. By a bit of trigonometry: px  = tan(a  b) qy where a is the angle of vertical to the normal on the circle at (x,y) and b is the angle of refracted path to the same normal. Snell's Law says that: n * sin(a) = m * sin(b) and from the relations sin(a) = x, cos(a) = y, plus: sin(b) = kx cos(b) = SQRT(1  k²x²) one can rewrite tan(a  b) = sin(a  b)/cos(a  b) in terms of (x,y): p  x (1  k²)x  =  q  y y + SQRT(k²  (k²x)²) If we think about both sides as a difference quotient (slope formula) involving (x,y) and one other point, this can be interpreted as saying that (p,q) lies on the straight line connecting the two points: (x,y) and (k²x, SQRT(k²  (k²x)²)) Note that while (x,y) lies on the upper semicircle, the latter point falls into the lower halfplane because of the negative ycoordinate. It's possible then to determine for (x,y) on the upper semicircle the yintercept of the refracted path is (0,f), where: f = ( k²y + SQRT(k²  (k²x)²) )/(1  k²) This "focal" f has the limit of k/(1k) as points (x,y) tend to (0,1) and another limit k/SQRT(1  k²) as (x,y) tends to (1,0) or (1,0). Any point (x,y) on the upper semicircle will have its refracted path through (0,f) for some f intermediate between these limits. Newton iterations  Given a differentiable function F(x) and an initial guess x for a root of F(x) = 0, we compute successive Newton approximations: x := x  F(x)/F'(x) Any of the three formulations could be used as the basis for a Newton iteration. The polynomial form would make it easy in some respects to code the computation of the polynomial and its derivative together in a common loop over the coefficients. However the polynomial form has two drawbacks in this application. It has artifact roots that we'd need to be on guard against, and its coefficients require a greater amount of calculation from (p,q) each time they change than other formulations. The subject of starting values will be taken up in the next section. For the moment we will concern ourselves with the calculations needed for the Newton iteration itself. The Fermat principle can be reformulated in a variety of ways, even into the Snell's law form (and vice versa). So we don't claim that the particular choice of equation made here is a necessary one. Clearing denominators and moving all terms to one side produces: F(x) = py  (q + kd)x = 0 where d = SQRT((px)² + (qy)²). Differentiating with respect to x: F'(x) = py'  (q + kd)  kxd' where y' = x/y and d' = (p  x + (q  y)y')/d. Alternatively we can use the Snell's law formulation in a similar way, clearing the denominators and grouping all terms on one side: G(x) = (p  x)(y + z)  (1  k²)(q  y)x = 0 where z = SQRT(k²  (k²x)²). Differentiating with respect to x: G'(x) = (p  x)(y' + z')  (y + z)  (1  k²)(q  y  y'x) where y' is as above and z' = (k²)²x/z. Starting values  Once a Newton iteration is available for a problem that requires modest computation for an exact derivative, as here, the main difficulty left is to settle upon initial guesses that are good enough to allow the asymptotic quadraticconvergence to quickly appear. Here we have a couple of reasonable starting values that bracket the true root which are nearly "free" in terms of computation required. The point (p,SQRT(1  p²)) would lie on the straight path through (p,q) if k were 1, ie. if there were no refraction. On the other hand ( p/SQRT(p²+q²), q/SQRT(p²+q²) ) would lie on the line through (p,q) and the origin, which is the limiting path (the normal line) if k were 0. It is easy to show that the functions F(x) and G(x) both have a change in sign between these two points on the upper semicircle, thus proving that a root exists between them in each case. A simple linear interpolation between these values may be the best choice, certainly in the case that (p,q) is close to either the semicircle or to the yaxis. One can construct narrower bounds using a different pair of lines passing through (p,q), at a cost of some additional computation. For example, if we take the limiting points (0,f) found earlier: f = k/(1  k) and k/SQRT(1  k²) then lines from these two yintercepts through (p,q) will intersect the unit semicircle in a pair of points that more tightly bracket the root we want. Another thought here is that instead of using linear interpolation between the respective values of F(x) or G(x) at such a bracketing pair to obtain an initial guess for the Newton iteration, we can fairly easily find the intersection of the two refracted paths from those two points in the lower halfplane, then project from that "false focus" through (p,q) to a "guess" on the unit semicircle. regards, mathtalkga  
 
 
 
 


Subject:
Re: Find roots of the polynomial: a*x^6+b*x^5+c*x^4+d*x^3+e*x^2+f
From: mathtalkga on 03 Jan 2005 09:14 PST 
Progress Report 1. The intuitive notion that (x,y) on the circle of radius r centered at the origin will depend continuously on the parameters (p,q), a point in the interior of the disk of that radius, leads to a mapping of the open disk into the circle. Such a mapping cannot be both continuous and onto because of topological considerations. 2. Two simplifications that can be readily introduced are to eliminate parameter r (by scaling x,y,p,q accordingly) and to combine m,n into a ratio: P = p/r, Q = q/r, X = x/r, Y = y/r Then X = cos(t), Y = sin(t) and: P*cot(t)  Q = (n/m) SQRT( (cos(t)  P)^2 + (sin(t)  Q)^2 ) regards, mathtalkga 
Subject:
Re: Find roots of the polynomial: a*x^6+b*x^5+c*x^4+d*x^3+e*x^2+f
From: mathtalkga on 03 Jan 2005 13:11 PST 
Ooops, I got the tangent confused with the cotangent in that last equation, and it should read: P*tan(t)  Q = (n/m) SQRT( (cos(t)  P)^2 + (sin(t)  Q)^2 ) since y/x = sin(t)/cos(t) = tan(t). regards, mathtalkga 
Subject:
Re: Find roots of the polynomial: a*x^6+b*x^5+c*x^4+d*x^3+e*x^2+f
From: johnswanstonega on 03 Jan 2005 17:00 PST 
Mathtalk, The analysis I've posted here deals with minimizing the time it takes a ray of light to travel from (x,r) to (p,q) according to Fermat's Principle and it leads to 6th order polynomials no matter which variable you solve for; x or y. Another analysis would be a geometrical one as follows: The equation of a circle in terms of an angle, a: x = r*sin(a) y = r*sin(b) Using the rule of similar angles of a parallel lines we know that a is also the angle that a ray of light traveling parallel to the yaxis will strike the surface of the semicircle. Using Snell's law: n*sin(a) = m*sin(b), where b will be the angle the light will refract at the boundary of the semicircle. Using the two point equation of a line through <x,y> and <p,q>: px = s * (qy) we need to find the slope s which in this case will be express as run/rise. Using geometrical reasoning it is known that s = tan(a  b). So dividing the line equation by (qy) gives: px  = s qy Substituting for 'x', 'y', and 's' gives: pr*sin(a)  = tan(a  b) qr*cos(a) solving Snell's law for b: b = arcsin(sin(a)/n) a final substitution gives: pr*sin(a)  = tan(a  arcsin(sin(a)/n) qr*cos(a) Which if solved for 'a' would also give me the solution I'm looking for. If it is impossible to solve this equation or the polynomial analytically then I'd accept a solution which describes and implements the "solve by factoring" algorithm used in symbolic algebra packages such as Mathmatica, Maple, Maxima, or MacSym. When I execute this function the computer comes back with all the roots of the polynomial including the complex roots. What I'd do then is use this to find the root that is the solution for my initial <p,q> point and then use the computed <x,y> as the initial guess for the next <p+delta, q+delta> as per your earlier suggestion. I would also accept an analysis of this equation or the polynomial similar to the one in this paper on page 19 (http://www.math.uab.edu/cl/cl1/details.ps) which guarantees that a numerical method will always converge to the right solution given a known starting value. Let me know if there's anything more I can do to help you. Thanks, John. 
Subject:
Re: Find roots of the polynomial: a*x^6+b*x^5+c*x^4+d*x^3+e*x^2+f
From: johnswanstonega on 03 Jan 2005 17:19 PST 
Mathtalk, Addressing your question about what happens when <p,q> = <0,0> Under this condition, the time it takes light to travel that path is reduced to: t = n/c * (ry) + m/c * r My original analysis minimized this function by taking the derivitave with respect to y and setting it to zero, then solving for y. In this case that approach doesn't work since taking the derivitive leads to dt/dy = n/c which is a constant. (Remember c here represents the speed of light in air). I don't know how to minimize this function ... maybe you do. John. 
Subject:
Re: Find roots of the polynomial: a*x^6+b*x^5+c*x^4+d*x^3+e*x^2+f
From: mathtalkga on 03 Jan 2005 18:16 PST 
To minimize: t = n/c * (ry) + m/c * r with respect to y on the interval [r,+r], take y = r. Here the minimum occurs at a boundary point, rather than at an interior "critical" point such as you might find by setting the derivative to zero. I have a minor correction to your alternative derivation; regarding Snell's Law: n*sin(a) = m*sin(b) we obtain: b = arcsin( (n/m) * sin(a) ) This suggests a natural limitation on angle a (equivalent to my angle t), that:  sin(a)  <= m/n in order that angle b is defined by this relation. regards, mathtalkga 
Subject:
Re: Find roots of the polynomial: a*x^6+b*x^5+c*x^4+d*x^3+e*x^2+f
From: mathtalkga on 04 Jan 2005 05:09 PST 
Hi, John: It appears that your problem has a left/right symmetry, so that: (x,y) the solution for (p,q) implies: (x,y) the solution for (p,q) but is not up/down symmetric in the corresponding sense. However the formulation: m*(p*y/xq) = n*sqrt(r^22*p*x+p^22*q*y+q^2) implies that p*(y/x)  q is nonnegative, so y/x >= q/p when p > 0. regards, mathtalkga 
Subject:
Re: Find roots of the polynomial: a*x^6+b*x^5+c*x^4+d*x^3+e*x^2+f
From: johnswanstonega on 04 Jan 2005 09:46 PST 
Mathtalk, Going back to the original problem: t = n/c * (ry) + m/c * sqrt((xp)^2+(yq)^2) and x^2 + y^2 = r^2 you can either substitute for y or x. Substituting for x gives an equation in y: t = n/c * (ry) + m/c * sqrt((sqrt(y^2+r^2)p)^2+(q+y)^2) Substituing for y gives and equation in x: t = n/c * (rsqrt(x^2+r^2)) + m/c * sqrt((xp)^2+(sqrt(r^2x^2)q)^2) My problem is to minimize either one (or both) of these equations. i.e. Take the derivitive with respect to the remaining variable, set it equal to zero, and solve for the roots. Doing it either way leads to a 6th order polynomial. Solving the equation for x leads to a polynomial without the x^1 term. Solving the equation for y leads to a polynomial with all the terms. This equation: m*(p*y/xq) = n*sqrt(r^22*p*x+p^22*q*y+q^2) is simply the derivitive of the equation in x set to 0 with y substituted back in. This was my first attempt and it turned out to be quite messy. I mention this because it seems like your focusing in on this equation and while solving it would answer my question, I wanted to point out the other equations as well since they might be easier to work with. John. 
Subject:
Re: Find roots of the polynomial: a*x^6+b*x^5+c*x^4+d*x^3+e*x^2+f
From: mathtalkga on 04 Jan 2005 10:21 PST 
Thanks, John. I'm posting a mixture of Comments here, both for the sake of giving a "progress" indicator and for the sake of your feedback. My understanding of the physics of the application is that, given a point (p,q) inside a "crystal ball" you want to know where to shine a light to illuminate that point. Outside the ball and inside the ball the light will travel in a straight line, but at the boundary (x,y) the light beam is refracted due to the changing speed of light. I can show a solution exists for p,q > 0, and by left/right symmetry we can drop the restriction that p > 0. However I'm expecting that there are values q < 0 which have no solution. The physical intuition is that as a light beam approach the vertical edge x = r, there may be points (p,q) "below" the path inside the ball that cannot be attained. regards, mathtalkga 
Subject:
Re: Find roots of the polynomial: a*x^6+b*x^5+c*x^4+d*x^3+e*x^2+f
From: mathtalkga on 05 Jan 2005 07:33 PST 
I haven't been able to show directly the equivalence of the two formulations, but they share a characteristic which I think can be used to show the equivalence indirectly. This is that (p,q) may be replaced in either equation by any point inside the circle on a line between (x,y) and (p,q). The Snell's Law formulation: (px)/(qy) = tan(a  b) isolates all the dependence on (p,q) on one side of the equation, as the angles a,b depend only on (x,y) and the ratio of refraction indexes. So adding a multiple of (px) to p and the same multiple of (qy) to q will leave the equation unchanged. The situation is only a little more complicated with the "original" formulation: p*(y/x)  q = (n/m)*SQRT( (px)^2 + (qy)^2 ) but the corresponding adjustments to (p,q) would proportionally change both sides of the equation. An indirect verification of equivalence can therefore be derived by showing that the slopes of the lines as (p,q) approaches (x,y) are the same for both equations, e.g. dp/dq = tan(a  b). * * * * * * * * * * * * * * * Closer to your original request for an exact solution, I observe that for the special case q = 0, a cubic equation can be formulated for x with coefficients in p, r, and k = (n/m): p*(y/x) = k*SQRT( p^2  2px + r^2 ) After multiplying both sides by x and squaring: (p^2)*(r^2  x^2) = (k^2)*(x^2)*(p^2  2px + r^2) Therefore an "exact solution" for x can be given in this case using Cardan's formula for roots of a cubic equation. See for example: http://answers.google.com/answers/threadview?id=433886 http://answers.google.com/answers/threadview?id=441538 * * * * * * * * * * * * * * * There would be a very simple exact solution if all of the lines (p,q) to (x,y) described in the first section of this Comment were to intersect at a common focus, say (0,f). For then one has only to intersect a line through (p,q) and (0,f) with the upper semicircle x^2 + y^2 = r^2 to find (x,y). regards, mathtalkga 
Subject:
Re: Find roots of the polynomial: a*x^6+b*x^5+c*x^4+d*x^3+e*x^2+f
From: johnswanstonega on 05 Jan 2005 11:17 PST 
Mathtalk, There would be a very simple exact solution if all of the lines (p,q) to (x,y) described in the first section of this Comment were to intersect at a common focus, say (0,f). For then one has only to intersect a line through (p,q) and (0,f) with the upper semicircle x^2 + y^2 = r^2 to find (x,y). I agree that if all the lines (p,q) to (x,y) were to intersect at a common point then the solution would be easy. Unfortunatly this is not the case for general values of n and m. Your observation that "(p,q) may be replaced in either equation by any point inside the circle on a line between (x,y) and (p,q)." could be expoited to reduce the dimensionality of the problem if it could be reformulated using this knowledge (something like a coordinate transform or an equation using this property as the basis functions). This reduction in dimensionality may then enable an algebriac solution. What do you think? John. 
Subject:
Re: Find roots of the polynomial: a*x^6+b*x^5+c*x^4+d*x^3+e*x^2+f
From: mathtalkga on 05 Jan 2005 11:56 PST 
Hi, John: Yes, in a special case one can reduce the dimensionality to one (e.g q = 0) and obtain a reduced degree polynomial problem with a "closed form" solution. However I still expect that with good initial approximations, the numerical method will converge just as fast as any "exact" solution. I'm thinking of "projecting" points onto the xaxis (q=0) using an approximate focal point (p=0) mostly as a technique for generating the initial approximations. I believe that for a "spherical" lens, paraxial rays do have a fairly tight "focus"; it's the rays which fall near x = r which are responsible for the majority of "spherical aberration". I think there is a critical angle: y/x = n/m at which refraction changes to reflection that the math "mirrors" the physical limitation in your model. regards, mathtalkga 
Subject:
Re: Find roots of the polynomial: a*x^6+b*x^5+c*x^4+d*x^3+e*x^2+f
From: mathtalkga on 07 Jan 2005 07:27 PST 
Hi, John: Here's my approach, still to be coded and tested. I'm assuming n/m < 1 is the case of real interest. Given an array of (p,q) points: . . . . . . . . . . . . . . . . . . we solve (assuming these lie on one side of the yaxis, since solutions are symmetric about that line) for the extreme corners and midpoint, finding an approximate common focal point where the three lines thru respective (x,y) and (p,q) nearly intersect. Then use this approximate focal point (p',q') to produce initial approximations (x',y') for each additional solution desired. Two Newton iterations should suffice to give all necessary accuracy; three should converge to machine precision (digits of precision will double with each iteration), assuming a fairly tight angular distribution reasonably close to the central axis. See my Request for Clarification above about the realistic size of the angular range for the block of (p,q) points. regards, mathtalkga 
Subject:
Re: Find roots of the polynomial: a*x^6+b*x^5+c*x^4+d*x^3+e*x^2+f
From: mathtalkga on 09 Jan 2005 16:21 PST 
Hi, John: A couple of quick notes: I found a way to show that the Snell's Law and Fermat's principle formulations are equivalent for q < y. I coded and tested Newton iterations for both versions. They have about the same operation count, and for n/m < 1, both seem to converge satisfactorily within 2 or 3 iterations if given starting values of only modest accuracy. Based on the limited amount of testing done so far, the method based on Fermat's principle seems to converge a little faster. regards, mathtalkga 
Subject:
Re: Find roots of the polynomial: a*x^6+b*x^5+c*x^4+d*x^3+e*x^2+f
From: mathtalkga on 20 Jan 2005 14:20 PST 
The polynomial for x or y turns out to be sextic (of sixth degree), but it seems reasonable that by choosing a different unknown the problem may be restated in a form of lesser degree. One suggestion of this has already been seen when, by restricting q = 0, one obtains a formulation in which x satisfies a cubic polynomial (whose coefficients depend algebraically on p) and in which p satisfies a quadratic polynomial (whose coefficients depend algebraically on x). I'm especially interested in a formulation in terms of the yintercept (0,f) of the line passing through (p,q), rather than in terms of point (x,y) where the line intercepts the (unit) circle. It seems to me that the left/right symmetry of the problem may produce a reduction in the degree of the problem. For example, suppose that the coordinate f of the yintercept satisfies an even polynomial of degree 6 (whose coefficients depend algebraically on p,q). The square of f must then satisfy a corresponding cubic equation, so an explicit formula for f (in terms of root extractions and field operations) could in principle be written down. If k = n/m, the dependence of f (which I've written with a negative sign above because the light path incident on the upper semicircle will have yintercept in the lower halfplane) on x,y can be stated: __________ k²y + k? 1  k²x² f =  (1  k²) This gives values of f that range from: ____ f = k/?1k² when (x,y) = (1,0) to: f = k/(1k) when (x,y) = (0,1) Even if this line of investigation doesn't produce the sought after closedform solution, it may still be useful in providing an "approximate focal point". That is, if we take an approximate yintercept (0,f) and project through (p,q) to the point on the unit circle (x,y), we may get a very satisfactory starting value for a Newton iteration. regards, mathtalkga 
Subject:
Re: Find roots of the polynomial: a*x^6+b*x^5+c*x^4+d*x^3+e*x^2+f
From: johnswanstonega on 27 Jan 2005 22:00 PST 
Mathtalk, I agree with your "gut" feeling that if I were just smart enough to know how to reformulate the problem, I could get an algebriac solution. I'm thinking that your "yintercept" idea is a good one. A coordinate system based on a point on the yaxis and an angle from vertical would seem to exactly encapsulate the free variables in this problem assuming of course that n/m < 1 (i.e. all things bend to the yaxis). John. 
If you feel that you have found inappropriate content, please let us know by emailing us at answerssupport@google.com with the question ID listed above. Thank you. 
Search Google Answers for 
Google Home  Answers FAQ  Terms of Service  Privacy Policy 