Lobachevsky Method

Back to MAIN

In the year 1834 Nikolaj Ivanovich Lobachevsky (1.Dec. 1792 - 24.Febr.1856) found a way to approximate the
roots of algebraic equations. This method has been discovered independently in the year 1837 by the German
mathematician Karl Heinrich Graeffe (7.Nov.1799 - 2.Dec.1873) and the French mathematician Germinal
Pierre Dandelin (12.Apr.1794 - 15.Febr.1847). We can see it today as a very convenient numerical method.

Before we start here, we should be aware of the fact that the equation:P(x) = xⁿ+a₁x^n-1+a₂x^n-2+ ... +a_n=0,
with the roots r₁,r₂, ... ,r_n and a given symetrical polynomial F[x₁,x₂,...,x_n], we can determine the number
s = F[r₁,r₂,...,r_n] without solving our original equation P(x) = 0.

As a consequence - in case of a given equation: P₁(x)=a₀xⁿ+a₁x^n-1+a₂x^n-2+ ... +a_n=0, and its roots
r₁,r₂, ... ,r_n we can find(without solving P₁) an other polynomial P_k(x)=a^k₀(x-r^k₁)(x-r^k₂)...(x-r^k_n), eg which has
roots which are exactly the squares, cubes etc. of the roots of P₁(x).

Just as an example, if the roots of the (above given)polynomial P(x) are x₁,x₂,x₃,x₄,x₅ then we can express the
sum of their 5^th powers as:(under the given index set {1,2,3,4,5})

S[x⁵_i] := x⁵₁+x⁵₂+x⁵₃+x⁵₄+x⁵₅ ; following this notation, we can now deduct:

S[x⁴_i]S[x_i] = S[x⁵_i] + S[[x⁴_ix_j];i < j] ==> S[x⁵_i] = S[x⁴_i]S[x_i]-S[x⁴_ix_j;i < j]= -a₁S[x⁴_i]-S[x⁴_ix_j;i < j], but
S[x³_i]S[x_ix_j;i < j] = S[x⁴_ix_j;i < j] + S[x³_ix_jx_k;i < j < k] ==> S[x⁴_ix_j;i < j] = a₂S[x³_i] - S[x³_ix_jx_k;i < j < k], but
S[x²_i]S[x_ix_jx_k;i < j < k] = S[x³_ix_jx_k;i < j < k] + S[x²_ix_jx_kx_l;i < j < k < l] .. which means that
S[x³_ix_jx_k;i < j < k] = -a₃S[x²_i] - S[x²_ix_jx_kx_l;i < j < k < l] .. and finally
S[x_i]S[x_ix_jx_kx_k;i < j < k < l] = S[x²_ix_jx_kx_l;i < j < j < l] + S[x_ix_jx_kx_lx_m;i < j < k < l < m] ..and so
S[x²x_ix_jx_kx_l;i < j < j < l] = a₄S[x_i] - 5a⁵ summarized as :

S[x⁵_i] = -a₁S[x⁴_i]-a₂S[x³_i]-a₃S[x²_i]-a₄S[x_i]-5a₅ and by setting q_k:=a^k₁+a^k₂+...+a^k_n we can reformulate it
as a general statement:

If the a_i are elementary symetrical functions of the variables x₁,x₂,...,x_n, then for all k < n+1 holds: q_k+a₁q_k-1+...+a_k-1q₁+ka_k = 0

..and in the same way we would be able to show that if k > n

q_k+a₁q_k-1+a₂q_k-2...+a_nq_k-n= 0

These both equations we know as Newton's recursion formula's . In the case n=4 we would obtain:

q₁+a₁=0, q₂+a₁q₁+2a₂=0, q₃+a₁q₂+a₂q₁+3a₃=0, q₄+a₁q₃+a₂q₂+a₃q₁+4a₄=0, q₅+a₁q₄+a₂q₃+a₃q₂+a₄q₁=0
q₆+a₁q₅+a₂q₄+a₃q₃+a₄q₂=0, q₇+a₁q₆+a₂q₅+a₃q₄+a₄q₃=0, e.t.c. and consequently..

q₁=-a₁, q₂= a²₁-2a₂, q₃=-a³₁+3a₁a₂-3a₃ q₄= a⁴₁-4a²₁a₂+4a₁a₃+2a²₂-4a₄
q₅=-a⁵₁+5a³₁a₂-5a²₁a₃-5a₁a²₂+5a₁a₄+5a₂a₃ e.t.c.

A small sparkling spin-off from all of this:

