<- previous index next ->
Finding Roots
First, a special case, the "roots of unity"
x +1=0 root is x=-1
x^2+1=0 roots are x=i x=-i i=sqrt(-1)
x^3+1=0 roots are x=-1 x=1/2 +i sqrt(3)/2 x=1/2 -i sqrt(3)/2
x^4+1=0 roots are x=+sqrt(2)/2 +i sqrt(2)/2
(all four combinations of +/-)
Note that all roots are on the unit circle in the complex plane.
Consider the Maple root finding for x^12+1=0 :
root12g.html
In general the roots of a polynomial are complex numbers.
A general purpose root finder for polynomials is difficult to
develop. Fortunately, some very smart people have written the
function, cpoly, that can take a general polynomial with complex
coefficients and find the complex roots. It also works with
real coefficients and real roots.
MatLab root finding for polynomials is just r = roots(p)
Where p is the vector of polynomial coefficients and
r is the vector of roots. Both r and p may be complex.
Finding roots of a polynomial seems easy:
given c_n x^n + c_n-1 x^n-1 + ... _ c_1 x + c_0 = 0
find a value, x=r, that satisfies the equation,
divide the polynomial by (x-r) and thus have a
polynomial of one degree lower. Repeat. Done.
Unfortunately there are a lot of pitfalls, solved by the codes below.
Standard techniques include dividing through by c_n,
if c_0 is zero, there is a root equal zero, divide through by x,
if highest power is 2, directly solve the quadratic equation.
The Fortran code, including the test program (driver) is
419.for in Fortran IV, compiles with g95
419_f.out is the output of many test cases
The Java code, including the test program (driver) is
c419.java in standard Java
c419.out is the output of many test cases
The Ada code, including the test program (driver) is
long_complex_polynomial_roots.ada in Ada
To understand how some iterative methods can converge quickly,
consider how to compute sqrt when you have no math library:
Basically y=sqrt(x) guess y, y_next = (x/y+y)/2, repeat.
sqrt.c
sqrt_c.out
A professional implementation of sqrt would reduce the input, x,
to a range from 0.25 to 1.0 by dividing by an even power of 2.
The sqrt(x^2n) is just x^n. The square root of a product is
the product of the square root of the factors. This can be easy
with IEEE floating point numbers because they are stored with
an exponent of 2.
So, since the sqrt iteration was basically Newtons method,
why not just use Newtons method to find roots?
here is how it works, even for complex roots.
f(x) = 0
guess x
x_next = x - f(x)/f'(x) f'(x) is the derivative of f(x)
repeat until |x_next-x| < epsilon
newton.c
newton_c.out
Notice the oscillation in the last few steps.
The problem can come from a bad guess that causes a big
oscillation. The problem can come from hitting an x
where f'(x) is zero. Thus, there have been many workarounds.
"cpoly" is one of the most general solutions.
Nonlinear Equations
Other examples: same equation with different constant.
First, expect x=1, we happen to have at least two possible solutions:
simple_nl.c
simple_nl_c.out
Second, expect x=2, have convergence to expected solution:
simple2_nl.c
simple2_nl_c.out
Systems of Nonlinear Equations
Three equations in three unknowns, not direct convergence.
Experiments show great difference with widely different
initial conditions. Two sets of equations with a choice of
initial condition. x1, x2, x3 are the variables.
equation_nl.c
equation_nl_c.out
inverse.c used by code above
In other languages:
equation_nl.adb
equation_nl_ada.out
inverse.adb
equation_nl.f90
equation_nl_f90.out
inverse.f90
equation_nl.java
equation_nl_java.out
inverse.java
A system of non linear equations with powers 1, 2, 3, -1 and -2.
This example has a lot of debug print and several checks.
Added is also, variable control of fctr, a multiplier on
how much correction is to be made each iteration.
The problem and method are described as comments in the code.
equation2_nl.adb
equation2_nl_ada.out
Same example as above, in various languages, with debug removed:
equation2_nl.f90
equation2_nl_f90.out
equation2_nl.c
equation2_nl_c.out
equation2_nl.java
equation2_nl_java.out
A polynomial nonlinear example that converges with no fctr control.
equation4_nl.c
equation4_nl_c.out
<- previous index next ->