<- previous    index    next ->

Lecture 16, Finding Roots and Nonlinear Equations 



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.
find_roots.m source code
find_roots_m.out output

Maple root finding for polynomials, not great:
find_roots.mws source code
find_roots_mws.out output
real and complex work some times.

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.

"cpoly" is a very old routine to compute roots of a polynomial
with complex coefficients c_0 ... c_n.

The Fortran code, including the test program (driver) for  cpoly  is
419.for in Fortran IV, compiles with g95
419_f.out is the output of many test cases

Another version, with an automatic conversion of Fortran to C.
The program  f2c  converts .f to .c  it may not be very pretty.
It also needs the library that comes with  f2c.
419.f driver and cpoly
419.c driver and cpoly


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.
    Do not guess zero for y, x/y does not work, use y = 1 as guess.


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.

Given: f(x) = 0
guess x
x_next = x - f(x)/f'(x)   f'(x) is the derivative of f(x)
                          avoid region where slope is zero
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:
given  f(x) = x + 2 x^2 + 3/x = 6 
thus   f'(x) = 1 + 4 x - 6/x^2
then  guess a value of x,  Newton iteration is
next  x = x - f(x)/f'(x)

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

in matrix form, the equations are: A * X = Y
| 1  2  3 |   |  x1  | = | 10.000 |
| 2  1  1 | * | x2^2 | = |  6.333 |
| 3  1  4 |   | 1/x3 | = |  8.333 |
     A      *    X     =     Y

from equations that have these derivatives
 f1(x1, x2, x3) = x1 + 2 x2^2 + 3/x3 = 10
 f1'(x1) = 1
 f1'(x2) = 4 x2
 f1'(x3) = -3/x3^2

 f2(x1, x2, x3) = 2 x1 + x2^2 + 1/x3 = 6.333
 f2'(x1) = 2
 f2'(x2) = 2 x2
 f2'(x3) = -1/x3^2

 f3(x1, x2, x3) = 3 x1 + x2^2 + 4/x3 = 8.333
 f3'(x1) = 3
 f3'(x2) = 2 x2
 f3'(x3) = -4/x3^2

create the Jacobian 
     |   1    |   | 1   4*x2  -3/x3^2 |
 A * | 2*x2   | = | 2   2*x2  -1/x3^2 | = Ja
     |-1/x3^2 |   | 3   2*x2  -4/x3^2 |
 A *     D      =        Ja
     deriv of X   wrt x1,x2,x3

The Newton iteration becomes
 next X = X - (A*X-Y)/Ja,  /Ja is times transpose inverse of Ja
 Ja is, in general, dependent on all variables

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

A more general routine for solving a system of nonlinear
simultaneous equations A x = y where there are n unknowns and
any unknown may be first, second, third power of negative one
or negative two power:
A is n by n+nonlinear terms.
y is n right hand sides of equations
x is n+nonlinear, initially the first guess at solution 
of the linear terms with space for non linear terms.

The caller sets up var1 with indices of first power.
The caller sets up var2 with indices of second power.
The caller sets up var3 with indices of third power.
The caller sets up vari1 with indices of negative first power.
The caller sets up vari2 with indices of negative second power.
Must have each term, then non linear, note: -1 means no term

a*x[0] + b*x[0]*x[1]*x[2] + c*x[1]*x[1]*x1[1]/(x[0]*x[2]) = d
n=3 variables, nonlinear=2 in this example

              x[0] x[1] x[2]   b    c
int var1[]  = { 0,   1,   2,   0,   1 }; // zero based subscript          
int var2[]  = {-1,  -1,  -1,   1,   1 }; // e.g. last is x3*x4            
int var3[]  = {-1,  -1,  -1,   2,   1 }; // possible cube or 3 term       
int vari1[] = {-1,  -1,  -1,  -1,   0 }; // zero based subscript          
int vari2[] = {-1,  -1,  -1,  -1,   2 }; // e.g. last is x3*x4            

simeq_newton5.java
test_simeq_newton5.java
test_simeq_newton5_java.out
inverse.java

simeq_newton5.py3
test_simeq_newton5.py3
test_simeq_newton5_py3.out



    <- previous    index    next ->

Other links

Go to top