MATLAB Central Newsreader - algebraic problem, no explicit solution

algebraic problem, no explicit solution

Michael Snyder — Mon, 18 Feb 2008 17:59:01 -0500

I have a well known physics equation (dipole field) and I’m
trying to solve for constant values. I can fill an array
with values and have Matlab plot the constant values very
easily. Wanted to solve the equation for a single constant
value and make it into a function that plots one curve at a
time.

i.e. k=.002 is a curved line and k=.03 is a different
curved line.

x1 y1 z1, x2 y2 z2 are fixed foci points.

k is a given value, I just want to solve for y in terms of x
and z.

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

>> syms x x1 x2 y y1 y2 z z1 z2 k K

>>
K=1/((x-x2)^2+(y-y2)^2+(z-z2)^2)^(1/2)-1/((x-x1)^2+(y-y1)^2+(z-z1)^2)^(1/2)-1/k

>> solve(K,y)

Warning: Explicit solution could not be found.

ans =

[ empty sym ]

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

How can a well known algebraic equation have no explicit
solution? It is frustrating because I know which parameters
plot which lines.

Also I’ve tried a log form and got the same results.

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

K =
log(((x-x2)^2+(y-y2)^2+(z-z2)^2)^(1/2)-((x-x1)^2+(y-y1)^2+(z-z1)^2)^(1/2))-1/2*log((x-x1)^2+(y-y1)^2+(z-z1)^2)-1/2*log((x-x2)^2+(y-y2)^2+(z-z2)^2)-log(1/k)

>> solve(K,y)

ans =

[ empty sym ]

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

Re: algebraic problem, no explicit solution

Roger Stafford — Tue, 19 Feb 2008 02:17:01 -0500

"Michael Snyder" <sirzerp@mathworks.com> wrote in message <fpch15$b8o
$1@fred.mathworks.com>...
> I have a well known physics equation (dipole field) and I’m
> trying to solve for constant values. I can fill an array
> with values and have Matlab plot the constant values very
> easily. Wanted to solve the equation for a single constant
> value and make it into a function that plots one curve at a
> time.
>
> i.e. k=.002 is a curved line and k=.03 is a different
> curved line.
>
> x1 y1 z1, x2 y2 z2 are fixed foci points.
>
> k is a given value, I just want to solve for y in terms of x
> and z.
>
> --------------------------------------
>
> >> syms x x1 x2 y y1 y2 z z1 z2 k K
>
> >>
> K=1/((x-x2)^2+(y-y2)^2+(z-z2)^2)^(1/2)-1/((x-x1)^2+(y-y1)^2+(z-z1)
^2)^(1/2)-1/k
>
> >> solve(K,y)
>
> Warning: Explicit solution could not be found.
>
> ans =
>
> [ empty sym ]
>
> --------------------------------------
>
> How can a well known algebraic equation have no explicit
> solution? It is frustrating because I know which parameters
> plot which lines.
>
> Also I’ve tried a log form and got the same results.
>
> --------------------------------------
>
> K =
> log(((x-x2)^2+(y-y2)^2+(z-z2)^2)^(1/2)-((x-x1)^2+(y-y1)^2+(z-z1)^2)^
(1/2))-1/2*log((x-x1)^2+(y-y1)^2+(z-z1)^2)-1/2*log((x-x2)^2+(y-y2)^2+
(z-z2)^2)-log(1/k)
>
> >> solve(K,y)
>
> ans =
>
> [ empty sym ]
-----------
  After playing with your equation a bit, it is clear that y can be expressed as
the solution to an eighth degree polynomial equation in terms of the other
quantities. As you probably know, there is no general solution to eighth
degree polynomial equations in terms of elementary functions, but, in theory,
the Symbolic Toolbox could have given the solution to you in terms of roots
of this polynomial. I don't know why it didn't succeed in this. My own
primitive version of the Symbolic Toolbox ran into a limitation of string
length since it is an early student version (1994.) In any event, it wouldn't
have given you the kind of explicit solution I think you were looking for
because in all likelihood such a solution does not exist. However, once you
derive this polynomial, you can use matlab to solve for y numerically using
the 'roots' function.

  I will get you started deriving the polynomial. First, it will simplify things if
