# Chapter 6 - Roots of non linear equations¶

## Example No. 6_01 Pg No. 126¶

In [1]:
from numpy import mat,shape
from math import sqrt

#Possible Initial guess values for roots

A = mat([[ 2],[-8] ,[2],[12]]) # Coefficients of x terms in the decreasing order of power
n = shape(A)#
x1 = -A[1]/A[0]#
print 'The largest possible root is  x1 =',x1
print 'No root can be larger than the value =',x1

x = sqrt((A[1]/A[0])**2 - 2*(A[2]/A[0])**2)

print '\nAll real roots lie in the interval (-%f,%f)\n'%(x,x)
print 'We can use these two points as initial guesses for the bracketing methods and one of them for open end methods'

The largest possible root is  x1 = [[4]]
No root can be larger than the value = [[4]]

All real roots lie in the interval (-3.741657,3.741657)

We can use these two points as initial guesses for the bracketing methods and one of them for open end methods


## Example No. 6_03 Pg No.¶

In [2]:
from numpy import mat,shape
from math import sqrt
#Evaluating Polynomial using Horner's rule

#Coefficients of x terms in the increasing order of power
A = mat([[6],[1],[-4],[1]])
x = 2
N=shape(A)
n,c = N[0],N[1]
p=[0,0,0]
p.append(A[n-1])
print 'p(4) =',p[n-1][0,0]
for i in range(1,n-1):
p[n-i]= p[n-i+1-1]*x + A[n-i-1]
print '\n p(%d)= %d\n'%(n-i,p[n-i])

print '\n f(%d) = p(1) = %d'%(x,p[0])

p(4) = 1

p(3)= -2

p(2)= 1

f(2) = p(1) = 0


## Example No. 6_04 Pg No. 132¶

In [3]:
from numpy import poly1d,polyval, sqrt

#Root of a Equation Using Bisection Method

#Coefficients in increasing order of power of x starting from 0
A = [-10 ,-4, 1]#
print 'First finding the interval that contains a root,this can be done by using Eq 6.10'
xmax = sqrt((A[1]/A[2])**2 - 2*(A[0]/A[2]))
print '\n Both the roots lie in the interval (-%d,%d) \n'%(xmax,xmax)
x = range(-6,7)
p= poly1d(A)# p = poly(A,'x'%('c'

fx = [p(xx) for xx in x]#
for i in range(0,12):
if fx[i]*fx[i] < 0:
break

print '\n The root lies in the interval (%d,%d)\n'%(x[i],x[i])

First finding the interval that contains a root,this can be done by using Eq 6.10

Both the roots lie in the interval (-6,6)

The root lies in the interval (5,5)



## Example No. 6_05 Pg No. 139¶

In [4]:
from numpy import poly1d,polyval

#False Position Method

#Coefficients of polynomial in increasing order of power of x
A = [-2 , -1,  1]
x1 = 1 #
x2 = 3 #
fx = poly1d(A)
for i in range(1,16):
print 'Iteration No. %d \n'%(i)
fx1 = fx(x1)
fx2 = fx(x2)
x0 = x1 - fx1*(x2-x1)/(fx2-fx1)
print 'x0 = %f \n'%(x0)#
fx0 = fx(x0)#
if fx1*fx0 < 0:
x2 = x0
else:
x1 = x0 #


Iteration No. 1

x0 = 1.000000

Iteration No. 2

x0 = 1.000000

Iteration No. 3

x0 = 1.000000

Iteration No. 4

x0 = 1.000000

Iteration No. 5

x0 = 1.000000

Iteration No. 6

x0 = 1.000000

Iteration No. 7

x0 = 1.000000

Iteration No. 8

x0 = 1.000000

Iteration No. 9

x0 = 1.000000

Iteration No. 10

x0 = 1.000000

Iteration No. 11

x0 = 1.000000

Iteration No. 12

x0 = 1.000000

Iteration No. 13

x0 = 1.000000

Iteration No. 14

x0 = 1.000000

Iteration No. 15

x0 = 1.000000



## Example No. 6_07 Pg No. 147¶

In [5]:
from numpy import poly1d,polyval,polyder
#Root of the Equation using Newton Raphson Method

#Coefficients of polynomial in increasing order of power of x
A = [ 2,  -3,  1]#
fx = poly1d(A)
dfx = polyder(fx)

