FORTRAN PROGRAMS FOR SOLVIMG NUMERICAL PROBLEMS Designed by T K Rajan Seln. Grade Lecturer in Mathematics Govt Victoria College, Palakkad
BISECTION METHOD
PROGRAM bs c To find zero of an equation by bisection method write(*,40) 40 format(1x,'To find a real root of an equation using Bisection',\) write(*,50) 50 format(' method.') write(*,*) write(*,*)'Enter numbers between which the root is to be found:' read(*,*)a,b IF((f(a).GT.0).AND.(f(b).LT.0))THEN w=a a=b b=w ENDIF write(*,*)'Input error value:' read(*,*)e 20 c=(a+b)/2 IF(abs(f(c)).LT.e) GOTO 30 IF(f(c).LT.0)THEN a=c ELSE b=c ENDIF GOTO 20 30 write(*,10)c 10 format(1x,'The approximate root of the function is:',F8.3) STOP END function f(c) f=c*c-5*c+6 return END
GAUSS ELIMINATION
PROGRAM gelmn c Gauss elimination method to solve a linear system of eqns dimension a(30,30),x(30) write(*,10) 10 format(1x,'To solve a linear system of equations using Gauss ',\) write(*,20) 20 format('Elimination method with pivoting(using subroutine).') write(*,*)'Enter the number of variables:' read(*,*)n write(*,*)'Enter the coefficients in the equations:' read(*,*)((a(i,j),j=1,n+1),i=1,n)
DO 30 k=1,n-1 call pivot(a,k,n) DO 90 i=k+1,n u=a(i,k)/a(k,k) DO 100 j=k,n+1 a(i,j)=a(i,j)-u*a(k,j) 100 continue 90 continue 30 continue IF(abs(a(n,n)).LE.(.00001))THEN write(*,*)'Ill conditioned equations.' STOP ENDIF x(n)=a(n,n+1)/a(n,n) DO 60 i=n-1,1,-1 sum=0 DO 70 j=i+1,n sum=sum+a(i,j)*x(j) 70 continue x(i)=(a(i,n+1)-sum)/a(i,i) 60 continue write(*,*)'Values of the variables are as follows:' DO 80 i=1,n write(*,110)x(i) 110 format(1x,F10.3) 80 continue STOP END subroutine pivot(a,k,n) dimension a(30,30) real mx integer p,q mx=abs(a(k,k))
p=k DO 40 m=k+1,n IF(abs(a(m,k)).GT.mx)THEN mx=abs(a(m,k)) p=m ENDIF 40 continue IF(mx.LE.(.00001))THEN write(*,*)'Ill-conditioned equations.' STOP ENDIF DO 50 q=k,n+1 temp=a(k,q) a(k,q)=a(p,q) a(p,q)=temp 50 continue return END
EULER’ METHOD
PROGRAM eu c Euler's method to solve a first order de dimension x(50),y(50) write(*,*)'To solve a differential equation using Euler method.' write(*,*)'Enter initial values:' read(*,*)a,c write(*,*)'Enter value at which result is to be found:' read(*,*)b write(*,*)'Enter the number of subintervals required:' read(*,*)n x(1)=a y(1)=c h=(b-a)/n write(*,*)'The function values at each step is given below.' write(*,*) write(*,*)' x Function value' write(*,10)x(1),y(1) 10 format(1x,F10.4,5x,F10.4) DO 20 i=1,n x(i+1)=x(i)+h y(i+1)=y(i)+h*f(x(i),y(i)) write(*,10)x(i+1),y(i+1) 20 continue STOP END
function f(x,y) f=(y-x)/(y+x) return END
REGULA FALSI METHOD
PROGRAM regfal c Regula falsi Method to a find a root of an eqn write(*,40) 40 format(1x,'To find the approximate root of a function using ') write(*,50) 50 format(\,'Method of False Position.') write(*,*)'Enter the numbers between which root is to be found:' read(*,*)a,b IF((f(a).GT.0).AND.(f(b).LT.0))THEN w=a a=b b=w ENDIF write(*,*)'Input error value:' read(*,*)e 20 c=(a*f(b)-b*f(a))/(f(b)-f(a)) IF(abs(f(c)).LT.e) GOTO 10 IF(f(c).LT.0)THEN a=c ELSE b=c ENDIF GOTO 20 10 write(*,30)c 30 format(1x,'The approximate root of the function is:',F8.3) STOP END function f(c) f=c**3-9*c+1 return END
GCD OF A SEQUENCE OF NUMBERS
PROGRAM ss integer x(30) write(*,*)'To find the gcd of a set of numbers.' write(*,*)'Enter number of numbers:' read(*,*)n write(*,*)'Enter the numbers:' read(*,*)(x(i),i=1,n) j=x(1) DO 10 i=2,n k=gg(j,x(i)) j=k 10 continue write(*,20)j 20 format(1x,' GCD= ',I4) STOP END function gg(a,b) integer a,b IF(a.GT.b)THEN j=a a=b b=j ENDIF 30 i=mod(b,a) IF(i.EQ.0)THEN gg=a return ELSE b=a a=i GOTO 30 ENDIF END
FIBONACCI SEQUENCE
PROGRAM fbnc c Fobonacci sequence integer f1,f2,f3 f1=1 f2=1 write(*,*)'To display fibonacci numbers.' write(*,*)'Enter the number of terms to be displayed:' read(*,*)n write(*,*)'The fibonacci sequence is : ' IF(n.EQ.1)THEN write(*,10)f1 10 format(1x,I6) ELSE write(*,10)f1 write(*,10)f2 DO 20 i=3,n f3=f1+f2 write(*,10)f3 f1=f2 f2=f3 20 continue ENDIF STOP END
GCD & LCM OF TWO NUMBERS
PROGRAM gcdlcm integer a,b,l write(*,*)'To find the gcd and lcm of two numbers.' write(*,*)'Enter 2 numbers:' read(*,*)a,b m=a n=b IF(a.GT.b)THEN j=a a=b b=j ENDIF 10 i=mod(b,a) IF(i.EQ.0)THEN write(*,20)a 20 format(1x,'GCD= ',I4) ELSE b=a a=i GOTO 10 ENDIF l=m*n/a write(*,30)l 30 format(1x,'LCM= ',I4) STOP END
LAGRANGE’S METHOD OF INTERPOLATION
PROGRAM ll c Lagrange's Method of interpolation dimension x(20),y(20) write(*,50) 50 format(1x,'To find the function value for a particular value ',\) write(*,60) 60 format('of x using Lagrange interpolation formula.') write(*,*) write(*,*)'Enter the number of known function values:' read(*,*)n write(*,*)'Enter value of x and corresponding function values:' DO 10 i=1,n read(*,*)x(i),y(i) 10 continue write(*,*)'Enter x for which function value is to be found:' read(*,*)a s=0 DO 20 i=1,n p=1 DO 30 j=1,n IF(i.NE.j) p=p*(a-x(j))/(x(i)-x(j)) 30 continue s=s+p*y(i) 20 continue write(*,*) write(*,40)a,s 40 format(1x,'The function value at x =',F7.2,' is:',F9.3) STOP END
NEWTON’S METHOD OF INTERPOLATION
PROGRAM newt c Newton's Mthod of interpolation dimension x(20),y(20,20) write(*,100) 100 format(1x,'To find the value of a function corresponding to a',\) write(*,110) 110 format(' value of x using Newtons divided difference method.') write(*,*) write(*,*)'Enter the number of known function values:' read(*,*)n write(*,*)'Enter the values of x and the corresponding function *values:' DO 10 i=1,n read(*,*)x(i),y(i,1) 10 continue write(*,*)'Enter the value of x for which function value is to be *found:' read(*,*)a k=0 DO 20 j=2,n k=k+1 DO 30 i=1,n-k y(i,j)=(y(i+1,j-1)-y(i,j-1))/(x(i+k)-x(i)) 30 continue 20 continue write(*,*)'The divided difference table is : ' write(*,*) DO 70 i=1,n write(*,90)x(i) DO 80 j=1,n-i+1 write(*,90)y(i,j) 90 format(1x,F10.3\) 80 continue write(*,*) write(*,*) 70 continue s=y(1,1) DO 40 j=2,n p=1 DO 50 i=1,j-1 p=p*(a-x(i)) 50 continue s=s+p*y(1,j) 40 continue write(*,60)s 60 format(1x,'The corresponding function value is:',F10.3) STOP END
NEWTON RAPHSON METHOD
PROGRAM nr c Newton Raphson Method to fin a root of an eqn integer count count=0 write(*,40) 40 format(1x,'To find a root of an equation using Newton Raphson',\) write(*,50) 50 format(' Method.') write(*,*)'Enter initial root and error value:' read(*,*)x,e 10 IF(g(x).EQ.0)THEN write(*,*)'Incorrect initial root.' STOP ENDIF y=x-(f(x)/g(x)) IF(abs(f(y)).LT.e) GOTO 20 count=count+1 IF(count.GT.500)THEN write(*,*)'An error has occured.' STOP ENDIF x=y GOTO 10 20 write(*,30)x 30 format(1x,'Approximate root= ',F10.3) STOP END function f(x) f=x*x-5*x+6 return END function g(x) g=2*x-5 return END
NUMERICAL DIFFERENTAION
PROGRAM ndff c Numerical Differentiation dimension x(20),y(20,20) write(*,70) 70 format(1x,'To find the value of derivative of a function at a',\) write(*,80) 80 format(' given value using Numerical Differentiation.') write(*,*) write(*,*)'Enter the number of known function values:' read(*,*)n write(*,*)'Enter values of x and corresponding function values:' read(*,*)(x(i),y(i,1),i=1,n) write(*,*)'Enter the value at which derivative is to be found:' read(*,*)a DO 10 j=2,n DO 10 i=1,n-j+1 y(i,j)=(y(i+1,j-1)-y(i,j-1))/(x(i+j-1)-x(i)) 10 continue v=y(1,2) DO 20 i=3,n s=0 DO 30 j=1,i-1 p=1 DO 40 k=1,i-1 IF(k.NE.j) p=p*(a-x(k)) 40 continue s=s+p 30 continue v=v+s*y(1,i) 20 continue write(*,50)a,v 50 format(1x,'The value of the derivative at x= ',F9.4,' is: ',F9.4) STOP END
PRODUCT OF TWO MATRICES
PROGRAM matsub c Product of two matrices dimension a(10,10),b(10,10),c(10,10) integer p,q write(*,*)'To find the product of two matrices.' write(*,*)'Enter the order of the first matrix:' read(*,*)m,n write(*,*)'Enter the order of the second matrix:' read(*,*)p,q IF(n.NE.p)THEN write(*,*)'The matrices are not conformable for multiplication.' STOP ENDIF call inpt(a,m,n) call inpt(b,p,q) call pro(a,b,c,m,n,q) call outpt(c,m,q) STOP END subroutine inpt(a,m,n) dimension a(10,10) write(*,10)m,n 10 format(1x,'Enter the matrix of order',I2,' *'I2,' :') read(*,*)((a(i,j),j=1,n),i=1,m) END subroutine pro(a,b,c,m,n,q) dimension a(10,10),b(10,10),c(10,10) integer q DO 20 i=1,m DO 30 j=1,q c(i,j)=0 DO 40 k=1,n c(i,j)=c(i,j)+a(i,k)*b(k,j) 40 continue 30 continue 20 continue END subroutine outpt(c,m,q) dimension c(10,10) integer q write(*,*)'The product matrix is: ' DO 50 i=1,m DO 60 j=1,q write(*,70)c(i,j) 70 format(3x,F6.2,4x\) 60 continue write(*,*) 50 continue END
RUNGE-KUTTA METHOD
PROGRAM rr c Rung- Kutta Method to solve a d.e. integer count real k1,k2 count=0 write(*,10) 10 format(1x,'To solve a differential equation using second ',\) write(*,20) 20 format('order Runge Kutta method.') write(*,*)'Enter the initial values:' read(*,*)x1,y1 write(*,*)'Enter value at which function value is to be found:' read(*,*)a write(*,*)'Enter the number of subintervals:' read(*,*)n h=(a-x1)/n write(*,*)'Values of x and corresponding function values are: ' write(*,*) write(*,30) 30 format(8x,' x ',5x,' f(x) ') write(*,40)x1,y1 40 format(1x,F10.4,F10.4) 50 k1=h*f(x1,y1) k2=h*f(x1+h,y1+k1) y2=y1+(k1+k2)/2 x2=x1+h write(*,40)x2,y2 count=count+1 IF(count.LT.n) THEN x1=x2 y1=y2 GOTO 50 ENDIF write(*,*) write(*,60)x2,y2 60 format(1x,'The value of the function at x=',F10.4,' is:',F10.4) STOP END function f(x,y) f=x+y return END
RUNGE-KUTTA METHOD OF ORDER 4
PROGRAM rrkk c Runge-Kutta method of order 4 integer count count=0 write(*,10) 10 format(1x,'To solve a differential equation using fourth ',\) write(*,20) 20 format('order Runge Kutta method.') write(*,*) write(*,*)'Enter the initial values:' read(*,*)x1,y1 write(*,*)'Enter value at which function value is to be found:' read(*,*)a write(*,*)'Enter the number of subintervals:' read(*,*)n h=(a-x1)/n write(*,*)'The function values are as follows: ' write(*,*) write(*,30) 30 format(8x,' x ',5x,' f(x) ') write(*,40)x1,y1 40 format(1x,F10.4,F10.4) 50 s1=f(x1,y1) s2=f(x1+h/2,y1+s1*h/2) s3=f(x1+h/2,y1+s2*h/2) s4=f(x1+h,y1+s3*h) y2=y1+(s1+2*s2+2*s3+s4)*h/6 x2=x1+h write(*,40)x2,y2 count=count+1 IF(count.LT.n)THEN x1=x2 y1=y2 GOTO 50 ENDIF write(*,*) write(*,60)x2,y2 60 format(1x,'The value at ',F10.4,' is :',F10.4) STOP END function f(x,y) f=x+y return END
SIMPSON’S RULE OF INTEGRATION
PROGRAM zz c Simpson's rule of integration integer count count=1 write(*,30) 30 format(1x,'To find the value of integral of a function using',\) write(*,40) 40 format(' Simpsons rule. ') write(*,*) write(*,*)'Enter the limits of integration:' read(*,*)a,b write(*,*)'Enter the number of subintervals(an even no.) :' read(*,*)n h=(b-a)/n x=a y=f(x) sum=y 10 IF(count.LT.n)THEN x=x+h y=f(x) count=count+1 IF(mod(count,2).EQ.0)THEN sum=sum+4*y ELSE sum=sum+2*y ENDIF GOTO 10 ENDIF x=x+h y=f(x) sum=sum+y sum=h*sum/3 write(*,20)sum 20 format(1x,'The value of the integral is: ',F10.4) STOP END function f(x) f=1/(1+x*x) return END
SIMPSON’S METHOD OD INTEGRATION (ARRAY GIVEN)
PROGRAM zz Simpson's method of integration dimension x(20),y(20) write(*,*)'To find the value of integral of a function using *Simpsons rule. ' write(*,*) write(*,*)'Enter the limits of integration:' read(*,*)a,b write(*,*)'Enter the number of subintervals(an even no.) :' read(*,*)n h=(b-a)/n x(1)=a y(1)=f(x(1)) DO 10 i=2,n+1 x(i)=x(i-1)+h y(i)=f(x(i)) 10 continue s=y(1)+y(n+1) DO 20 i=2,n,2 s=s+4*y(i) 20 continue DO 30 i=3,n-1,2 s=s+2*y(i) 30 continue s=h*s/3 write(*,40)s 40 format(1x,'The value of the integral is: ',F10.4) STOP END c
function f(x) f=sin(x) return END
TRAPEZOIDAL RULE OF INTEGRATION PROGRAM trp c Trapezoidal rule of integration integer count count=1 write(*,30) 30 format(1x,'To find the value of the integral of a function ',\) write(*,40) 40 format('using trapezoidal rule.') write(*,*) write(*,*)'Enter the limits of integration: ' read(*,*)a,b IF(b.LT.a)THEN c=a a=b b=c ENDIF write(*,*)'Enter the number of subintervals(an even no.) :' read(*,*)n h=(b-a)/n x=a y=f(x) sum=y 10 IF(count.LT.n)THEN x=x+h y=f(x) sum=sum+2*y count=count+1 GOTO 10 ENDIF x=x+h y=f(x) sum=sum+y sum=h*sum/2 write(*,*) write(*,20)sum 20 format(1x,'The value of the integral is :',F10.4) STOP END function f(x) f=1/(1+x*x) return END
TRAPEZOIDAL RULE OF INTEGRATION (ARRAY )
PROGRAM ttrr c Trapezoidal rule of integration dimension x(20),y(20) write(*,40) 40 format(1x,'To find the integral of a given function using ',\) write(*,50) 50 format('trapezoidal rule.') write(*,*)'Enter the limits of integration:' read(*,*)a,b write(*,*)'Enter the number of subintervals:' read(*,*)n IF(a.GT.b)THEN c=a a=b b=c ENDIF h=(b-a)/n x(1)=a y(1)=f(x(1)) write(*,*)x(1),y(1) DO 10 i=2,n+1 x(i)=x(i-1)+h y(i)=f(x(i)) write(*,*)x(i),y(i) 10 continue s=y(1)+y(n+1) DO 20 i=2,n s=s+2*y(i) 20 continue v=h*s/2 write(*,30)v 30 format(1x,'The value of the integral is: ',F10.4) STOP END function f(x) f=sin(x) return END
TRIANGULAR FACTORISATION TO SOLVE A LINEAR SYSTEM
PROGRAM tri c Triangular factorisation to solve a linear system of eqns real l(20,20) dimension a(20,20),b(20),u(20,20),x(20),y(20) write(*,130) 130 format(1x,'To solve a system of equations using Method of ') write(*,140) 140 format(\,'Triangularisation.') write(*,*)'Enter the number of equations:' read(*,*)n write(*,*)'Enter the coefficients of each of the equations:' read(*,*)((a(i,j),j=1,n),b(i),i=1,n) DO 10 i=1,n DO 20 j=i+1,n l(i,j)=0 20 continue DO 30 j=1,i-1 u(i,j)=0 30 continue u(1,i)=a(1,i) l(i,i)=1 IF(i.GT.1) l(i,1)=a(i,1)/u(1,1) 10 continue DO 40 i=2,n DO 50 j=2,i-1 s=0 DO 60 k=1,j-1 s=s+l(i,k)*u(k,j) 60 continue l(i,j)=(a(i,j)-s)/u(j,j) 50 continue DO 70 j=i,n s=0 DO 80 k=1,i-1 s=s+l(i,k)*u(k,j) 80 continue u(i,j)=a(i,j)-s 70 continue 40 continue
write(*,*)((l(i,j),i=1,n),j=1,n) write(*,*)((u(i,j),i=1,n),j=1,n) y(1)=b(1) DO 90 i=2,n sum=0 DO 100 j=1,n-1 sum=sum+l(i,j)*y(j) 100 continue y(i)=b(i)-sum 90 continue write(*,*)(y(i),i=1,n) x(n)=y(n)/u(n,n) DO 110 i=n-1,1,-1 sum=0 DO 120 j=i+1,n sum=sum+u(i,j)*x(j) 120 continue x(i)=(y(i)-sum)/u(i,i) 110 continue write(*,*)'The solution system is: ' DO 150 i=1,n write(*,160)x(i) 160 format(1x,F6.3) 150 continue STOP END