you call

P = (x-x1)^2+(y-y1)^2+(z-z1)^2
Q = (x-x2)^2+(y-y2)^2+(z-z2)^2

This gives us

1/sqrt(P) - 1/sqrt(Q) = 1/k

1/sqrt(P) = 1/k + 1/sqrt(Q)

Square both sides:

1/P = 1/k^2 + 2/(k*sqrt(Q)) + 1/Q

2/(k*sqrt(Q)) = 1/P - 1/Q - 1/k^2

Again square both side and transpose:

1/k^4-2/(k^2*P)-2/(k^2*Q)+1/P^2-2/(P*Q)+1/Q^2 = 0

Multiply by k^4*P^2*Q^2:

P^2*Q^2 - 2*k^2*P*Q*(P+Q) + k^4*(P-Q)^2 = 0

If you now substitute the above expressions in y in place of P and Q, and
multiply it all out, you will see the eighth degree polynomial in y exhibited in
all its glory. I leave that to you.

  To solve it numerically, you need to collect like powers of y and use the
'roots' function for specific values of (x1,y1,z1), (x2,y2,z2), x, z, and k. In
each case there should be eight roots altogether, but not all of them satisfy
your original equation. It will be your problem to sort out the real-valued
solutions that do.

Roger Stafford

Re: algebraic problem, no explicit solution

Roger Stafford — Tue, 19 Feb 2008 05:38:02 -0500

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in
message <fpde6t$a55$1@fred.mathworks.com>...
> .....SNIP.....
> If you now substitute the above expressions in y in place of P and Q, and
> multiply it all out, you will see the eighth degree polynomial in y exhibited
> in all its glory. I leave that to you.
> .....SNIP.....
> Roger Stafford
--------
Michael, when I said, "I leave that to you", in reference to finding
expressions for each of the coefficients of that eighth-order polynomial, I
didn't realize how much labor it would involve doing it by hand. (You must
have thought I was rather cruel.) I would recommend using your Symbolic
Toolbox for that task instead. Use 'subs' to substitute in
(x-x1)^2+(y-y1)^2+(z-z1)^2 for P and (x-x2)^2+(y-y2)^2+(z-z2)^2 for Q,
or better still, substitute (y-y1)^2+A for P and (y-y2)^2+B for Q where A and
B stand for (x-x1)^2+(z-z1)^2 and (x-x2)^2+(z-z2)^2, respectively. Then
use 'collect(-,y)' and 'simple' to get expressions for the coefficients of the
polynomial in y from P^2*Q^2-2*k^2*P*Q*(P+Q)+k^4*(P-Q)^2. You can
expect these coefficient expressions to be rather long even with
simplification. You can use the resulting nine coefficient expressions for
substituting in numerical values so that 'roots' can find specific solutions to
the polynomial. Note that you will have to use the original equation

1/sqrt(P) - 1/sqrt(Q) = 1/k

to eliminate false solutions, as well as discarding the complex-valued ones.

Roger Stafford

Re: algebraic problem, no explicit solution

Michael Snyder — Tue, 19 Feb 2008 15:49:03 -0500

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

>
> to eliminate false solutions, as well as discarding the
complex-valued ones.
>
> Roger Stafford
>

Re: algebraic problem, no explicit solution

roberson@ibd.nrc-cnrc.gc.ca (Walter Roberson) — Tue, 19 Feb 2008 16:42:13 -0500

In article <fpetpf$8i8$1@fred.mathworks.com>,
Michael Snyder <sirzerp@mathworks.com> wrote:

>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));

"solvefor is deprecated. Please use the solve command"

>You will get 41 pages of result with 25 rootof causes.

In Maple 11, I get approximately 17200 characters of output.
If I put spaces around the '+' and automatically reformat,
it comes out about 300 lines. Only 4 RootOf clauses are
involved.
--
"The shallow murmur, but the deep are dumb." -- Sir Walter Raleigh

Re: algebraic problem, no explicit solution

Roger Stafford — Tue, 19 Feb 2008 19:10:22 -0500

"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