x=[0]
f=[]
df=[]
for i in range(1,11):
f.append(fx(x[i-1]))
if f[i-1] != 0:
df.append(dfx(x[i-1]))
x.append(x[i-1] - f[i-1]/df[i-1])
print 'x%d = %f\n'%(i+1,x[i])#
else:
print 'Since f(%f) = 0, the root closer to the point x = 0 is %f \n'%(x[i-1],x[i-1] )
break


x2 = 1.000000

Since f(1.000000) = 0, the root closer to the point x = 0 is 1.000000



## Example No. 6_08 Pg No. 151¶

In [6]:
from numpy import poly1d,polyval,polyder

#Root of the Equation using Newton Raphson Method

#Coefficients of polynomial in increasing order of power of x
A = [ 6,  1 , -4 , 1 ]
fx = poly1d(A)
dfx = polyder(fx)
f=[];df=[]
x = [5.0] #
for i in range(1,7):
f.append(fx(x[i-1]))
if f[i-1] != 0:
df.append(dfx(x[i-1]))
x.append(x[i-1] - f[i-1]/df[i-1])
print 'x%d = %f\n'%(i+1,x[i])

print 'From the results we can see that number of correct digits approximately doubles with each iteration'

x2 = 3.342105

x3 = 2.248631

x4 = 1.535268

x5 = 1.079139

x6 = 0.797330

x7 = 0.632716

From the results we can see that number of correct digits approximately doubles with each iteration


## Example No. 6_09 Pg No. 153¶

In [7]:
from __future__ import division
from numpy import poly1d
#Root of the equation using SECANT Method

#Coefficients of polynomial in increasing order of power of x
A = [ -10,  -4,  1]
x1 = 4 #
x2 = 2 #
fx = poly1d(A)
for i in range(1,7):
print '\n For Iteration No. %d\n'%(i)
fx1 = fx(x1)
fx2 = fx(x2)
x3 = x2 - fx2*(x2-x1)/(fx2-fx1) #
print '\n x1 = %f\n x2 = %f \n fx1 = %f \n fx2 = %f \n x3 = %f \n'%(x1,x2,fx1,fx2,x3) #
x1 = x2#
x2 = x3#

print 'This can be still continued further for accuracy'

 For Iteration No. 1

x1 = 4.000000
x2 = 2.000000
fx1 = -175.000000
fx2 = -47.000000
x3 = 1.265625

For Iteration No. 2

x1 = 2.000000
x2 = 1.265625
fx1 = -47.000000
fx2 = -20.080566
x3 = 0.717818

For Iteration No. 3

x1 = 1.265625
x2 = 0.717818
fx1 = -20.080566
fx2 = -7.023891
x3 = 0.423122

For Iteration No. 4

x1 = 0.717818
x2 = 0.423122
fx1 = -7.023891
fx2 = -2.482815
x3 = 0.261999

For Iteration No. 5

x1 = 0.423122
x2 = 0.261999
fx1 = -2.482815
fx2 = -0.734430
x3 = 0.194317

For Iteration No. 6

x1 = 0.261999
x2 = 0.194317
fx1 = -0.734430
fx2 = -0.154860
x3 = 0.176233

This can be still continued further for accuracy


## Example No. 6_11 Pg No. 161¶

In [8]:
from numpy import poly1d
from __future__ import division
#Fixed point method

#Coefficients of polynomial in increasing order of power of x
A = [ -2,  1,  1 ]#
B = [ 2, 0, -1 ]#
gx = poly1d(B)
x = [0] ##initial guess x0 = 0
for i in range(2,11):
x.append (gx(x[i-2]))
print '\n x%d = %f\n'%(i-1,x[i-1])
if (x[i-1]-x[(i-2)]) == 0:
print '\n%f is root of the equation,since x%d - x%d = 0 \n'%(x[i-1],i-1,i-2)
break

#Changing initial guess x0 = -1
x[0] = -1 #
for i in range(2,11):
x[i-1]= gx(x[i-2])
print '\nx%d = %f\n'%(i-1,x[i-1])
if (x[i-1]-x[i-2]) == 0:
print '\n %f is root of the equation,since x%d - x%d = 0'%(x[i-1],i-1,i-2)
break

 x1 = -1.000000

x2 = 1.000000

x3 = 1.000000

1.000000 is root of the equation,since x3 - x2 = 0

x1 = 1.000000

x2 = 1.000000

1.000000 is root of the equation,since x2 - x1 = 0


## Example No. 6_12 Pg No. 162¶