if we have P₁(x)=a₀xⁿ+a₁x^n-1+a₂x^n-2+ ... +a_n=0, and our roots r₁,r₂,...,r_n, then r²₁+r²₂+...+r²_n=a²₁-2a₂
which means that if a²₁ <= 2a₂, then at least one pair of roots must be complex.(we assume also a_n not 0,
and n > 2)

One important remaining point here:

Lets once look again at our P₁(x)=a₀xⁿ+a₁x^n-1+a₂x^n-2+ ... +a_n=0, with the roots r₁,r₂,...,r_n, and that means :
P₁(x)=a₀(x-r₁)(x-r₂) ... (x-r_n), and therefore P₁(-x)=(-1)ⁿa₀(x+r₁)(x+r₂) ... (x+r_n). What we look for is
a polynomial Q₂(x) having r²₁,r²₂,...,r²_n as roots - without solving P₁(x)=0, of course.
But (-1)ⁿP₁(x)P₁(-x)=a²₀(x²-r²₁)(x²-r²₂) ... (x²-r²_n)=Q₂(z) with z=x², and by decomposing P₁(x)=H₁(x)+xH₂(x),
where H₁(x) consists of even exponents of x, and xH₂(x) of the odd ones, we can conclude that
P₁(-x)=H₁(x)-xH₂(x), eg:

Q₂(x)=(-1)ⁿ[H²₁(x)-x²H²₂(x)]

But now we are better prepared for our little rendez-vous with an idea which now celebrates its 170^th birthday.

As a nature of the thing, we can ask us the question, if the polynomials Q₂,Q₄,Q₈,Q₁₆..(having the 2ⁿth power
of the original polynomial as roots) can help us to find r₁,r₂,...,r_n in a numerical way. The answer to that
question is yes. Lets assume, we know all r_i and lets for a moment also assume they are all real and different,
so after re-indexing we have: |r₁| > |r₂| > |r₃| > ... > |r_n-1| > |r_n| .. so if we take Q_2^k = a_0,kxⁿ+a_1,kx^n-1+...+a_n,k=0

and the known relations we see that:

-a_1,k/a_0,k = r^{2^k}₁+r^{2^k}₂+ ... +r^{2^k}_n = r^{2^k}₁[1+(r^{2^k}₂/r^{2^k}₁)+...+(r^{2^k}_n/r^{2^k}₁)]
+a_2,k/a_0,k = r^{2^k}₁r^{2^k}₂+r^{2^k}₁r^{2^k}₃ + ... + r^{2^k}_n-1r^{2^k}_n = r^{2^k}₁r^{2^k}₂(1+r^{2^k}₃/r^{2^k}₂ + ... + r^{2^k}_n-1r^{2^k}_n/r^{2^k}₁r^{2^k}₂)
-a_3,k/a_0,k = r^{2^k}₁r^{2^k}₂r^{2^k}₃+r^{2^k}₁r^{2^k}₂r^{2^k}₄ + ... + r^{2^k}_n-2r^{2^k}_n-1r^{2^k}_n = r^{2^k}₁r^{2^k}₂r^{2^k}₃(1+(r^{2^k}₄/r^{2^k}₃) + ... + (r^{2^k}_n-2r^{2^k}_n-1r^{2^k}_n/ r^{2^k}₁r^{2^k}₂r^{2^k}₃))

....
..etc..
....

(-1)ⁿa_n,k/a_0,k = r^{2^k}₁r^{2^k}₂ ... r^{2^k}_n-1r^{2^k}_n

.. with other words...

-a_1,k/(a_0,kr^{2^k}₁)=1+(r^{2^k}₂/r^{2^k}₁)+...+(r^{2^k}_n/r^{2^k}₁)
-a_2,k/(a_1,kr^{2^k}₂)=(1+r^{2^k}₃/r^{2^k}₂ + ... + r^{2^k}_n-1r^{2^k}_n/r^{2^k}₁r^{2^k}₂) / 1+(r^{2^k}₂/r^{2^k}₁)+...+(r^{2^k}_n/r^{2^k}₁)
-a_3,k/(a_2,kr^{2^k}₃)=(1+(r^{2^k}₄/r^{2^k}₃) + ... + (r^{2^k}_n-2r^{2^k}_n-1r^{2^k}_n/ r^{2^k}₁r^{2^k}₂r^{2^k}₃))/1+r^{2^k}₃/r^{2^k}₂ + ... + r^{2^k}_n-1r^{2^k}_n/r^{2^k}₁r^{2^k}₂

..etc

