#example 2.1
#bisection method
#page 24
import math
def f(x):
return math.pow(x,3)-x-1
x1=1
x2=2 #f(1) is negative and f(2) is positive
d=0.0001 #for accuracy of root
c=1
print "Succesive approximations \t x1\t \tx2\t \tm\t \tf(m)\n"
while abs(x1-x2)>d:
m=(x1+x2)/2.0
print " \t%f\t%f\t%f\t%f\n" %(x1,x2,m,f(m))
if f(m)*f(x1)>0.0:
x1=m
else:
x2=m
c=c+1 # to count number of iterations
print "the solution of equation after %i iteration is %g" %(c,m)
#example 2.2
#bisection method
#page 25
import math
def f(x):
return math.pow(x,3)-2*x-5
x1=2
x2=3 #f(2) is negative and f(3) is positive
d=0.0001 #for accuracy of root
c=1
print "Succesive approximations \t x1\t \tx2\t \tm\t \tf(m)\n"
while abs(x1-x2)>d:
m=(x1+x2)/2.0
print " \t%f\t%f\t%f\t%f\n" %(x1,x2,m,f(m))
if f(m)*f(x1)>0:
x1=m
else:
x2=m
c=c+1;# to count number of iterations
print "the solution of equation after %i iteration is %0.4g" %(c,m)
#example 2.3
#bisection method
#page 26
import math
def f(x):
return math.pow(x,3)+math.pow(x,2)+x+7
x1=-3
x2=-2 #f(-3) is negative and f(-2) is positive
d=0.0001 #for accuracy of root
c=1
print "Succesive approximations \t x1\t \tx2\t \tm\t \tf(m)\n"
while abs(x1-x2)>d:
m=(x1+x2)/2.0
print " \t%f\t%f\t%f\t%f\n" %(x1,x2,m,f(m))
if f(m)*f(x1)>0:
x1=m
else:
x2=m
c=c+1 # to count number of iterations
print "the solution of equation after %i iteration is %0.4g" %(c,m)
#example 2.4
#bisection method
#page 26
import math
def f(x):
return x*math.exp(x)-1
x1=0
x2=1 #f(0) is negative and f(1) is positive
d=0.0005 #maximun tolerance value
c=1
print "Succesive approximations \t x1\t \tx2\t \tm\t \ttol\t \tf(m)\n"
while abs((x2-x1)/x2)>d:
m=(x1+x2)/2.0 #tolerance value for each iteration
tol=((x2-x1)/x2)*100
print " \t%f\t%f\t%f\t%f\t%f\n" %(x1,x2,m,tol,f(m))
if f(m)*f(x1)>0:
x1=m
else:
x2=m
c=c+1 # to count number of iterations
print "the solution of equation after %i iteration is %0.4g" %(c,m)
#example 2.5
#bisection method
#page 27
import math
def f(x):
return 4*math.exp(-x)*math.sin(x)-1
x1=0
x2=0.5 #f(0) is negative and f(1) is positive
d=0.0001 #for accuracy of root
c=1
print "Succesive approximations \t x1\t \tx2\t \tm\t \t \tf(m)\n"
while abs(x2-x1)>d:
m=(x1+x2)/2.0
print " \t%f\t%f\t%f\t%f\n" %(x1,x2,m,f(m))
if f(m)*f(x1)>0:
x1=m
else:
x2=m
c=c+1 # to count number of iterations
print "the solution of equation after %i iteration is %0.3g" %(c,m)
#example 2.6
#false position method
#page 28
import math
def f(x):
return x**3-2*x-5
a=2.0
b=3.0 #f(2) is negative and f(3)is positive
d=0.00001
print "succesive iterations \ta\t b\t f(a)\t f(b)\t\ x1\n"
for i in range(1,25):
x1=b*f(a)/(f(a)-f(b))+a*f(b)/(f(b)-f(a))
if(f(a)*f(x1))>0:
b=x1
else:
a=x1
if abs(f(x1))<d:
break
print " \t%f %f %f %f %f\n" %(a,b,f(a),f(b),x1)
print "the root of the equation is %f" %(x1)
#example 2.7
#false position method
#page 29
def f(x):
return x**2.2-69
a=5.0
b=6.0 #f(5) is negative and f(6)is positive
d=0.00001
print "succesive iterations \ta\t b\t f(a)\t f(b)\t\ x1\n"
for i in range(1,25):
x1=b*f(a)/(f(a)-f(b))+a*f(b)/(f(b)-f(a));
if(f(a)*f(x1))>0:
b=x1
else:
a=x1
if abs(f(x1))<d:
break
print " \t%f %f %f %f %f\n" %(a,b,f(a),f(b),x1)
print "the root of the equation is %f" %(x1)
#example 2.8
#false position method
#page 29
import math
def f(x):
return 2*x-log10(x)-7
a=3.0
b=4.0 #f(3) is negative and f(4)is positive
d=0.00001
print "succesive iterations \ta \t b\t f(a)\t f(b)\t x1\n"
for i in range(1,25):
x1=b*f(a)/(f(a)-f(b))+a*f(b)/(f(b)-f(a))
if(f(a)*f(x1))>0:
b=x1
else:
a=x1
if abs(f(x1))<d:
break
print " \t%f %f %f %f %f\n" %(a,b,f(a),f(b),x1)
print "the root of the equation is %0.4g" %(x1)
#example 2.9
#false position method
#page 30
import math
def f(x):
return 4*math.exp(-x)*math.sin(x)-1
a=0.0
b=0.5 #f(0) is negative and f(0.5)is positive
d=0.00001
print "succesive iterations \ta\t b\t f(a)\t f(b)\t\ x1\n"
for i in range(1,25):
x1=b*f(a)/(f(a)-f(b))+a*f(b)/(f(b)-f(a))
if(f(a)*f(x1))>0:
b=x1
else:
a=x1
if abs(f(x1))<d:
break
print " \t%f %f %f %f %f\n" %(a,b,f(a),f(b),x1)
print "the root of the equation is %f" %(x1)
#example 2.10
#iteration method
#page 33
import math
def f(x):
return 1/(math.sqrt(x+1))
x1=0.75
x2=0.0
n=1
d=0.0001 #accuracy opto 10^-4
c=0 #to count no of iterations
print "successive iterations \t\x01\tf(x1)\n"
while abs(x1-x2)>d:
print " \t%f %f\n" %(x1,f(x1))
x2=x1
x1=f(x1)
c=c+1
print " the root of the eqaution after %i iteration is %0.4g" %(c,x1)
#example 2.11
#iteration method
#page34
import math
def f(x):
return cos(x)/2.0+3.0/2.0
x1=1.5 # as roots lies between 3/2 and pi/2
x2=0
d=0.0001 # accuracy opto 10^-4
c=0 # to count no of iterations
print "successive iterations \t\x01\tf(x1)\n"
while abs(x2-x1)>d:
print " \t%f %f\n" %(x1,f(x1))
x2=x1
x1=f(x1)
c=c+1
print " the root of the eqaution after %i iteration is %0.4g" %(c,x1)
#example 2.12
#iteration method
#page 35
import math
def f(x):
return math.exp(-x)
x1=1.5 # as roots lies between 0 and 1
x2=0
d=0.0001 # accuracy opto 10^-4
c=0 # to count no of iterations
print "successive iterations \t x1 \t f(x1)\n"
while abs(x2-x1)>d:
print " \t%f %f\n" %(x1,f(x1))
x2=x1
x1=f(x1)
c=c+1
print " the root of the eqaution after %i iteration is %0.4g" %(c,x1)
#example 2.13
#iteration method
#page 35
import math
def f(x):
return 1+math.sin(x)/10
x1=1.0 # as roots lies between 1 and pi evident from graph
x2=0
d=0.0001 # accuracy opto 10^-4
c=0 # to count no of iterations
print "successive iterations \t x1 \t f(x1)\n"
while abs(x2-x1)>d:
print " \t%f %f\n" %(x1,f(x1))
x2=x1
x1=f(x1)
c=c+1
print " the root of the eqaution after %i iteration is %0.4g" %(c,x1)
#example 2.14
#aitken's process
#page 36
import math
def f(x):
return 1.5+math.cos(x)/2.0
x0=1.5
y=0
e=0.0001
c=0
print "successive iterations \t x0 \t x1 \t x2 \t x3 \t y\n"
for i in range(1,10):
x1=f(x0)
x2=f(x1)
x3=f(x2)
y=x3-((x3-x2)**2)/(x3-2*x2+x1)
d=y-x0
x0=y
if abs(f(x0))<e:
break
c=c+1
print " \t%f %f %f %f %f\n" %(x0,x1,x2,x3,y)
print "the root of the equation after %i iteration is %f" %(c,y)
#example 2.15
#newton-raphson method
#page 39
import math
def f(x):
return x**3-2*x-5
def f1(x):
return 3*x**2-2 # first derivative of the function
x0=2.0 # initial value
d=0.0001
c=0
n=1
print "successive iterations \t x0 \t f(x0) \t f1(x0)\n"
while n==1:
x2=x0
x1=x0-(f(x0)/f1(x0))
x0=x1
print " \t%f \t%f \t%f \n" %(x2,f(x1),f1(x1))
c=c+1
if abs(f(x0))<d:
break
print "the root of %i iteration is:%f" %(c,x0)
#example 2.16
#newton-raphson method
#page 40
import math
def f(x):
return x*math.sin(x)+math.cos(x)
def f1(x):
return x*math.cos(x) #first derivation of the function
x0=math.pi # initial value
d=0.0001
c=0
n=1
print "successive iterations \tx0\t f(x0)\t f1(x0)\n"
while n==1:
x2=x0
x1=x0-(f(x0)/f1(x0))
x0=x1
print " \t%f \t%f \t%f\n" %(x2,f(x1),f1(x1))
c=c+1
if abs(f(x0))<d:
break
print "the root of %i iteration is:%f" %(c,x0)
#example 2.17
#newton-raphson method
#page 40
import math
def f(x):
return x*math.exp(x)-1
def f1(x):
return math.exp(x)+x*math.exp(x) #first derivative of the function
x0=0 # initial value
d=0.0001
c=0
n=1
print "successive iterations \tx0\t f(x0)\t f1(x0)\n"
while n==1:
x2=x0
x1=x0-(f(x0)/f1(x0))
x0=x1
print " \t%f \t%f \t%f\n" %(x2,f(x1),f1(x1))
c=c+1
if abs(f(x0))<d:
break
print "the root of %i iteration is:%f" %(c,x0)
#example 2.18
#newton-raphson method
#page 41
import math
def f(x):
return math.sin(x)-x/2.0
def f1(x):
return math.cos(x)-0.5
x0=math.pi/2.0 # initial value
d=0.0001
c=0
n=1
print "successive iterations \t x0 \t f(x0)\t f1(x0)\n"
while n==1:
x2=x0
x1=x0-(f(x0)/f1(x0))
x0=x1
print " \t%f\t%f\t%f\n" %(x2,f(x1),f1(x1))
c=c+1
if abs(f(x0))<d:
break
print "the root of %i iteration is: %0.4g" %(c,x0)
#example 2.19
#newton-raphson method
#page 41
import math
def f(x):
return 4*math.exp(-x)*math.sin(x)-1
def f1(x):
return math.cos(x)*4*math.exp(-x)-4*math.exp(-x)*math.sin(x)
x0=0.2 # initial value
d=0.0001
c=0
n=1
print "successive iterations \t x0 \t f(x0)\t f1(x0)\n"
while n==1:
x2=x0
x1=x0-(f(x0)/f1(x0))
x0=x1
print " \t%f \t%f \t%f\n" %(x2,f(x1),f1(x1))
c=c+1
if abs(f(x0))<d:
break
print "the root of %i iteration is: %0.3g" %(c,x0)
#example 2.20
#generalized newton-raphson method
#page 42
import math
def f(x):
return x**3-x**2-x+1
def f1(x):
return 3*x**2-2*x-1
def f2(x):
return 6*x-2
x0=0.8 # initial value to finf double root
n=1
print "successive iterations \t x0 \t x1\t x2\n"
while n==1:
x1=x0-(f(x0)/f1(x0));
x2=x0-(f1(x0)/f2(x0));
if abs(x1-x2)<0.000000001:
x0=(x1+x2)/2.0
break
else:
x0=(x1+x2)/2;
print " %f\t %f\t %f\n" %(x0,x1,x2)
print "\n \nthe double root is at: %f" %(x0)
#ramanujan's method
#example 2.21
#page 45
def f(x):
return 1-(13.0/12.0)*x-(3.0/8.0)*x**2+(1.0/24.0)*x**3
a1=13.0/12.0
a2=-3.0/8.0
a3=1.0/24.0
b1=1
b2=a1
b3=a1*b2+a2*b1
b4=a1*b3+a2*b2+a3*b1
b5=a1*b4+a2*b3+a3*b2
b6=a1*b5+a2*b4+a3*b3
b7=a1*b6+a2*b5+a3*b4
b8=a1*b7+a2*b6+a3*b5
b9=a1*b8+a2*b7+a3*b6
print "\n\n%f" %(b1/b2)
print "\n%f" %(b2/b3)
print "\n%f" %(b3/b4)
print "\n%f" %(b4/b5)
print "\n%f" %(b5/b6)
print "\n%f" %(b6/b7)
print "\n%f" %(b7/b8)
print "\n%f" %(b8/b9)
print "\n it appears as if the roots are converging at 2"
#ramanujan's method
#example 2.22
#page 46
def f(x):
return x+x**2+(x**3)/2.0+(x**4)/6.0+(x**5)/24.0
a1=1.0
a2=1.0
a3=1.0/2.0
a4=1.0/6.0
a5=1.0/24.0
b1=1
b2=a2
b3=a1*b2+a2*b1
b4=a1*b3+a2*b2+a3*b1
b5=a1*b4+a2*b3+a3*b2
b6=a1*b5+a2*b4+a3*b3
print "\n%f" %(b1/b2)
print "\n%f" %(b2/b3)
print "\n%f" %(b3/b4)
print "\n%f" %(b4/b5)
print "\n%f" %(b5/b6)
print "\n it appears as if the roots are converging at around %f" %(b5/b6)
#ramanujan's method
#example 2.23
#page 47
from __future__ import division
def f(x):
return 1-2*((3*x/2.0+(x**2)/4.0-(x**4)/48.0+(x**6)/1440.0)-(x**8)*2/80640.0)
a1=3/2
a2=1/4
a3=0
a4=1/48
a5=0
a6=1/1440
a7=0
a8=-1/80640
b1=1
b2=a1
b3=a1*b2+a2*b1
b4=a1*b3+a2*b2+a3*b1
b5=a1*b4+a2*b3+a3*b2
b6=a1*b5+a2*b4+a3*b3
b7=a1*b6+a2*b5+a3*b4
b8=a1*b7+a2*b6+a3*b5
b9=a1*b8+a2*b7+a3*b6
print "\n%f" %(b1/b2)
print "\n%f" %(b2/b3)
print "\n%f" %(b3/b4)
print "\n%f" %(b4/b5)
print "\n%f" %(b5/b6)
print "\n%f" %(b6/b7)
print "\n%f" %(b7/b8)
print "\n it appears as if the roots are converging at around %f" %(b7/b8)
#ramanujan's method
#example 2.24
#page 47
import math
def f(x):
return 1-(x-x**2.0/math.factorial(2.0)**2.0+x**3.0/math.factorial(3.0)**2.0-x**4.0/math.factorial(4.0)**2.0)
a1=1
a2=-1/math.factorial(2.0)**2.0
a3=1/math.factorial(3.0)**2.0
a4=-1/math.factorial(4.0)**2.0
a5=-1/math.factorial(5.0)**2.0
a6=1/math.factorial(6.0)**2.0
b1=1
b2=a1
b3=a1*b2+a2*b1
b4=a1*b3+a2*b2+a3*b1
b5=a1*b4+a2*b3+a3*b2
print "\n\n%f" %(b1/b2)
print "\n\n%f" %(b2/b3)
print "\n%f" %(b3/b4)
print "\n%f" %(b4/b5)
print "\n it appears as if the roots are converging at around %f" %(b4/b5)
#example 2.25
#secant method
#page 49
import math
from __future__ import division
def f(x):
return x**3-2*x-5
x1=2
x2=3 # initial values
n=1
c=0
print "successive iterations \t x1 \t x2 \t x3 \t f(x3)\n"
while n==1:
x3=(x1*f(x2)-x2*f(x1))/(f(x2)-f(x1))
print " \t%f \t%f \t%f \t%f\n" %(x1,x2,x3,f(x3))
if f(x3)*f(x1)>0:
x2=x3;
else:
x1=x3
if abs(f(x3))<0.000001:
break
c=c+1
print "the root of the equation after %i iteration is: %f" %(c,x3)
#example 2.26
#secant method
#page 50
import math
from __future__ import division
def f(x):
return x*math.exp(x)-1
x1=0
x2=1 # initial values
n=1
c=0
print "successive iterations \t x1 \t x2 \t x3 \t f(x3)\n"
while n==1:
x3=(x1*f(x2)-x2*f(x1))/(f(x2)-f(x1))
print " \t%f \t%f \t%f \t%f\n" %(x1,x2,x3,f(x3))
if f(x3)*f(x1)>0:
x2=x3
else:
x1=x3
if abs(f(x3))<0.0001:
break
c=c+1
print "the root of the equation after %i iteration is: %0.4g" %(c,x3)
# example 2.27
#mulller's method
#page 52
from __future__ import division
import math
def f(x):
return x**3-x-1
x0=0
x1=1
x2=2 # initial values
n=1
c=0
print "successive iterations \t x0 \t x1 \t x2 \t f(x0)\t f(x1)\t f(x2)\n"
while n==1:
c=c+1
y0=f(x0)
y1=f(x1)
y2=f(x2)
h2=x2-x1
h1=x1-x0
d2=f(x2)-f(x1)
d1=f(x1)-f(x0)
print " \t%f\t %f\t %f\t %f\t %f\t %f\n" %(x0,x1,x2,f(x0),f(x1),f(x2))
A=(d2/h2-d1/h1)/(h1+h2)
B=d2/h2+A*h2
S=math.sqrt(B**2-4*A*f(x2))
x3=x2-(2*f(x2))/(B+S)
E=abs((x3-x2)/x2)*100
if E<0.003:
break
else:
if c==1:
x2=x3
if c==2:
x1=x2
x2=x3
if c==3:
x0=x1
x1=x2
x2=x3
if c==3:
c=0
print "the required root is : %0.4f" %(x3)
#graeffe's method
#example 2.28
#page 55
import math
from __future__ import division
def f(x):
return x**3-6*(x**2)+11*x-6
print "the equation is:\n"
A=[1, 14, 49, 36] #coefficients of the above equation
print "%0.4g\n" %(math.sqrt(A[3]/A[2]))
print "%0.4g\n" %(math.sqrt(A[2]/A[1]))
print "%0.4g\n" %(math.sqrt(A[1]/A[0]))
print "the equation is:\n"
B=[1, 98, 1393, 1296]
print "%0.4g\n" %((B[3]/B[2])**(1/4))
print "%0.4g\n" %((B[2]/B[1])**(1/4))
print "%0.4g\n" %((B[1]/B[0])**(1/4))
print "It is apparent from the outputs that the roots converge at 1 2 3"
#quadratic factor by lin's--bairsttow method
#example 2.29
#page 57
from numpy import matrix
from __future__ import division
def f(x):
return x**3-x-1
a=[-1, -1, 0, 1]
r1=1
s1=1
b4=a[3]
def f3(r):
return a[2]-r*a[3]
def f2(r,s):
return a[1]-r*a[2]+r**2*a[3]-s*a[3]
def f1(r,s):
return a[0]-s*a[2]+s*r*a[3]
A=matrix([[1,1],[2,-1]])
C=matrix([[0],[1]])
X=A.I*C
X1=[[ 0.33333333],[-0.33333333]]
dr=X1[0][0]
ds=X1[1][0]
r2=r1+dr
s2=s1+ds
#second pproximation
r1=r2
s1=s2
b11=f1(r2,s2)
b22=f2(r2,s2)
h=0.001
dr_b1=(f1(r1+h,s1)-f1(r1,s1))/h
ds_b1=(f1(r1,s1+h)-f1(r1,s1))/h
dr_b2=(f2(r1+h,s1)-f2(r1,s1))/h
ds_b2=(f2(r1,s1+h)-f2(r1,s1))/h
A=matrix([[dr_b1,ds_b1],[dr_b2,ds_b2]])
C=matrix([[-f1(r1,s1)],[-f2(r1,s2)]])
X=A.I*C
r2=r1+X[0][0]
s2=s1+X[1][0]
print "roots correct to 3 decimal places are : %0.3f %0.3f" %(r2,s2)
#method of iteration
#example 2.31
#page 62
from __future__ import division
def f(x,y):
return (3*y*x**2+7)/10
def g(x,y):
return (y**2+4)/5
h=0.0001
x0=0.5
y0=0.5
f1_dx=(f(x0+h,y0)-f(x0,y0))/h
f1_dy=(f(x0,y0+h)-f(x0,y0))/h
g1_dx=(g(x0+h,y0)-g(x0,y0))/h
g1_dy=(g(x0+h,y0)-g(x0,y0))/h
if (f1_dx+f1_dy<1) and (g1_dx+g1_dy<1):
print "coditions for convergence is satisfied\n\n"
print "X \t Y\t\n\n"
for i in range(0,10):
X=(3*y0*x0**2+7)/10
Y=(y0**2+4)/5
print "%f\t %f\t\n" %(X,Y)
x0=X
y0=Y
print "\n\n CONVERGENCE AT (1 1) IS OBVIOUS"
#newton raphson metho
#example 2.32
#page 65
import numpy
def f(x,y):
return 3*y*x**2-10*x+7
def g(y):
return y**2-5*y+4
hh=0.0001
x0=0.5
y0=0.5 #initial values
f0=f(x0,y0)
g0=g(y0)
df_dx=(f(x0+hh,y0)-f(x0,y0))/hh
df_dy=(f(x0,y0+hh)-f(x0,y0))/hh
dg_dx=(g(y0)-g(y0))/hh
dg_dy=(g(y0+hh)-g(y0))/hh
d=[[df_dx,df_dy],[dg_dx,dg_dy]]
D1=numpy.linalg.det(d)
dd=[[-f0,df_dy],[-g0,dg_dy]]
h=numpy.linalg.det(dd)/D1
ddd=[[df_dx,-f0],[dg_dx,-g0]]
k=numpy.linalg.det(ddd)/D1;
x1=x0+h
y1=y0+k
f0=f(x1,y1)
g0=g(y1)
df_dx=(f(x1+hh,y1)-f(x1,y1))/hh
df_dy=(f(x1,y1+hh)-f(x1,y1))/hh
dg_dx=(g(y1)-g(y1))/hh
dg_dy=(g(y1+hh)-g(y1))/hh
dddd=[[df_dx,df_dy],[dg_dx,dg_dy]]
D2=numpy.linalg.det(dddd)
ddddd=[[-f0,df_dy],[-g0,dg_dy]]
h=numpy.linalg.det(ddddd)/D2
d6=[[df_dx,-f0],[dg_dx,-g0]]
k=numpy.linalg.det(d6)/D2
x2=x1+h
y2=y1+k
print " the roots of the equation are x2=%f and y2=%f" %(x2,y2)
#newton raphson method
#example 2.33
#page 66
import math
import numpy
def f(x,y):
return x**2+y**2-1
def g(x,y):
return y-x**2
hh=0.0001
x0=0.7071
y0=0.7071 #initial values
f0=f(x0,y0)
g0=g(x0,y0)
df_dx=(f(x0+hh,y0)-f(x0,y0))/hh
df_dy=(f(x0,y0+hh)-f(x0,y0))/hh
dg_dx=(g(x0+hh,y0)-g(x0,y0))/hh
dg_dy=(g(x0,y0+hh)-g(x0,y0))/hh
D1=numpy.linalg.det([[df_dx,df_dy],[dg_dx,dg_dy]])
h=numpy.linalg.det([[-f0,df_dy],[-g0,dg_dy]])/D1
k=numpy.linalg.det([[df_dx,-f0],[dg_dx,-g0]])/D1
x1=x0+h
y1=y0+k
f0=f(x1,y1)
g0=g(x1,y1)
df_dx=(f(x1+hh,y1)-f(x1,y1))/hh
df_dy=(f(x1,y1+hh)-f(x1,y1))/hh
dg_dx=(g(x1+hh,y1)-g(x1,y1))/hh
dg_dy=(g(x1,y1+hh)-g(x1,y1))/hh
D2=numpy.linalg.det([[df_dx,df_dy],[dg_dx,dg_dy]])
h=numpy.linalg.det([[-f0,df_dy],[-g0,dg_dy]])/D2
k=numpy.linalg.det([[df_dx,-f0],[dg_dx,-g0]])/D2
x2=x1+h
y2=y1+k
print "the roots of the equation are x2=%f and y2=%f " %(x2,y2)
#newton raphson method
#example 2.34
#page 67
import math
import numpy
def f(x,y):
return math.sin(x)-y+0.9793
def g(x,y):
return math.cos(y)-x+0.6703
hh=0.0001
x0=0.5
y0=1.5 #initial values
f0=f(x0,y0)
g0=g(x0,y0)
df_dx=(f(x0+hh,y0)-f(x0,y0))/hh
df_dy=(f(x0,y0+hh)-f(x0,y0))/hh
dg_dx=(g(x0+hh,y0)-g(x0,y0))/hh
dg_dy=(g(x0,y0+hh)-g(x0,y0))/hh
d1=[[df_dx,df_dy],[dg_dx,dg_dy]]
D1=numpy.linalg.det(d1)
d2=[[-f0,df_dy],[-g0,dg_dy]]
h=numpy.linalg.det(d2)/D1
d3=[[df_dx,-f0],[dg_dx,-g0]]
k=numpy.linalg.det(d3)/D1
x1=x0+h
y1=y0+k
f0=f(x1,y1)
g0=g(x1,y1)
df_dx=(f(x1+hh,y1)-f(x1,y1))/hh
df_dy=(f(x1,y1+hh)-f(x1,y1))/hh
dg_dx=(g(x1+hh,y1)-g(x1,y1))/hh
dg_dy=(g(x1,y1+hh)-g(x1,y1))/hh
d4=[[df_dx,df_dy],[dg_dx,dg_dy]]
D2=numpy.linalg.det(d4)
h=numpy.linalg.det([[-f0,df_dy],[-g0,dg_dy]])/D2
k=numpy.linalg.det([[df_dx,-f0],[dg_dx,-g0]])/D2
x2=x1+h
y2=y1+k
print "the roots of the equation are x2=%0.4f and y2=%0.4f" %(x2,y2)