In [9]:
#Fixed point method

A = [ -5,  0,  1 ]#
def g(x):
x = 5.0/x
return x
x = [1] #
print '\n x0 = %f \n'%(x[0])
for i in range(2,6):
x.append(g(i))
print ' x%d = %f \n'%(i-1,x[i-1])

#Defining g(x) in different way
def g(x):
x = x**2 + x - 5
return x
x=[0]
print '\n x0 = %f \n'%(x[0])
for i in range(2,6):
x.append(g(i))
print ' x%d = %f \n'%(i-1,x[i-1])

#Third form of g(x)
def g(x):
x = (x + 5/x)/2
return x
x=[1]
print '\n x0 = %f \n'%(x[0])
for i in range(2,8):
x.append(g(i))
print ' x%d = %f \n'%(i-1,x[i-1])

 x0 = 1.000000

x1 = 2.500000

x2 = 1.666667

x3 = 1.250000

x4 = 1.000000

x0 = 0.000000

x1 = 1.000000

x2 = 7.000000

x3 = 15.000000

x4 = 25.000000

x0 = 1.000000

x1 = 2.250000

x2 = 2.333333

x3 = 2.625000

x4 = 3.000000

x5 = 3.416667

x6 = 3.857143



## Example No. 6_13 Pg No. 169¶

In [10]:
#Solving System of non-linear equations using FIXED POINT METHOD

print ' x**2 - y**2 = 3 \n x**2 + x*y \n'
def f(x,y):
x = y + 3/(x+y)
return x
def g(x):
y = (6-x**2)/x
return y
x=[1]
y=[1]
print '\n x0 = %f \n y0 = %f \n'%(x[0],y[0])
for i in range(2,5):
x.append(f((i-1),(i-1)))
y.append(g((i-1)))
print '\n x%d = %f \n y%d = %f \n'%(i-1,x[i-1],i-1,y[i-1])

 x**2 - y**2 = 3
x**2 + x*y

x0 = 1.000000
y0 = 1.000000

x1 = 2.500000
y1 = 5.000000

x2 = 2.750000
y2 = 1.000000

x3 = 3.500000
y3 = -1.000000



## Example No. 6_14 Pg No. 172¶

In [11]:
#Solving System of Non-linear equations using Newton Raphson Method

print 'x**2 + x*y = 6 \n x**2 - y**2 = 3 \n'#
def F(x,y):
f = x**2 + x*y - 6
return f
def G(x,y):
g = x**2 - y**2 -3
return g
def dFx(x,y):
f1 = 2*x + y
return f1
def dFy(x,y):
f2 = y
return f2
def dGx(x,y):
g1 = 2*x
return g1
def dGy(x,y):
g2 = -2*y
return g2
x=[1]
y=[1]

for i in range(2,4):
Fval = F(i,i)
Gval = G(i,i)
f1 = dFx(i-1,i-1)
f2 = dFy(i-1,i-1)
g1 = dGx(i-1,i-1)
g2 = dGy(i-1,i-1)
D =  f1*g2 - f2*g1

x.append(x[i-2] - (Fval*g2 - Gval*f2)/D )
y.append(y[i-2] - (Gval*f1 - Fval*g1)/D )
print '\n x%d = %f \n y%d = %f \n'%(i-1,x[i-1],i-1,y[i-1])

x**2 + x*y = 6
x**2 - y**2 = 3

x1 = 0.875000
y1 = -0.625000

x2 = -0.437500
y2 = -2.687500



## Example No. 6_15 Pg No. 176¶

In [12]:
from __future__ import division
from numpy import poly1d, arange
#Synthetic Division

a = [-9 ,15 ,-7, 1]
b=[0,0,0,0]
for i in arange(3,0,-1):
b[i]= a[i] + b[i]*3
print 'b%d = %f\n'%(i,b[i-1])
print 'Thus the polynomial is'
print poly1d(b)

b3 = 0.000000

b2 = 0.000000

b1 = 0.000000

Thus the polynomial is
2
15 x - 7 x + 1


## Example No. 6_16 Pg No. 187¶

In [13]:
#Quadratic factor of a polynomial using Bairstow's Method

a = [ 10, 1 ,0 ,1]#
n = len(a)#
u = 1.8 #
v = -1 #
b=[];c=[]
for nn in range(n):
b.append(0)
c.append(0)
b[n-1] = a[n-1]
b[n-2] = a[n-2] + u*b[n-1]
c[n-1] = 0
c[n-2] = b[n-1]

