Path: news.mathworks.com!not-for-mail
From: "Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid>
Newsgroups: comp.soft-sys.matlab
Subject: Re: algebraic problem, no explicit solution
Date: Tue, 19 Feb 2008 19:10:22 +0000 (UTC)
Organization: The MathWorks, Inc.
Lines: 75
Message-ID: <fpf9iu$mh5$1@fred.mathworks.com>
References: <fpch15$b8o$1@fred.mathworks.com> <fpde6t$a55$1@fred.mathworks.com> <fpdpvq$133$1@fred.mathworks.com> <fpetpf$8i8$1@fred.mathworks.com>
Reply-To: "Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid>
NNTP-Posting-Host: webapp-05-blr.mathworks.com
Content-Type: text/plain; charset="ISO-8859-1"
Content-Transfer-Encoding: 8bit
X-Trace: fred.mathworks.com 1203448222 23077 172.30.248.35 (19 Feb 2008 19:10:22 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Tue, 19 Feb 2008 19:10:22 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 1187260
Xref: news.mathworks.com comp.soft-sys.matlab:452478


"Michael Snyder" <sirzerp@mathworks.com> wrote in message <fpetpf$8i8
$1@fred.mathworks.com>...
> Thank you Mr. Stafford.
> 
> I'm working on it today.  Yes I am using CAS.  Both Matlab
> and Maple.
> 
> Want to see something funny?  Type this command into Maple.
> 
> solvefor[y](1/k=1/((x-x1)^2+(y-y1)^2+(z-z1)^2)^(1/2)-1/((x-x2)^2+(y-
y2)^2+(z-z2)^2)^(1/2));
> 
> y =
> (-2*x*x1+x1^2+y1^2-2*z*z1+z1^2+2*x*x2-x2^2-y2^2+2*z*z2-z2^2
+2*k*Root\
> Of(4*x^2*y2^2-4*y1*z2^2*y2-4*x*x1*z1^2
+4*y1^2*x^2-8*x^2*y2*y1-8*y1*z^2*y\
> 2-4*y1*x2^2*y2-4*z1^2*y2*y1+4*x^2*x2^2-8*z*z1*x*x2
+2*x1^2*y2^2-4*x*x1*y1\
> ^2-4*x1^2*y2*y1-8*x*x1*z*z2+8*z*z1*y2*y1-4*_Z^7*k^3+
(-24*z1^3*k^2*z+12*x\
> ...
> 
> You will get 41 pages of result with 25 rootof causes.
> 
> Trying to solve a single line well known equation produces
> 41 pages!  That probably is why Matlab didn't return a result.  
> 
> BTW, I believe I want the complex roots.  I'm writing a
> paper, and some of it deals with the complex plane.    
> 
> Thanks again,
> 
> Michael Snyder
---------
  Michael, the whole problem becomes considerably easier if you rotate and 
translate your axes so that the midpoint of the line segment with endpoints 
(x1,y1,z1) and (x2,y2,z2) falls on the origin and the line segment is parallel to 
the new z axis.  (After all, any self-respecting dipole would want to be viewed 
in this esthetically pleasing symmetric manner.)  Then the problem can be 
expressed in cylindrical coordinates, z and r.  Your equation can be written

 1/sqrt((z-a)^2+r^2) - 1/sqrt((z+a)^2+r^2) = 1/k

and it is clear that we have axial symmetry about the z-axis.  The resulting 
eighth order polynomial has only even powers of z, so it can be solved as a 
fourth order polynomial in z^2.  Also the expressions for its five non-zero 
coefficients are considerably shorter.

  Here is the solution I got using my primitive Symbolic Toolbox (4a).

c8 = 1;
c6 = -4*(k^2-r^2+a^2);
c4 = 4*(a^2-3*r^2)*k^2+2*(3*a^4-2*a^2*r^2+3*r^4);
c2 = 16*a^2*k^4+4*(a^2-3*r^2)*(a^2+r^2)*k^2-4*(a^2-r^2)*(a^2+r^2)^2;
c0 = (a^2+r^2)^3*(a^2+r^2-4*k^2);
R = roots([c8,c6,c4,c2,c0]);
z = [sqrt(R);-sqrt(R)];

  As I mentioned earlier, not all of these eight roots will be solutions.  When I 
put in random values for a, r, and k, generally only two out of the eight z's 
were valid solutions, but as would be expected even these were not always 
real-valued.

  By the way, the solution you get using the 'collect' trick in the general case 
with arbitrary dipole orientation, though it is certainly complicated, has 
expressions for the nine coefficients that would all easily fit on one page.  
The entire polynomial occupies just 17 lines altogether, each of which is no 
more than 80 characters long.  I would say this method is far better than 
using 'solve'.  I would not want to have done this by hand, but the result is 
quite practical to carry out with matlab computations using 'roots'.

Roger Stafford