Xsinx.f90, Example of carefully programming f(x) = x - sinx. Chapter 3: Locating Roots of Equations. Bisection.f90, Bisection method.
Let’s understand the secant method in numerical analysis and learn how to implement secant method in C programming with an explanation, output, advantages, disadvantages and much more.
What is Secant Method?
The secant method is a root-finding method that uses a succession of the roots of secant lines to find a better approximation of root.
The secant method algorithm is a root bracketing method and is the most efficient method of finding the root of a function. It requires two initial guesses which are the start and end interval points. This method is very similar to the Regula Falsi method.
The secant method is a Quasi-Newton method and it is seen to be faster in execution in some scenarios as compared to Newton-Raphson method and the False Position Method well.
The secant algorithm does not ensure convergence. The rate of convergence of secant method algorithm is 1.618, which is really fast comparatively. The order of convergence of secant method is superlinear.
The equation used in the following secant method c programs are as follows:
- x3 – 5x + 3
- x3 – 3x – 8
Secant Method Formula
Secant Method Algorithm
2 4 6 8 | Step1:Assume an equationf(x)=0(which must be defined) Step3:Do wherei=1,2,3,4,..... Step4:Evaluated Result |
Convergence criteria for Secant Method
- Fixing apriori the total number of iterations (limit).
- Testing the condition (xi+1−xi), is less than some tolerance limit (epsilon).
Advantages
- It evaluates one function at every iteration as compared with Newton’s method which evaluates two.
- The convergence rate is very fast.
- The secant method rule does not use the derivatives of a function.
Disadvantages
- The convergence may not always happen.
- Newton’s method generalizes much efficiently to new methods for solving simultaneous systems of nonlinear equations as compared to the Secant method.
Note: This secant method in C programming is compiled with GNU GCC compiler using CodeLite IDE on Microsoft Windows 10 operating system.
Method 1: C Program For Secant Method using Do While Loop
2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 | { floatstart_interval,end_interval,midpoint,allowed_error; printf('nEnter Value for Start Interval Point:t'); printf('nEnter Value for End Interval Point:t'); printf('nEnter Value for Allowed Error:t'); printf('nEnter Limit for Iterations:t'); do if(equation(start_interval)equation(end_interval)) printf('Please enter different start and end intervalsn'); } midpoint=(start_interval *equation(end_interval)-end_interval *equation(start_interval))/(equation(end_interval)-equation(start_interval)); end_interval=midpoint; printf('nIteration: %dtRoot: %fn',count,midpoint); if(countlimit) break; }while(fabs(equation(midpoint))>allowed_error); return0; { } |
Output
Method 2: Implement Secant Method in C Programming using While Loop
2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 | #include<math.h> #define EQUATION(x) (x)*(x)*(x) - 3*(x) - 8 intmain() floatinterval_start,interval_end,midpoint; printf('nEnter Value for Start Interval Point:t'); printf('nEnter Value for End Interval Point:t'); printf('nEnter value for Allowed Error:t'); printf('n'); printf('nStartttEndttMidPointtfunction(Start)tfunction(End)'); while(temp>allowed_error) a=EQUATION(interval_start); midpoint=interval_end-((b *(interval_end-interval_start))/(b-a)); printf('n%ft%ft%ft%ft%f',interval_start,interval_end,midpoint,a,b); interval_end=midpoint; { } { } printf('nnRoot of the Equation: %fn',midpoint); |
Output
If you have any doubts about the implementation of Secant method in C programming, let us know about it in the comment section. Find more about it on WikiPedia.
NUMERICAL METHODS C PROGRAMS |
---|
Newton-Raphson Method C Program |
Weddle’s Rule Algorithm C Program |
Euler’s Method C Program |
Bisection Method C Program |
Gauss Seidel Method C Program |
Simpson’s 3/8th Rule C Program |
Picard’s Method C Program |
Regula Falsi Method C Program |
Bisection Method Algorithm and Flowchart |
Simpson’s 1/3rd Rule C Program |
Trapezoidal Rule C Program |
The root of a one dimensional equation is the valueof x for which the equation y=f(x)=0.One might ask what this has to do with chemistry or science, etc.There are many examples of chemical problemsin which one needs to find the roots of equations.Some of the chemical examples include
·Weakacid/base equilibria; titrations
·Thevapor pressure equation, ln P = A + (B/T) + C ln T
·Equationof state problems
·Equilibriumcalculations
·Quantummechanics
·Kinetics
The philosophyof finding roots of equations is to use a series of approximations anditeration to solve the problems, in short, to use brute force.Why do we use this approach?A computer can do about 106operations/sec or about 3.6 x 109 operations/hr.A talented person might be able to do 1operation/sec or about 4 x 108 operations in a lifetime.In short, one hour of computer time is worthabout 1 human lifetime in mathematical problem solving.Thus the speed of computers makes the bruteforce approach an effective strategy.
The simplest method that is relatively foolproof isto plot the function and see where it crosses the x axis. Depending on howaccurately you need to know the solution, this method should give you theanswer you need relatively simply.
Newton’s method consists of four steps:
- Guess a value of the root, x.
- Calculate f(x) and the derivative f’(x), i.e., df(x)/dx
- Iterate
Graphically this can be depicted as:
One of the problems with the Newton method is thatone needs a good first guess, especially for multidimensional problems.The following illustrations indicate some ofthe difficulties that can occur when one has a poor guess and some pathologicalbehavior of the equation.
The following derivation demonstrates theconvergence of the Newton method.
What if the derivatives are not known?
Onecan approximate the derivatives using numerical methods (Lesson 11)Thus we have :
Ash goes to zero, this is the derivative.BUT the computer arithmetic won’t work as h goes to zero.Can get a sufficiently accurateapproximation in many cases by taking h = 10-3x.This should work for a seven significantfigure variable.
Two other methods of getting rootsare based on having two point initial guesses about the root, which hopefullybracket the true value.They are calledthe secant and regula falsi approaches.They are shown schematically below.
Inthe secant method, two points 1 and 2 are chosen.One draws a line between them.That line intersects the x axis and the projection of that point back onthe curve gives the point 3.One drawsa line from 2 to 3.That lineintersects the x axis at a point which is projected to give the point 4.The procedure is now repeated using thepoints 3 and 4.The generatingalgorithm for this method is
In the regula falsi method (method of false position) ,one draw a line between the two initial points 1 and 2, to give a intersectionpoint with the x axis.That point isprojected to the curve to give the point 3.A line is drawn between 3 and point 2.The point at which it intersects the x axis, is projected to give point4.Blah, blah, blah…The generatingalgorithm for the regula falsi method is
Problems with the regula falsi method include the effectof roundoff errors and slow converging nature of the algorithm.
Bisection
Bisection is the fail-safe method of root-finding.One picks two points that bracket theroot.One picks a value of x halfway inbetween the two point and test whether the root is to the right or left of thispoint.One then uses one of theoriginal endpoints and the new middle point to start the cycle over again.
Computer Methods
FORTRAN90
These methods are most transparently displayed in FORTRANwhere one can follow the program logic and the generating algorithms.The first program is xscrsho, a program toimplement the idiot’s method, i.e., a program to graph a function to show itsroots..The function being graphed isthe zero order Bessel function.
PROGRAMxscrsho
Cdriver forroutine scrsho
REAL bessj0
EXTERNALbessj0
write(*,*)'Graph of the Bessel Function J0:'
callscrsho(bessj0)
END
C(C) Copr.1986-92 Numerical Recipes Software
SUBROUTINEscrsho(fx)
INTEGERISCR,JSCR
REAL fx
EXTERNAL fx
PARAMETER(ISCR=60,JSCR=21)
INTEGERi,j,jz
REALdx,dyj,x,x1,x2,ybig,ysml,y(ISCR)
CHARACTER*1scr(ISCR,JSCR),blank,zero,yy,xx,ff
SAVEblank,zero,yy,xx,ff
DATAblank,zero,yy,xx,ff/' ','-','l','-','x'/
1continue
write (*,*) 'Enter x1,x2 (= to stop)'
read (*,*)x1,x2
if(x1.eq.x2)return
do 11j=1,JSCR
scr(1,j)=yy
scr(ISCR,j)=yy
11continue
do 13i=2,ISCR-1
scr(i,1)=xx
scr(i,JSCR)=xx
do 12j=2,JSCR-1
scr(i,j)=blank
12continue
13continue
dx=(x2-x1)/(ISCR-1)
ybig=0.
ysml=ybig
do 14i=1,ISCR
y(i)=fx(x)
if(y(i).lt.ysml) ysml=y(i)
if(y(i).gt.ybig) ybig=y(i)
x=x+dx
14continue
if(ybig.eq.ysml) ybig=ysml+1.
dyj=(JSCR-1)/(ybig-ysml)
jz=1-ysml*dyj
do 15i=1,ISCR
scr(i,jz)=zero
j=1+(y(i)-ysml)*dyj
scr(i,j)=ff
15continue
write(*,'(1x,1pe10.3,1x,80a1)') ybig,(scr(i,JSCR),i=1,ISCR)
do 16j=JSCR-1,2,-1
write(*,'(12x,80a1)') (scr(i,j),i=1,ISCR)
16continue
write(*,'(1x,1pe10.3,1x,80a1)') ysml,(scr(i,1),i=1,ISCR)
write(*,'(12x,1pe10.3,40x,e10.3)') x1,x2
goto 1
END
C(C) Copr.1986-92 Numerical Recipes Software
FUNCTIONbessj0(x)
REAL bessj0,x
REAL ax,xx,z
DOUBLEPRECISION p1,p2,p3,p4,p5,q1,q2,q3,q4,q5,r1,r2,r3,r4,r5,r6,
*s1,s2,s3,s4,s5,s6,y
SAVEp1,p2,p3,p4,p5,q1,q2,q3,q4,q5,r1,r2,r3,r4,r5,r6,s1,s2,s3,s4,
*s5,s6
DATAp1,p2,p3,p4,p5/1.d0,-.1098628627d-2,.2734510407d-4,
*-.2073370639d-5,.2093887211d-6/, q1,q2,q3,q4,q5/-.1562499995d-1,
*.1430488765d-3,-.6911147651d-5,.7621095161d-6,-.934945152d-7/
DATAr1,r2,r3,r4,r5,r6/57568490574.d0,-13362590354.d0,
*651619640.7d0,-11214424.18d0,77392.33017d0,-184.9052456d0/,s1,s2,
*s3,s4,s5,s6/57568490411.d0,1029532985.d0,9494680.718d0,
*59272.64853d0,267.8532712d0,1.d0/
if(abs(x).lt.8.)then
y=x**2
bessj0=(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))))/(s1+y*(s2+y*(s3+y*
*(s4+y*(s5+y*s6)))))
else
ax=abs(x)
z=8./ax
y=z**2
xx=ax-.785398164
bessj0=sqrt(.636619772/ax)*(cos(xx)*(p1+y*(p2+y*(p3+y*(p4+y*
*p5))))-z*sin(xx)*(q1+y*(q2+y*(q3+y*(q4+y*q5)))))
endif
return
END
C(C) Copr.1986-92 Numerical Recipes Software
The second program implements a solution using Newton’smethod.
PROGRAMxrtnewt
Cdriver forroutine rtnewt
INTEGERN,NBMAX
REAL X1,X2
PARAMETER(N=100,NBMAX=20,X1=1.0,X2=50.0)
INTEGER i,nb
REALbessj0,rtnewt,root,xacc,xb1(NBMAX),xb2(NBMAX)
EXTERNALfuncd,bessj0
nb=NBMAX
callzbrak(bessj0,X1,X2,N,xb1,xb2,nb)
write(*,'(/1x,a)') 'Roots of BESSJ0:'
write(*,'(/1x,t19,a,t31,a/)') 'x','F(x)'
do 11 i=1,nb
xacc=(1.0e-6)*(xb1(i)+xb2(i))/2.0
root=rtnewt(funcd,xb1(i),xb2(i),xacc)
write(*,'(1x,a,i2,2x,f12.6,e16.4)') 'Root ',i,root,bessj0(root)
11continue
END
SUBROUTINEfuncd(x,fn,df)
REALbessj0,bessj1,df,fn,x
fn=bessj0(x)
df=-bessj1(x)
return
END
C(C) Copr.1986-92 Numerical Recipes Software
FUNCTIONrtnewt(funcd,x1,x2,xacc)
INTEGER JMAX
REALrtnewt,x1,x2,xacc
EXTERNALfuncd
PARAMETER(JMAX=20)
REAL df,dx,f
rtnewt=.5*(x1+x2)
do 11j=1,JMAX
callfuncd(rtnewt,f,df)
dx=f/df
rtnewt=rtnewt-dx
if((x1-rtnewt)*(rtnewt-x2).lt.0.)pause
*'rtnewtjumped out of brackets'
if(abs(dx).lt.xacc) return
11continue
pause 'rtnewtexceeded maximum iterations'
END
C(C) Copr.1986-92 Numerical Recipes Software
FUNCTIONbessj0(x)
REAL bessj0,x
REAL ax,xx,z
DOUBLEPRECISION p1,p2,p3,p4,p5,q1,q2,q3,q4,q5,r1,r2,r3,r4,r5,r6,
*s1,s2,s3,s4,s5,s6,y
SAVEp1,p2,p3,p4,p5,q1,q2,q3,q4,q5,r1,r2,r3,r4,r5,r6,s1,s2,s3,s4,
*s5,s6
DATAp1,p2,p3,p4,p5/1.d0,-.1098628627d-2,.2734510407d-4,
*-.2073370639d-5,.2093887211d-6/, q1,q2,q3,q4,q5/-.1562499995d-1,
*.1430488765d-3,-.6911147651d-5,.7621095161d-6,-.934945152d-7/
DATAr1,r2,r3,r4,r5,r6/57568490574.d0,-13362590354.d0,
*651619640.7d0,-11214424.18d0,77392.33017d0,-184.9052456d0/,s1,s2,
*s3,s4,s5,s6/57568490411.d0,1029532985.d0,9494680.718d0,
*59272.64853d0,267.8532712d0,1.d0/
if(abs(x).lt.8.)then
y=x**2
bessj0=(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))))/(s1+y*(s2+y*(s3+y*
*(s4+y*(s5+y*s6)))))
else
ax=abs(x)
z=8./ax
y=z**2
bessj0=sqrt(.636619772/ax)*(cos(xx)*(p1+y*(p2+y*(p3+y*(p4+y*
*p5))))-z*sin(xx)*(q1+y*(q2+y*(q3+y*(q4+y*q5)))))
endif
return
END
C(C) Copr.1986-92 Numerical Recipes Software
FUNCTIONbessj1(x)
REAL ax,xx,z
DOUBLEPRECISION p1,p2,p3,p4,p5,q1,q2,q3,q4,q5,r1,r2,r3,r4,r5,r6,
*s1,s2,s3,s4,s5,s6,y
SAVEp1,p2,p3,p4,p5,q1,q2,q3,q4,q5,r1,r2,r3,r4,r5,r6,s1,s2,s3,s4,
*s5,s6
DATAr1,r2,r3,r4,r5,r6/72362614232.d0,-7895059235.d0,
*242396853.1d0,-2972611.439d0,15704.48260d0,-30.16036606d0/,s1,s2,
*s3,s4,s5,s6/144725228442.d0,2300535178.d0,18583304.74d0,
*99447.43394d0,376.9991397d0,1.d0/
DATAp1,p2,p3,p4,p5/1.d0,.183105d-2,-.3516396496d-4,
*-.2002690873d-3,.8449199096d-5,-.88228987d-6,.105787412d-6/
if(abs(x).lt.8.)then
y=x**2
bessj1=x*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))))/(s1+y*(s2+y*(s3+
*y*(s4+y*(s5+y*s6)))))
else
ax=abs(x)
z=8./ax
y=z**2
xx=ax-2.356194491
bessj1=sqrt(.636619772/ax)*(cos(xx)*(p1+y*(p2+y*(p3+y*(p4+y*
*p5))))-z*sin(xx)*(q1+y*(q2+y*(q3+y*(q4+y*q5)))))*sign(1.,x)
endif
END
C(C) Copr.1986-92 Numerical Recipes Software
SUBROUTINEzbrak(fx,x1,x2,n,xb1,xb2,nb)
INTEGER n,nb
REALx1,x2,xb1(nb),xb2(nb),fx
EXTERNAL fx
INTEGER i,nbb
REALdx,fc,fp,x
nbb=0
x=x1
fp=fx(x)
do 11 i=1,n
x=x+dx
fc=fx(x)
if(fc*fp.lt.0.) then
nbb=nbb+1
xb1(nbb)=x-dx
xb2(nbb)=x
if(nbb.eq.nb)goto 1
endif
fp=fc
11continue
1continue
nb=nbb
return
END
C(C) Copr.1986-92 Numerical Recipes Software
The next program implements the secant method of findingroots.
PROGRAMxrtsec
Cdriver forroutine rtsec
INTEGERN,NBMAX
REAL X1,X2
INTEGER i,nb
REALbessj0,rtsec,root,xacc,xb1(NBMAX),xb2(NBMAX)
EXTERNALbessj0
nb=NBMAX
callzbrak(bessj0,X1,X2,N,xb1,xb2,nb)
write(*,'(/1x,a)') 'Roots of BESSJ0:'
write(*,'(/1x,t19,a,t31,a/)')'x','F(x)'
do 11 i=1,nb
xacc=(1.0e-6)*(xb1(i)+xb2(i))/2.0
root=rtsec(bessj0,xb1(i),xb2(i),xacc)
write(*,'(1x,a,i2,2x,f12.6,e16.4)') 'Root ',i,root,bessj0(root)
11continue
END
C(C) Copr. 1986-92Numerical Recipes Software
FUNCTIONrtsec(func,x1,x2,xacc)
INTEGER MAXIT
REALrtsec,x1,x2,xacc,func
EXTERNAL func
PARAMETER(MAXIT=30)
INTEGER j
REALdx,f,fl,swap,xl
fl=func(x1)
f=func(x2)
if(abs(fl).lt.abs(f))then
rtsec=x1
xl=x2
swap=fl
fl=f
f=swap
else
xl=x1
rtsec=x2
endif
do 11j=1,MAXIT
dx=(xl-rtsec)*f/(f-fl)
xl=rtsec
fl=f
rtsec=rtsec+dx
f=func(rtsec)
if(abs(dx).lt.xacc.or.f.eq.0.)return
11continue
pause 'rtsecexceed maximum iterations'
END
C(C) Copr.1986-92 Numerical Recipes Software
SUBROUTINEzbrak(fx,x1,x2,n,xb1,xb2,nb)
INTEGER n,nb
EXTERNAL fx
INTEGER i,nbb
REALdx,fc,fp,x
nbb=0
x=x1
dx=(x2-x1)/n
fp=fx(x)
do 11 i=1,n
x=x+dx
fc=fx(x)
if(fc*fp.lt.0.) then
nbb=nbb+1
xb2(nbb)=x
if(nbb.eq.nb)goto 1
endif
fp=fc
11continue
1continue
nb=nbb
return
END
C(C) Copr.1986-92 Numerical Recipes Software
FUNCTIONbessj0(x)
REAL bessj0,x
REAL ax,xx,z
DOUBLEPRECISION p1,p2,p3,p4,p5,q1,q2,q3,q4,q5,r1,r2,r3,r4,r5,r6,
*s1,s2,s3,s4,s5,s6,y
SAVEp1,p2,p3,p4,p5,q1,q2,q3,q4,q5,r1,r2,r3,r4,r5,r6,s1,s2,s3,s4,
*s5,s6
DATAp1,p2,p3,p4,p5/1.d0,-.1098628627d-2,.2734510407d-4,
*-.2073370639d-5,.2093887211d-6/, q1,q2,q3,q4,q5/-.1562499995d-1,
*.1430488765d-3,-.6911147651d-5,.7621095161d-6,-.934945152d-7/
DATAr1,r2,r3,r4,r5,r6/57568490574.d0,-13362590354.d0,
*651619640.7d0,-11214424.18d0,77392.33017d0,-184.9052456d0/,s1,s2,
*s3,s4,s5,s6/57568490411.d0,1029532985.d0,9494680.718d0,
*59272.64853d0,267.8532712d0,1.d0/
if(abs(x).lt.8.)then
y=x**2
bessj0=(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))))/(s1+y*(s2+y*(s3+y*
*(s4+y*(s5+y*s6)))))
else
ax=abs(x)
z=8./ax
y=z**2
![For For](/uploads/1/2/6/3/126364695/798837569.jpg)
xx=ax-.785398164
bessj0=sqrt(.636619772/ax)*(cos(xx)*(p1+y*(p2+y*(p3+y*(p4+y*
*p5))))-z*sin(xx)*(q1+y*(q2+y*(q3+y*(q4+y*q5)))))
endif
return
END
C(C) Copr.1986-92 Numerical Recipes Software
The next program implements regula falsi method.
PROGRAMxrtflsp
Cdriver forroutine rtflsp
INTEGERN,NBMAX
REAL X1,X2
PARAMETER(N=100,NBMAX=20,X1=1.0,X2=50.0)
INTEGER i,nb
REALbessj0,rtflsp,root,xacc,xb1(NBMAX),xb2(NBMAX)
EXTERNALbessj0
nb=NBMAX
callzbrak(bessj0,X1,X2,N,xb1,xb2,nb)
write(*,'(/1x,a)') 'Roots of BESSJ0:'
write(*,'(/1x,t19,a,t31,a/)') 'x','F(x)'
do 11 i=1,nb
root=rtflsp(bessj0,xb1(i),xb2(i),xacc)
write(*,'(1x,a,i2,2x,f12.6,e16.4)') 'Root ',i,root,bessj0(root)
11continue
END
C(C) Copr.1986-92 Numerical Recipes Software
FUNCTIONrtflsp(func,x1,x2,xacc)
INTEGER MAXIT
REALrtflsp,x1,x2,xacc,func
EXTERNAL func
PARAMETER(MAXIT=30)
INTEGER j
REALdel,dx,f,fh,fl,swap,xh,xl
fl=func(x1)
fh=func(x2)
if(fl*fh.gt.0.) pause 'root must be bracketed in rtflsp'
if(fl.lt.0.)then
xl=x1
xh=x2
else
xl=x2
xh=x1
swap=fl
fl=fh
fh=swap
endif
dx=xh-xl
do 11j=1,MAXIT
rtflsp=xl+dx*fl/(fl-fh)
f=func(rtflsp)
del=xl-rtflsp
xl=rtflsp
fl=f
else
del=xh-rtflsp
xh=rtflsp
fh=f
endif
dx=xh-xl
if(abs(del).lt.xacc.or.f.eq.0.)return
11continue
pause 'rtflspexceed maximum iterations'
END
C(C) Copr.1986-92 Numerical Recipes Software
FUNCTIONbessj0(x)
REAL bessj0,x
REAL ax,xx,z
DOUBLEPRECISION p1,p2,p3,p4,p5,q1,q2,q3,q4,q5,r1,r2,r3,r4,r5,r6,
*s1,s2,s3,s4,s5,s6,y
*s5,s6
DATAp1,p2,p3,p4,p5/1.d0,-.1098628627d-2,.2734510407d-4,
*-.2073370639d-5,.2093887211d-6/, q1,q2,q3,q4,q5/-.1562499995d-1,
*.1430488765d-3,-.6911147651d-5,.7621095161d-6,-.934945152d-7/
DATAr1,r2,r3,r4,r5,r6/57568490574.d0,-13362590354.d0,
*651619640.7d0,-11214424.18d0,77392.33017d0,-184.9052456d0/,s1,s2,
*s3,s4,s5,s6/57568490411.d0,1029532985.d0,9494680.718d0,
*59272.64853d0,267.8532712d0,1.d0/
if(abs(x).lt.8.)then
y=x**2
bessj0=(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))))/(s1+y*(s2+y*(s3+y*
*(s4+y*(s5+y*s6)))))
else
ax=abs(x)
z=8./ax
y=z**2
xx=ax-.785398164
bessj0=sqrt(.636619772/ax)*(cos(xx)*(p1+y*(p2+y*(p3+y*(p4+y*
*p5))))-z*sin(xx)*(q1+y*(q2+y*(q3+y*(q4+y*q5)))))
endif
return
END
C(C) Copr.1986-92 Numerical Recipes Software
SUBROUTINEzbrak(fx,x1,x2,n,xb1,xb2,nb)
INTEGER n,nb
REALx1,x2,xb1(nb),xb2(nb),fx
EXTERNAL fx
INTEGER i,nbb
REALdx,fc,fp,x
nbb=0
x=x1
dx=(x2-x1)/n
fp=fx(x)
do 11 i=1,n
x=x+dx
fc=fx(x)
if(fc*fp.lt.0.) then
nbb=nbb+1
xb1(nbb)=x-dx
xb2(nbb)=x
if(nbb.eq.nb)goto 1
endif
fp=fc
11continue
1continue
nb=nbb
return
END
C(C) Copr.1986-92 Numerical Recipes Software
The final program implements the bisection method.
PROGRAMxrtbis
Cdriver forroutine rtbis
INTEGERN,NBMAX
REAL X1,X2
PARAMETER(N=100,NBMAX=20,X1=1.0,X2=50.0)
INTEGER i,nb
REALbessj0,rtbis,root,xacc,xb1(NBMAX),xb2(NBMAX)
EXTERNALbessj0
nb=NBMAX
callzbrak(bessj0,X1,X2,N,xb1,xb2,nb)
write(*,'(/1x,a)') 'Roots of BESSJ0:'
write(*,'(/1x,t19,a,t31,a/)') 'x','F(x)'
do 11 i=1,nb
xacc=(1.0e-6)*(xb1(i)+xb2(i))/2.0
root=rtbis(bessj0,xb1(i),xb2(i),xacc)
write(*,'(1x,a,i2,2x,f12.6,e16.4)')'Root ',i,root,bessj0(root)
11continue
END
C(C) Copr.1986-92 Numerical Recipes Software
FUNCTIONrtbis(func,x1,x2,xacc)
INTEGER JMAX
REALrtbis,x1,x2,xacc,func
EXTERNAL func
PARAMETER(JMAX=40)
INTEGER j
REALdx,f,fmid,xmid
fmid=func(x2)
f=func(x1)
if(f*fmid.ge.0.) pause 'root must be bracketed in rtbis'
if(f.lt.0.)then
rtbis=x1
dx=x2-x1
else
rtbis=x2
dx=x1-x2
endif
do 11j=1,JMAX
dx=dx*.5
xmid=rtbis+dx
fmid=func(xmid)
if(fmid.le.0.)rtbis=xmid
if(abs(dx).lt.xacc .or. fmid.eq.0.) return
11continue
pause 'toomany bisections in rtbis'
END
C(C) Copr.1986-92 Numerical Recipes Software
FUNCTIONbessj0(x)
REAL bessj0,x
REAL ax,xx,z
DOUBLEPRECISION p1,p2,p3,p4,p5,q1,q2,q3,q4,q5,r1,r2,r3,r4,r5,r6,
*s1,s2,s3,s4,s5,s6,y
SAVEp1,p2,p3,p4,p5,q1,q2,q3,q4,q5,r1,r2,r3,r4,r5,r6,s1,s2,s3,s4,
*s5,s6
DATAp1,p2,p3,p4,p5/1.d0,-.1098628627d-2,.2734510407d-4,
*-.2073370639d-5,.2093887211d-6/, q1,q2,q3,q4,q5/-.1562499995d-1,
*.1430488765d-3,-.6911147651d-5,.7621095161d-6,-.934945152d-7/
DATAr1,r2,r3,r4,r5,r6/57568490574.d0,-13362590354.d0,
*651619640.7d0,-11214424.18d0,77392.33017d0,-184.9052456d0/,s1,s2,
*s3,s4,s5,s6/57568490411.d0,1029532985.d0,9494680.718d0,
*59272.64853d0,267.8532712d0,1.d0/
if(abs(x).lt.8.)then
y=x**2
*(s4+y*(s5+y*s6)))))
else
ax=abs(x)
z=8./ax
y=z**2
xx=ax-.785398164
bessj0=sqrt(.636619772/ax)*(cos(xx)*(p1+y*(p2+y*(p3+y*(p4+y*
*p5))))-z*sin(xx)*(q1+y*(q2+y*(q3+y*(q4+y*q5)))))
endif
return
END
C(C) Copr.1986-92 Numerical Recipes Software
SUBROUTINEzbrak(fx,x1,x2,n,xb1,xb2,nb)
INTEGER n,nb
REALx1,x2,xb1(nb),xb2(nb),fx
EXTERNAL fx
REALdx,fc,fp,x
nbb=0
x=x1
dx=(x2-x1)/n
fp=fx(x)
do 11 i=1,n
x=x+dx
fc=fx(x)
if(fc*fp.lt.0.) then
nbb=nbb+1
xb1(nbb)=x-dx
xb2(nbb)=x
if(nbb.eq.nb)goto1
endif
fp=fc
11continue
1continue
nb=nbb
return
END
C(C) Copr.1986-92 Numerical Recipes Software
Mathematica
In Mathematica, root finding is done in a black boxmanner, using canned functions. Thethree relevant functions are Solve, which finds all roots; Nsolve, which givesa numerical approximation to all rootsand FindRootwhich searches forthe root of an equation in the vicinity of an initial guess.These functions are summarized in the tablebelow.
The following example shows the use of the FindRootfunction
Excel
Excel has two black box functions that can be used to findroots, Goal Seek and Solver.Thefollowing example uses Goal Seek to solve cubic equations.
Consider a cubic equation of the ax3 + bx2+cx +d =0.The cells b3-b7 contain thecoefficients a,b,c,d.Thesecoefficients have been set to 2,1,-246 and 360.Initially the cells e4-e6 are set equal to the formula=a*d4^3+b*d4^2+c*d4+d.Move to cell d4and use Goal Seek to find the first solution by varying d4 until e4 has a zerovalue.Move to e5 and use Goal Seek tovary d5 until e5 is zero.Repeat withe6 and d6.The cells d4:d6 shouldcontain the solutions 10,1.5,-12.
To do the same procedure using Solver, restore thespreadsheet to the original values you had before using Goal Seek.Move to cell e4 and select Solver from theTools menu.Make the target cell $e$4.set value =0 and select this option.Select $d$4 as the cell to vary.Click Solve.Repeat for e5, d5 and e6 and d6.