for i in range(n-2,0,-1):
b[i-1]= a[i-1]+ u*b[i] + v*b[i+1]
c[i-1]= b[i] + u*c[i] + v*c[i+1]

for i in range(n,0,-1):
print 'b%d = %f \n'%(i-1,b[i-1])

for i in range(n,0,-1):
print 'c%d = %f \n'%(i-1,b[i-1])

D = c[1]*c[1] - c[0]*c[2]
du = -1*(b[1]*c[1] - c[0]*c[2])/D
dv = -1*(b[0]*c[1] - b[1]*c[0])/D
u = u + du #
v = v + du #
print '\n D = %f \n du = %f \n dv = %f \n u = %f\n v = %f \n'%(D,du,dv,u,v)


b3 = 1.000000

b2 = 1.800000

b1 = 3.240000

b0 = 14.032000

c3 = 1.000000

c2 = 1.800000

c1 = 3.240000

c0 = 14.032000

D = 4.240000
du = -0.694340
dv = -5.250566
u = 1.105660
v = -1.694340



## Example No. 6_17 Pg No. 197¶

In [14]:
from math import sqrt
#Solving Leonard's equation using MULLER'S Method

def f(x):
y = x**3 + 2*x**2 + 10*x - 20
return y
x1 = 0 #
x2 = 1 #
x3 = 2 #
for i in range(1,11):
f1 = f(x1)
f2 = f(x2)
f3 = f(x3)
h1 = x1-x3 #
h2 = x2-x3 #
d1 = f1 - f3 #
d2 = f2 - f3 #
D = h1*h2*(h1-h2)#
a0 = f3 #
a1 = (d2*h1**2 - d1*h2**2)/D #
a2 = (d1*h2 - d2*h1)/D #
if abs(-2*a0/( a1 + sqrt( a1**2 - 4*a0*a2 ) )) < abs( -2*a0/( a1 - sqrt( a1**2 - 4*a0*a2 ) )):
h4 = -2*a0/(a1 + sqrt(a1**2 - 4*a0*a2))#
else:
h4 = -2*a0/(a1 - sqrt(a1**2 - 4*a0*a2))

x4 = x3 + h4 #
print '\n x1 = %f\n x2 = %f\n x3 = %f\n f1 = %f\n f2 = %f\n f3 = %f\n h1 = %f\n h2 = %f\n d1 = %f\n d2 = %f\n a0 = %f\n a1 = %f\n a2 = %f\n h4 = %f\n x4 = %f\n '%(x1,x2,x3,f1,f2,f3,h1,h2,d1,d2,a0,a1,a2,h4,x4) #
relerr = abs((x4-x3)/x4)#
if relerr <= 0.00001:
print 'root of the polynomial is x4 = %f'%(x4)
break

x1 = x2 #
x2 = x3 #
x3 = x4 #

 x1 = 0.000000
x2 = 1.000000
x3 = 2.000000
f1 = -20.000000
f2 = -7.000000
f3 = 16.000000
h1 = -2.000000
h2 = -1.000000
d1 = -36.000000
d2 = -23.000000
a0 = 16.000000
a1 = 28.000000
a2 = 5.000000
h4 = -0.645934
x4 = 1.354066

x1 = 1.000000
x2 = 2.000000
x3 = 1.354066
f1 = -7.000000
f2 = 16.000000
f3 = -0.309679
h1 = -0.354066
h2 = 0.645934
d1 = -6.690321
d2 = 16.309679
a0 = -0.309679
a1 = 21.145451
a2 = 6.354066
h4 = 0.014581
x4 = 1.368647

x1 = 2.000000
x2 = 1.354066
x3 = 1.368647
f1 = 16.000000
f2 = -0.309679
f3 = -0.003394
h1 = 0.631353
h2 = -0.014581
d1 = 16.003394
d2 = -0.306286
a0 = -0.003394
a1 = 21.103381
a2 = 6.722713
h4 = 0.000161
x4 = 1.368808

x1 = 1.354066
x2 = 1.368647
x3 = 1.368808
f1 = -0.309679
f2 = -0.003394
f3 = -0.000001
h1 = -0.014742
h2 = -0.000161
d1 = -0.309678
d2 = -0.003392
a0 = -0.000001
a1 = 21.096136
a2 = 6.091521
h4 = 0.000000
x4 = 1.368808

root of the polynomial is x4 = 1.368808