-a_n,k/(a_n-1,kr^{2^k}_n)=1/1+(r^{2^k}₂/r^{2^k}₁)+...+(r^{2^k}_n/r^{2^k}₁)

But as |r^{2^k}_i+1/r^{2^k}_i| < 1, all -a_i+1,k/(a_i,kr^{2^k}_i+1) must converge to 1 if i --> infinity. That means for "larger" values k,
we have

r^{2^k}_i+1 ~ -a_i+1,k/a_i,k

It is absolutely clear, that we could have taken any exponent - not just only powers of 2, we took it for practical
reasons

One possible solution using k4 can be:

/some preliminaries

mod:{x-y xbar x}
dix:{x@&~(!#x) in y}
f:{dix[x;&mod[(!#x);2]]}

xexp:{exp y*log x}
xe2:{xexp[2f;x]}
abs:{x|-x}
rot:{,/|(0;mod[x;#y])_y}

q:{+/rot'[-!#x;(,/'(s#'x) *\:+,x),\:(s+1-2*s:#x)#0]}

/lobachevsky .. just the correct sign needs to be verified after run .. a trivial K exercise
lb:{xexp[(abs@ - %/(1,-1) _\: x {f@ (q@ 1f*x*~b) - q@ 1f*x*b:mod . ((!#x);2)}/1f*y);%xe2 x]}

Sample 1: x⁵-2x⁴-19x³+13x²+36x-5 .. lb . (4;1 -2 -19 13 36 -5) returns:
x₁=5.001294, x₂=3.548031, x₃= 1.675143, x₄= 1.258001, x₅=0.133711
..the exakt solutions are (6 digits)
x₁=5.000000, x₂= -3.548948, x₃=1.674056, x₄= -1.258819, x₅=0.133711

..we see here that the polynomials having the 16^th power as roots of the original polynomial(eg.. we handle
here polynomials of degree 80) bring us actually very close to the real world..as we said, just ignoring the sign.
The precision (although very convincing in a normal case) is not the real strength of this method. The advantage
of the brilliant Lobachevsky approach is clearly, that we get all solutions in one go.

By definition, the assumption of strictly decreasing absolute values of the roots excludes the search for complex
solutions

If we are interested in the search for the complex roots - we slightly change our approach .. Just as an
illustrative display for a cubic polynomial(..and real coefficients), having one pair of complex solutions:

We assume that, r₁=Re^it and r₂=Re^-it and |r₁|=|r₂|>|r₃|,
so: r^k₁+r^k₂+r^k₃=R^k(e^ikt+e^-ikt)+r^k₃=R^k(2cos kt + r^k₃R^-k)=-a_1,k/a_0,k
and -a_2,k/a_0,k=r^k₁r^k₂+r^k₁r^k₃+r^k₂r^k₃=R^2k[1+2(r₃/R)^kcos kt] and -a_3,k/a_0,k = r^k₁r^k₂r^k₃ = R^2kr^k₃

(.. in case of |r₁|=|r₂|<|r₃| we would take -a_1,k/a_0,k=r^k₃[1+2(R/r₃)^kcos kt]
and -a_2,k/a_0,k=r^k₃R^k(2cos kt+(R/r₃)^k] and -a_3,k/a_0,k = R^2kr^k₃..)

So..-a_3,k/a_2,k ~ r^k₃ as (r₃/R)^k goes to 0 if k goes to infinity, and -a_1,k/a_0,k ~ r^k₃, and -a_2,k/a_0,k~R^2k

Sample 2: P(x) = x³+4x²-2x-20=0, by taking Q₁₆(z)=z³-8.446003*10⁷z²+1.000553*10¹⁶z-6.5536*10²⁰
, which means R~3.162332 and r₃~1.999931,(the correct figures are r₃=2 and R=3.162278)
Based on the knowledge of r₃=1.999931 and that r₁+r₂+r₃=-4 follows...that cos t = 0.5R^-1(-5.999931)
So, we obtain (as sign is not determined) t = 2.819756 resp. t = 0.3218366 ... a verification would tell us,
that the first solution is the correct one., So r₁=-3+i and r₂=-3-i.

One friendly hint here: in case of cubic polynomial with only real solutions, all polynomials Q_2^k(x)
must have positive roots, which means their coefficients must have the sign pattern ( + - + - ), which
is not the case in our example, consequently our cubic polynomial must have complex solutions.
It is by definition clear, that polynomials with real solutions having the same absolute value cannot be
handled by this method but those are just exceptional cases.

Back to MAIN