# multiple eqn Newton Raphson method.
# developed bu Dagc.

from numpy import*
import math

print "Welcome to the NR-Multiple Root Calculator."
print "It works with only 3 simultaneous nonlinear eqn systems."
print "Please type your eqn's into the source code."
print "And call the 'multipleNR' function to work.."

def determinant(a):
    # with cofactor method, 3x3
        part1 = a[0,0]*(-1)**(0+0)*(a[1,1]*a[2,2]-a[1,2]*a[2,1])
        part2 = a[0,1]*(-1)**(0+1)*(a[1,0]*a[2,2]-a[1,2]*a[2,0])
        part3 = a[0,2]*(-1)**(0+2)*(a[1,0]*a[2,1]-a[1,1]*a[2,0])
        result = part1 + part2 + part3
        return result


#Since there are three functions:

def functionEqn1(x1,x2,x3):
    result = 3*x1-cos(x2*x3)-0.5
    return result

def functionEqn2(x1,x2,x3):
    result = x1**2 - 81*(x2+0.1)**2 + sin(x3) + 1.06
    return result

def functionEqn3(x1,x2,x3):
    result = e**(-x1*x2)+20*x3+(10*pi-3)/3
    return result


def multipleNR(x1, x2, x3, Es):

    Ea1 = 1.0    #assume
    Ea2 = 1.0
    Ea3 = 1.0

    counter = 0
    while (Ea1 > Es or Ea2 > Es or Ea3 > Es):
        # required derivatives
        delta = 10**(-5)
        d1x = (functionEqn1(x1+delta,x2,x3)-functionEqn1(x1,x2,x3))/delta
        d1y = (functionEqn1(x1,x2+delta,x3)-functionEqn1(x1,x2,x3))/delta
        d1z = (functionEqn1(x1,x2,x3+delta)-functionEqn1(x1,x2,x3))/delta

        d2x = (functionEqn2(x1+delta,x2,x3)-functionEqn2(x1,x2,x3))/delta
        d2y = (functionEqn2(x1,x2+delta,x3)-functionEqn2(x1,x2,x3))/delta
        d2z = (functionEqn2(x1,x2,x3+delta)-functionEqn2(x1,x2,x3))/delta

        d3x = (functionEqn3(x1+delta,x2,x3)-functionEqn3(x1,x2,x3))/delta
        d3y = (functionEqn3(x1,x2+delta,x3)-functionEqn3(x1,x2,x3))/delta
        d3z = (functionEqn3(x1,x2,x3+delta)-functionEqn3(x1,x2,x3))/delta

        # values of functions
        f1 = functionEqn1(x1,x2,x3)
        f2 = functionEqn2(x1,x2,x3)
        f3 = functionEqn3(x1,x2,x3)


        A = matrix([[d1x,d1y,d1z],[d2x,d2y,d2z],[d3x,d3y,d3z]])
        b = matrix([[-f1 + x1*d1x + x2*d1y + x3*d1z],[-f2 + x1*d2x + x2*d2y + x3*d2z],[-f3 + x1*d3x + x2*d3y + x3*d3z]])
        #print b

        denominator = determinant(A)
        #print denominator

        acolumn1 = matrix([[A[0,0]],
                           [A[1,0]],
                           [A[2,0]]])
        acolumn2 = matrix([[A[0,1]],
                           [A[1,1]],
                           [A[2,1]]])
        acolumn3 = matrix([[A[0,2]],
                           [A[1,2]],
                           [A[2,2]]])
##        print acolumn1
##        print acolumn2
##        print acolumn3

        x1matrice = matrix([[b[0,0],acolumn2[0,0],acolumn3[0,0]],
                            [b[1,0],acolumn2[1,0],acolumn3[1,0]],
                            [b[2,0],acolumn2[2,0],acolumn3[2,0]]])
        x2matrice = matrix([[acolumn1[0,0],b[0,0],acolumn3[0,0]],
                            [acolumn1[1,0],b[1,0],acolumn3[1,0]],
                            [acolumn1[2,0],b[2,0],acolumn3[2,0]]])
        x3matrice = matrix([[acolumn1[0,0],acolumn2[0,0],b[0,0]],
                            [acolumn1[1,0],acolumn2[1,0],b[1,0]],
                            [acolumn1[2,0],acolumn2[2,0],b[2,0]]])

        xn1 = determinant(x1matrice)/denominator
        xn2 = determinant(x2matrice)/denominator
        xn3 = determinant(x3matrice)/denominator


        Ea1 = float(abs((xn1-x1)/xn1))
        Ea2 = float(abs((xn2-x2)/x2))
        Ea3 = float(abs((xn3-x3)/xn3))
##        print Ea1
##        print Ea2
##        print Ea3

        x1 = xn1
        x2 = xn2
        x3 = xn3


        counter = counter + 1



    print "The roots are:"
    print "x1: ",x1
    print "x2: ", x2
    print "x3: ", x3
    print "respectively."
    print "Computation is completed in ",counter," iterations."