Fortran Programs For Physics Students

FORTRAN PROGRAMS for Physics Students
(With THEORY and ALGORITHM)
Abhijit Kar Gupta
Dr. Abhijit Kar Gupta, Dept. of Physics, PBC, email: kg.abhi@gmail.com Page 1
Simple Programs
Area of a Triangle: Heron’s Formula
Theory:
For a triangle, if two sides and and the angle between them are given, the third side,
√ . The area of the triangle, √ , where
.
Algorithm:
1. Input: , and (in deg.)
2. Covert the angle into radian and find the third side from formula
3. Calculate ‘s’. Calculate area from formula. Write in a file.
FORTRAN Program:
C Area of a Triangle: Heron’s Formula

C
write(*,*)'Give Two sides and the adjacent angle'
read(*,*)a,b,theta
x = 3.14*theta/180 ! Converting into radian
c = sqrt(a**2+b**2-2*a*b*cos(x))
s =(a+b+c)/2.
area = sqrt(s*(s-a)*(s-b)*(s-c))
write(*,*)'Area =', area
stop
end
Prime Numbers:
Theory:
Prime numbers are those which are not divisible by any number except 1 and the same
number. For example, 3, 5, 7…. Here we check prime numbers between two numbers
Algorithm:
1. Do Loop starts from 2 up to n-1
2. Inside the loop, check if the number ‘n’ is divisible by any number
3. Set a counter to check how many times there are no divisors.
4. If the count = n-2, we find this is a prime number.
FORTRAN Program: C To find Prime Numbers
C
write(*,*)'Prime Numbers'
do 10 n=1,1000
k=0
do i=2,n-1
ii=mod(n,i)
if(ii.eq.0)go to 10
if(ii.ne.0)k=k+1
enddo
if(k.eq.n-2)write(*,*)n
10 continue
stop
end
Square Number:
Theory:
This is to check if a given number is a perfect square, for example, . So we take a
number and calculate its square root, the result will be an integer when perfect square, else a
real number. If the result is a real number (and not an integer), the integer part squared again,
will not match with the original number. In case of perfect square number, this matches.
Algorithm:
1. Input: the number
2. Take square root and keep the integer part.
3. Square the Integer part (in step 2).
4. Check if the squared number (in step 3) matches with the given number.
FORTRAN Program:
C Check if a number is a perfect square

C
write(*,*)'Give a number'
read(*,*)n
nsqrt = sqrt(real(n))
nsq = nsqrt*nsqrt
if(nsq.eq.n)then
write(*,*)'Perfect Square'
else
write(*,*)'NOT a Perfect Square'
endif
stop
end
TO CHECK LEAP YEAR:
Theory:
Earth takes 365 days, 5 hours, 48 minutes, and 45 seconds – to revolve once around the Sun.
But our Gregorian calendar considers 365 days in a year. So a leap day is regularly added (to the
shortest month February) in every 4 years to bring it in sync with the Solar year. After 100
years our calendar does have a round off by around 24 days.
In the Gregorian calendar three criteria must be taken into account to identify leap years:
 The year can be evenly divided by 4;
 If the year can be evenly divided by 100, it is NOT a leap year, unless;
 The year is also evenly divisible by 400. Then it is a leap year.
This means that in the Gregorian calendar, the years 2000 and 2400 are leap years,
while 1800,1900, 2100, 2200, 2300 and 2500 are NOT leap years.
Algorithm:
1. Input: the number
2. Check if it is not divisible by 100 but divisible by 4
3. Check if it is divisible by 400
4. For steps 2 and 3, the Number is Leap year, otherwise it is not.
C To Check Leap Year
integer y
write(*,*) 'enter year'
read(*,*) y
if(mod(y,100).ne.0.and.mod(y,4).eq.0)
then
write(*,*) 'leap year'
elseif(mod(y,400).eq.0) then
write(*,*) 'leap year'
else
write(*,*) 'not a leap year'
endif
end
Stirling’s Formula:
Theory:
Stirling’s formula to find out logarithm of factorial of a large number:
We shall verify this. But since usually would be a very large number for a big number, and so
the computer would not be able to handle it and store beyond some limit. Thus we use
∑ . [Here we avoid computing factorial!]
Algorithm:
1. Start a do-loop to sum the logarithm of the integers from 1 up to the given number.
2. Calculate the value from Stirling’s formula
3. Compare the results from (1) and (2)
FORTRAN Program:
C Stirling’s Formula
C
open(1,file='stir.d')
do n=10,10000,10
sum=0.0
do i=1,n
sum=sum+log(real(i))
enddo
app=n*log(real(n))-n
diff=sum-app
write(*,*)n,sum,app
write(1,*)n,diff/sum
enddo
stop
end
Maximum and Minimum from a list of Numbers:
Algorithm:
1. Call the first number to be maximum (or minimum).
2. Compare the next number and find the max (or min) between the two. Then redefine
the max (or min).
3. Repeat step 2
FORTRAN Program:
C To find Max and Min numbers from a list.

C
write(*,*)’How many numbers?’
read(*,*)n
write(*,*)’Give the numbers’
read(*,*)a
amax=a
amin=a
do i=2,n
read(*,*)a
if(a.ge.amax)amax=a
if(a.le.amin)amin=a
enddo
write(*,*)’Max = ‘, amax, ‘Min = ‘, amin
stop
end
Note: If you want to retain the numbers then you have to store them through ‘dimension’.
C To find Max, Min and the Range.

C
dimension x(100)
write(*,*)'How many numbers?'
read(*,*)n
write(*,*)’Give the numbers'
read(*,*)(x(i),i=1,n) !Implied do-loop
amax=x(1)
do i=1,n
if(x(i).gt.amax)amax=x(i)
enddo
C
amin=x(1)
do i=1,n
if(x(i).lt.amin)amin=x(i)
enddo
C
range=amax-amin
write(*,*)'Max, Min, Range'
write(*,*)amax, amin, range
stop
end
Sorting of Numbers from a List: (Ascending and Descending)
Theory:
The task is to read from a list of numbers and compare between first two numbers and decide
which is bigger (or smaller) and then compare with the next and so on. To do this we have to
store the one number into a temporary variable and then orient the two.
Algorithm:
1. Read the numbers from a list
2. Compare the first two numbers. If the 1st one is greater than the 2nd, store the 2nd
number into a temporary variable.
3. Now call the 2nd number equal to 1st number.
4. Assign the 1st number equal to the stored number (previous 2nd number).
5. Repeat 1-4 until the end of the list. We get the ascending order.
FORTRAN Program:
C Ascending Order
C
dimension a(100)
write(*,*)'How many numbers?'
read(*,*)n
write(*,*)'Give the numbers'
read(*,*)(a(i),i=1,n)
do i=1,n-1
if(a(i).gt.a(i+1))then
temp=a(i+1)
a(i+1)=a(i)
a(i)=temp
endif
enddo
write(*,*)(a(i),i=1,n)
stop
end
*For descending order, you just have to change ‘gt’ (greater than) in the if-statement to ‘lt’
(less than).
Digits of a decimal number:
Theory:
We have to read a number and find out its digits at different places. A decimal number is found
by multiplying by the powers of 10 with the digit in place and then summing. So the digits of a
decimal number can be found out by dividing the number by 10 and taking the remainder.
Next time, divide the number by 10 and take the integer part and repeat the steps. For
example, take the number 358, divide 358 by 10, the remainder is 8. Next time, consider
and take the integer part which is 35 and repeat the steps.
Algorithm:
1. Input: The Number, number of digits (n)
2. Start a do loop starting from 1 to n.
3. Find the remainder (Modulo 10).
4. Divide the number by 10 and retain the integer part for the next step.
5. Sum the digits as obtained.
FORTRAN Program:
C To find the digits and their sum

C
write(*,*)'Give the number'
read(*,*)n
nsum=0
write(*,*)’Digits in reverse order’
nsum=0
10 k=mod(n,10)
n=n/10.
nsum=nsum+k
write(*,*)k
if(n.gt.9)go to 10
write(*,*)n
write(*,*)'Sum of digits = ', nsum
stop
end
Factorial of a Number:
Theory: Factorial of a number, .
Algorithm:
1. Start a do loop for taking the product of numbers in increasing (or decreasing order)
2. Write the product which is the factorial of the given number
FORTRAN Program:
C Factorial of a Number
write(*,*)’Give the number’
read(*,*)n
nfac = 1
do i=2,n
nfac = nfac*i
enddo
write(*,*)’Factorial = ‘, nfac
stop
end
Strong Number:
Theory:
A strong number is something where the sum of factorials of its digits equals to the number
itself. For example: We here find out strong numbers between 1 and some
number.
Algorithm:
1. Start a do-loop of numbers from 1 to some big number
2. The number variable is renamed when used.
3. Find each digit of the number and find factorial of that digit.
4. Sum the factorial number
5. Compare the sum with the original number
[Note: Here we have defined a function subroutine to calculate factorial as that has to be
done repeatedly.]
FORTRAN Program:
C Strong Number
C
do k=1,100000
n=k
nsum=0
10 i=mod(n,10)
nsum=nsum+nfac(i)
n=n/10.
if(n.gt.9)go to 10
nsum=nsum+nfac(n)
if(nsum.eq.k)write(*,*)'Strong Number ', k
enddo
stop
end
C
C Subroutine
C
function nfac(i)
nfac=1
do k=2,i
nfac=nfac*k
enddo
return
end
*In the above program, we have defined a function in a ‘subroutine’ for our purpose. This is just
to make our program simpler to handle.
Armstrong Number:
Theory:
Armstrong number is a number such that the sum of its digits raised to the third power is equal
to the number itself. For example, 153 is an Armstrong number, since .
Algorithm:
1. Set a counter.
2. For a three digit number, three nested do loop start for three integers.
3. Construct a number from three integers.
4. Take cubes of the integers and sum them inside the innermost loop.
5. Check if the sum matches with the constructed number.
6. Count how many Armstrong numbers are there between 0 and 999.
FORTRAN Program:
C Armstrong Numbers between 0 and 999

C
n=0
do i=0,9
do j=0,9
do k=0,9
ijk=i*100+j*10+k
i3j3k3=i**3+j**3+k**3
if(ijk.eq.i3j3k3)then
write(*,*)ijk
n=n+1
endif
enddo
enddo
enddo
write(*,*)'Total Armstrong Numbers = ', n
stop
end
C Armstrong Numbers between two Numbers

C
Integer a, b
write(*,*)’Give a and b’
read(*,*)a,b
do n=a,b
k=n
nsum=0
10 i=mod(k,10)
k=k/10.
nsum=nsum+i**3
if(k.gt.9)go to 10
nsum=nsum+k**3
if(nsum.eq.n)’Armstrong Number ‘, n
enddo
stop
end
Palindrome Number:
Theory:
A Palindrome is some word or number or a phrase which when read from the opposite end,
appears the same (reflection symmetry, example, 34543, symmetry around 5 in the middle).
We here read a decimal number and find out its digits at different places. After finding the
digits we construct the number taking the digits in the reverse order. For example, given the
number 456, we construct the new number (this is not
palindrome number).
Algorithm:
1. Input: The Number, number of digits (n)
2. Store the Number.
3. Start a do loop starting from 1 to n.
4. Find the remainder (Modulo 10).
5. Divide the number by 10 and retain the integer part for the next step.
6. Do the sum after multiplying the digit with the appropriate power of 10 for that
place.
7. Check the sum with the stored number: If equal, Palindrome.
FORTRAN Program:
C Palindrome Number
C
write(*,*)'Give the number'
read(*,*)num
write(*,*)'Number of digits = ?'
read(*,*)n
num0=num ! To store the Number
nsum=0
do i=1,n
k=mod(num,10)
num=num/10.
nsum=nsum+nd*10**(n-i)
write(*,*)nd
endd
if(nsum.eq.num0)write(*,*)'Palindrome'
if(nsum.ne.num0)write(*,*)'Not Palindrome'
stop
end
Fibonacci Numbers:
Theory:
Fibonacci numbers are the series of numbers where each number is the sum of two previous
numbers: 0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89….
The rule:
Algorithm: (Numbers up to 610, for example)

1. Input: Two initial numbers 0 and 1.
2. Start a do-loop. Obtain the 3rd number as the sum of first two numbers.
3. Now redefine the 2nd number as the 1st number and the 3rd number as the 2nd
number.
4. Obtain the 3rd number as the sum of newly defined 1st and 2nd numbers and repeat
inside the loop. Write down the numbers on the screen as they are obtained.
5. Check if the new number in the series is 610, then exist the loop and write the sum.
FORTRAN Program:
C Fibonacci Numbers up to 610 (or any no. you like)

C
n1=0
n2=1
n3=n2+n1
write(*,*)’Fibonacci Numbers’
write(*,*)n1
write(*,*)n2
write(*,*)n3
do while(n3.lt.610)
n1=n2
n2=n3
n3=n2+n1
write(*,*)n3
enddo
write(*,*)'Golden Ratio= ', float(n3)/n2
stop
end
Golden Ratio = 1.6180…
GCD and LCM:
Theory:
This is to find Greatest Common Divisor (GCD) and Lowest Common Multiplier (LCM) of two
positive integers by Euclid method. First, we divide the bigger number by the smaller number
and consider the residue. If the residue is zero, then the smaller number is the GCD. If this is
not, then we take the nonzero residue as the smaller number and previous smaller number to
be the bigger number. This is continued until we get a zero residue. Then at the last step, the
divisor is called the GCD. LCM is just obtained by dividing the product of two original numbers
by GCD.
Algorithm:
1. Arrange the two numbers: bigger number = m, smaller number = n

2. Find residue by dividing “m” by “n” (through system function)
3. If the residue in 2 is zero, GCD = n
4. If the residue in 2 is not zero, redefine: m = n, n= residue
5. Repeat the steps 2 - 4.
6. Calculate LCM = m*n/GCD
C GCD & LCM

C
write(*,*)'Give m, n'
read(*,*)m,n
nprod=m*n
if(m.lt.n)then
k=m
m=n
n=k
endif
10 i=mod(m,n)
if(i.eq.0)then
write(*,*)'GCD = ', n
else
m=n
n=i
go to 10
endif
write(*,*)'LCM = ', nprod/n
stop
end
Sum of an Infinite Series
Exponential Series:
Theory:
Taylor series expansion:
x2 x3 x4
exp( x)  1  x     .............
2! 3! 4!
Algorithm:
1. Input: , Tolerance value

2. Initial values: Sum = 0, Term = 1.0
3. Do Loop starts (n = 1, 2…)
4. New Term = Previous Term
5. Sum the terms.
6. Check: Absolute value of a term when less than tolerance, come out.
7. Write the sum.
FORTRAN Program:
C Exponential Series with tolerance value

C
write(*,*)'Give x and tolerance value'
read(*,*)x,tol
sum=0.0
term=1
do n=1,100
if(abs(term).lt.tol)go to 10
sum=sum+term
term=term*x/n
end do
10 write(*,*)’Calculated value, Actual value’
write(*,*)sum,exp(x)
stop
end
Cosine Series:
Theory:
Taylor series expansion:
Algorithm:
1. Input: , tolerance value

2. Initial values: Sum = 0, Term
3. Do Loop starts (n = 1, 2….)
4. New Term = Previous term ( )
5. Sum the terms
7. Write the sum.
FORTRAN Program:
C Sum of Cosine Series with tolerance value

C
write(*,*)'Give x and tolerence value'
read(*,*)x,tol
sum=0.0
term=1.0
do n=1,100
sum=sum+term
term=term*(-x**2/(2*n*(2*n-1)))
enddo
10 write(*,*)’Calculated value, Actual Value’
write(*,*)sum,cos(x)
stop
end
Sine Series:
Theory:
Taylor series expansion of:
Algorithm:
1. Input: , tolerance value

2. Initial values: Sum = 0, Term
3. Do Loop starts (n = 1, 2….)
4. New Term = Previous term ( )
5. Sum the terms
7. Write the sum
FORTRAN Program:
C Sum of Sine Series with tolerance value

C
write(*,*)'Give x and tolerence value'
read(*,*)x,tol
sum=0.0
term=x
do n=1,100
sum=sum+term
term=term*(-x**2/(2*n*(2*n+1)))
enddo
10 write(*,*)’Calculated value, Actual Value’
write(*,*)sum,sin(x)
stop
end
Conversion of Numbers
From Decimal to Binary: (Integer)
Theory:
Base 10 system -> Base 2 system.
Continue downwards, dividing each new quotient by two and writing the remainders to the right of each
dividend. Stop when the quotient is 1.
For example, a decimal number, 156 (may be written as 15610 ) is converted to a binary number
10011100 (may be written as 100111002)
2)156 0 Algorithm: (For Integer Numbers)

2)78 0
2)39 1 1. Input: a number (integer)
2)19 1 2. Divide by 2 (Modulo)
2)9 1 3. Store the remainder in memory.
2)4 0 4. Take the Quotient for the next step.
2)2 0 5. Repeat steps 2 & 3 until the Quotient is 1 or 0.
1 6. When the Quotient is 1 0r 0, store this as remainder.
7. Write the stored remainders (either 0 or 1) in reverse order
FORTRAN Program:
C From Decimal to Binary

C
dimension ib(100)
write(*,*)'Give a Decimal Number'
read(*,*)n
do i=1,100
if(n.le.1)then
ib(i)= n
go to 10
endif
ib(i)= mod(n,2)
n=n/2
enddo
10 write(*,*)'The Binary Number’
write(*,*)(ib(k),k=i,1,-1) !Implied do-loop
stop
end
Decimal to Binary (For Real Numbers, fractional part considered):
Theory:
For the conversion of fractional part, we have to multiply by 2 and collect the real part of the
product.
Algorithm:
1. Input: a number (real number) 8. Now take the fractional part,
2. Take the integer part. Divide by 2 multiply by 2, store the integer part
3. Store the remainder in memory. of real number obtained.
4. Take the Quotient for the next step. 9. Subtract the above integer from the
5. Repeat steps 3 & 4 until the new real number and obtain a new
Quotient is 1 or 0. fraction.
6. When the Quotient is 1 or 0, store 10. Repeat the steps from 8-10 until the
this as remainder. fraction becomes 0 or stop after
7. Write the stored remainders (either some steps for recurring fraction.
0 or 1) in reverse order. 11. Write the integers
C From Decimal to Binary

C
dimension ib(100), ibc(100)
read(*,*)x
intx=x ! Integer part
fracx=x-intx ! Fractional part
do i=1,1000 ! Loop for Integer part
if(intx.le.1)then
ib(i)=intx
go to 10
endif
ib(i)= mod(intx,2)
intx=intx/2
enddo
10 do j=1,1000 ! Loop for Fractional part
fracx=fracx*2
ibc(j)=fracx
fracx=fracx-ibc(j)
if(fracx.eq.0.0.or.j.eq.8)go to 20
enddo
20 write(*,*)'Binary Number = ',
(ib(k),k=i,1,-1),'.',(ibc(k),k=1,j)
stop
end
From Decimal to Octal:
Theory:
Base 10 system -> Base 8 system.
We have to divide by 8 and we stop when the quotient is less than 8.
The octal (base 8) numeral system has 7 possible values, represented as 0, 1, 2, 3, 4, 5, 6 and 7
for each place-value.
Continue downwards, dividing each new quotient by 8 and store the remainders. Stop when
the quotient is less than 8. For example, a decimal number, 156 (may be written as 15610 ) is
converted to a Octal number 234 (may be written as 2348)
8)156 4 Algorithm: (For Integer Numbers)

8)19 3
2 1. Input: a number (integer)
2. Divide by 8 (Modulo)
3. Store the remainder in memory.
4. Take the Quotient for the next step.
5. Repeat steps 2 & 3 until the Quotient is less than 8.
6. When the Quotient less than 8, store this as remainder.
7. Write the stored remainders (between 0 & 7) in reverse
order
FORTRAN Program:
C From Decimal to Octal

C
dimension ib(100)
read(*,*)n
do i=1,100
if(n.le.7)then
ib(i)=n
go to 10
endif
ib(i)= mod(n,8)
n=n/8
enddo
10 write(*,*)'The Octal Number’
write(*,*)(ib(k),k=i,1,-1) ! Implied do-loop
stop
end
Decimal to Octal (For Real Numbers, fractional part considered):
Theory:
For the conversion of fractional part, we have to multiply by 2 and collect the real part of the
product.
Algorithm:
1. Input: a number (real number) 8. Now take the fractional part,
2. Take the integer part. Divide by 8 multiply by 8, store the integer part
3. Store the remainder in memory. real number obtained.
4. Take the Quotient for the next step. 9. Subtract the above integer from the
5. Repeat steps 3 & 4 until the new real number and obtain a new
Quotient is less than 8. fraction.
6. When the Quotient is less than 8, 10. Repeat the steps from 8-10 until the
store this as remainder. fraction becomes 0 or stop after
7. Write the stored remainders some steps if recurring.
(between 0 and 7) in reverse order 11. Write the stored integers
C From Decimal to Octal

C
Dimension io(100), ioc(100)
read(*,*)x
intx=x ! Integer part
fracx=x-intx ! Fractional part
do i=1,1000 ! Loop for integer part
if(intx.le.7)then
io(i)=intx
go to 30
endif
io(i)=mod(intx,8)
intx=intx/8
enddo
30 do j=1,1000 ! Loop for fractional part
fracx=fracx*8
ioc(j)=fracx
fracx=fracx-ioc(j)
if(fracx.eq.0.0.or.j.ge.8)go to 40
enddo
40 write(*,*)'Octal Number = ',
(io(k),k=i,1,-1),'.',(ioc(k),k=1,j)
stop
end
From Binary to Decimal:
Theory:
To convert a binary number: List the powers of 2 from right to left. Start at 2 0, next increase the
exponent by 1 for each power. Stop when the amount of elements in the list is equal to the
amount of digits in the binary number. For example, the number 10011 has 5 digits, so the
corresponding decimal number = = 19
Algorithm:
1. Input: Binary number, total number of digits
2. Define the binary number as a ‘character’ so as to read each digit in it.
3. Find each digit of the number (either 0 or 1) and its position (n) counted from right to
left.
4. The digit (0 or 1) found from step 3 is multiplied by
5. Obtain the sum in a do-loop.
6. Write the sum which is the decimal number.
FORTRAN Program:
C
C From Binary to Decimal
C
write(*,*)'Give a Binary Number'
read(*,*)number
write(*,*)'Total Number of Digits?'
read(*,*)n
nsum=0
do i=1,n
nd=mod(number,10)
number=number/10.
nsum=nsum+nd*2**(i-1)
enddo
write(*,*)'Decimal Number = ',nsum
stop
end
From Octal to Decimal:
Theory:
To convert a octal number: List the powers of 8 from right to left. Start at 8 0, next increase the
exponent by 1 for each power. Stop when the amount of elements in the list is equal to the
amount of digits in the binary number. For example, the number 234 has 3 digits, so the
corresponding decimal number = = 156
Algorithm:
1. Input: Octal number, total number of digits
2. Define the Octal number as a ‘character’ so as to read each digit in it.
3. Find each digit of the number (0, 1, 2, 3, 4, 5, 6 or 7) and its position (n) counted
from right to left.
4. The digit (0 to 7) found from step 3 is multiplied by
5. Obtain the sum in a do-loop.
6. Write the sum which is the decimal number.
FORTRAN Program:
C
C From Octal to Decimal
C
write(*,*)'Give an Octal Number'
read(*,*)number
write(*,*)'Total Number of Digits?'
read(*,*)n
nsum=0
do i=1,n
nd=mod(number,10)
number=number/10.
nsum=nsum+nd*8**(i-1)
enddo
write(*,*)'Decimal Number = ',nsum
stop
end
Matrices
Addition of Two Matrices:
Theory:
For the addition of two matrices, the rows and columns of the two matrices have to be equal.
The th element of a matrix is added with the same th element of another to obtain
the th element of the new matrix: .
Algorithm:
1. Input and output Matrix elements are stored in memory.
2. Input: Dimension of Matrices, Matrix elements
3. Do loop for two indices, i and j, loops are nested.
4. The -th element of a matrix is added with the same th element of another to
obtain the th element of the new matrix: .
5. Write down the matrix elements with implied do loop.
FORTRAN Program:
C Addition of Two Matrices

C
dimension a(10,10),b(10,10),c(10,10)
write(*,*)'Rows & Columns of Matrices'
read(*,*)m,n
write(*,*)'Matrix elements of A?'
do i=1,m
read(*,*)(a(i,j),j=1,n)
enddo
write(*,*)'Matrix elements of B?'
do i=1,m
read(*,*)(b(i,j),j=1,n)
enddo
do i=1,m
do j=1,n
c(i,j)= a(i,j)+b(i,j)
enddo
enddo
Write(*,*)'The added Matrix'
do i=1,m
write(*,*)(c(i,j),j=1,n)
enddo
stop
end
Product of Two Matrices:
Theory:
∑ , for the Product of two matrices, the number of columns of the
first matrix has to be the same as the number of rows of the 2nd matrix.
Algorithm:
1. Input and output Matrix elements are stored in memory.
2. Input: Dimensions of Matrices, Matrix elements
3. Do loop for three indices: i, j, k, loops are nested
4. Inside the k-loop, the th element of 1st matrix is multiplied with the th
nd
element of 2 matrix, and then sum is taken.
5. Write down the matrix elements with implied do loop.
FORTRAN Program:
C Product of Two Matrices
C
dimension a(10,10),b(10,10),c(10,10)
write(*,*)'Rows & Columns of Matrix A'
read(*,*)m,l
write(*,*)'Matrix elements of A?'
do i=1,m
read(*,*)(a(i,j),j=1,l)
enddo
write(*,*)'Rows & Columns of Matrix B'
read(*,*)l,n
write(*,*)'Matrix elements of B?'
do i=1,l
read(*,*)(b(i,j),j=1,n)
enddo
do i=1,m
do j=1,n
c(i,j)= 0.0
do k=1,l
c(i,j)= c(i,j)+a(i,k)*b(k,j)
enddo
enddo
enddo
Write(*,*)'The Product Matrix'
Do i=1,m
write(*,*)(c(i,j),j=1,n)
enddo
stop
end
Transpose of a Matrix:
Theory:
For the transpose of a matrix, the rows of the original matrix become the columns of the
transposed matrix and vice versa. So what we do is that we read the matrix elements and find
the elements of the transposed matrix:
Algorithm:
1. Matrix elements are stored in memory
2. Input: rows and columns, matrix elements
3. Do loop for i and j:
4. Transpose matrix elements written in implied do loop
FORTRAN Program:
C Transpose of a Matrix
C
dimension a(10,10),b(10,10)
write(*,*)'No. of rows and columns?'
read(*,*)m,n
write(*,*)'Matrix elements?'
do i=1,m
read(*,*)(a(i,j),j=1,n)
enddo
C
do i=1,n
do j=1,m
b(i,j)=a(j,i)
enddo
enddo
C
Write(*,*)'The Transposed Matrix'
Do i=1,n
write(*,*)(b(i,j),j=1,m)
enddo
stop
end
Numerical Analysis
Integration by Simpson’s 1/3rd Rule:
Theory:
h ab 
b
 f ( x)dx  3  f (a)  4 f 
a
  f (b)
2  
[See APPENDIX I for the detailed theory and calculations.]
Algorithm:
1. The function is defined in a subroutine.

2. Input: Lower limit , Upper limit
3. Calculate :
4. Evaluate: , , )
5. Calculate the sum: ( )
6. Multiply step 5 by
FORTRAN Program:
C Simpson’s 1/3 Rule
C
write(*,*)'Lower limit, Upper Limit’
read(*,*)a,b
h=(b-a)*0.5
x=(a+b)*0.5
sum=f(a)+4*f(x)+f(b)
s=h/3*sum
write(*,*)'Value of Integral= ', s
stop
end
C
C Subroutine for function
function f(x)
f=x**2
return
end
Integration by Simpson’s 1/3rd Rule (Modified):
Theory: ∫

h
 f ( x0 )  4 f ( x1 )  f ( x3 )  ...  f ( xn1 )  2 f ( x2 )  f ( x4 )  ...  f ( xn2 )  f ( xn )
3
Algorithm:
1. The function is defined in a subroutine.

2. Input: Lower limit , Upper limit , sub-intervals .
3. Calculate:
4. Evaluate:
5. Do Loop over n starts. Evaluate the function: for every odd sub-interval multiply by 4,
for every even sub-interval multiply by 2 and then add all of them.
6. Finally, multiply the sum by
FORTRAN Program:
C Simpson’s 1/3 Rule (Modified)

C
write(*,*)'Lower limit, Upper Limit,sub-intervals'
read(*,*)a,b,n
h=(b-a)/n
sum=f(a)+f(b)
d=4
do k=1,n-1
x=a+k*h
sum=sum+d*f(x)
d=6-d
enddo
s=h/3*sum
write(*,*)'Value of Integral= ', s
stop
end
C
C Subroutine for function
C
function f(x)
f=x**2
return
end
Least Square Method: Linear Fit
∑ ∑ ∑ ∑ ∑ ∑ ∑
Theory: Slope = , Intercept =
∑ ∑ ∑ ∑
[See APPENDIX I for the detailed theory and calculations.]
Algorithm:
1. Input: x, y 4. Few data points generated with

2. Inside the Do Loop, run the sums to obtained slope and intercept, the
evaluate: ∑ , ∑ , ∑ , ∑ straight line drawn is superposed
3. Calculation of slope and intercept over the original data points.
C Least Square fit

C
open(1,file='xy.dat')
open(2,file='fit.dat')
write(*,*)'Number of Points?'
read(*,*)n
sumx=0.0
sumy=0.0
sumsqx=0.0
sumxy=0.0
write(*,*)'Give data in the form: x,y'
do i=1,n
read(*,*)x,y
write(1,*)x,y
sumx=sumx+x
sumy=sumy+y
sumsqx=sumsqx+x*x
sumxy=sumxy+x*y
enddo
deno=n*sumsqx-sumx*sumx
slope=(n*sumxy-sumx*sumy)/deno
b=(sumsqx*sumy-sumx*sumxy)/deno
write(*,*)'Slope, Intercept= ',slope,b
C
write(*,*)'Give lower and upper limits of X'
read(*,*)xmin, xmax
x=xmin
dx=(xmax-xmin)/2.0
do i=1,3
y=slope*x+b
write(2,*)x,y
x=x+dx
enddo
stop
end
Note: The following programs are based on solving simple first order differential equations,
basically by Euler’s scheme. Runge-Kutta method of different orders can be used when necessary.
But we avoided that here for this set of elementary programs.
Projectile Motion
Theory:
Consider a particle projected at an angle from origin (

), with the horizontal direction being
x-axis and the vertical direction y-axis. The force on the particle, ⃗ ⃗. Components:
. For the purpose of computation, we take ; Mass,
. Projection angle, and initial speed are given. All the quantities are in arbitrary
units.
Algorithm:
1. Input: initial velocity , angle of projection . The interval of time is set.

2. Initial values: interval of time , Force components
3. Components of initial velocity: ,
4. Change in velocity components: ,
5. New velocity components: ,
6. Change in position coordinates: ,
7. New position coordinates: ,
8. Write the position coordinates in a file and then plot to see.
The trajectory of the particle projected at an angle (indicated in the figs.) and with some initial
velocity
FORTRAN Program:
C Projectile Motion of a particle

C
open(1,file='proj.dat')
write(*,*)'Angle (in degree) and initial speed?'
read(*,*)theta,v0
dt=0.01
g=980
theta=3.14/180*theta
fx=0
fy=-g
vx=v0*cos(theta)
vy=v0*sin(theta)
x=0
y=0
do while(y.gt.0.0)
dvx=fx*dt
dvy=fy*dt
vx=vx+dvx
vy=vy+dvy
dx=vx*dt
dy=vy*dt
x=x+dx
y=y+dy
write(1,*)x,y
enddo
stop
end
Motion under Central Force
Theory:
Consider the Central Force, and Newton’s 2nd law of motion: . For
computational purpose, we take , . Thus we have, . We solve for and

then therefrom derive the position. Initial values for the velocity components, and are
provided. We choose the time interval ( appropriately. All the quantities are in arbitrary
units.
Algorithm:
1. Initial position coordinates are given. Initial velocity components and are
given.
2. Time interval is chosen.
3. Components of Force: ,
4. √
5. Change in velocity components: ,
6. New velocity components: ,
7. New position coordinates: ,
8. Write the position coordinates in a file and then plot to see.
FORTRAN Program:
C Motion under a Central Force

C
open(1,file='central.dat')
dt=0.01
vx=0.0
vy=1.0
x=1.0
y=0.0
ncount=0
do i=1,10000
r2=x*x+y*y
r=sqrt(r2)
f=-1/r2
fx=x/r*f
fy=y/r*f
vx=vx+fx*dt
vy=vy+fy*dt
x=x+dt*vx
y=y+dt*vy
n=mod(i,100)
if(n.eq.0)write(1,*)x,y
enddo
stop
end
The Elliptical path of a particle under Central Force
Harmonic Oscillator
Theory:
Consider the Oscillator equation in the form: . For computational purpose, we

choose mass and the force constant a positive number (all in arbitrary units). The
equation is solved to obtain the velocity and therefrom we obtain the position of the particle.
Algorithm:
1. Give the initial position ( ), initial time ( ) and velocity ( ). Set the time interval .
2. Value of spring constant is chosen appropriately, usually a small positive number.
3.
4. Updating of velocity:
5. Change in position:
6. Updating of position:
7. Updating of time:
8. Write the time and position in a file and then plot to see.
FORTRAN Program:
C Harmonic Oscillator
C
open(1,file='harmonic.dat')
x=0.0
t=0
v=1.0
dt=0.01
k=5
do i=1,1000
dv=-k*x*dt
v=v+dv
dx=v*dt
x=x+dx
t=t+dt
write(1,*)t,x,v
enddo
stop
end
The Sine Curve obtained from the solution of Linear Harmonic Oscillator Equation.
**Equation for SHM [ ] is solved numerically in the above with the simple Euler scheme.
One can add a damping term (proportional to velocity) and tune the parameters to obtain well known
results and phase diagrams. However, it would be interesting to play with the parameters and gain
some insight through such numeric. See APPENDIX II.
Fourth Order Runge-Kutta (RK4) Method:
For the first order differential equation, with the initial value at , the
following steps are performed in order to find the value at the next grid point.
For the grid point , we compute
Then the value of is given by .
RK4 for Van der Pole Oscillator:
̈ ̇
Fortran Program:
open(1,file=’shmrk4.dat')
C
t1=0
x1=0
y1=0.1
h=0.01
C
write(1,*)t1,x1,y2
do i=1,10000
a1=f(t1,x1,y1)
b1=g(t1,x1,y1)
a2=f(t1+0.5*h,x1+0.5*h*a1,y1+0.5*h*b1)
b2=g(t1+0.5*h,x1+0.5*h*a1,y1+0.5*h*b1)
a3=f(t1+0.5*h,x1+0.5*h*a2+y1+0.5*h*b2)
b3=g(t1+0.5*h,x1+0.5*h*a2,y1+0.5*h*b2)
a4=f(t1+h,x1+h*a3+y1+h*b3)
b4=g(t1+h,x1+h*a3,y1+h*b3)
C
x2=x1+h/6*(a1+2*a2+2*a3+a4)
y2=y1+h/6*(b1+2*b2+2*b3+b4)
t2=t1+h
write(1,*)x2,y2,t2
x1=x2
y1=y2
t1=t2
enddo
c
stop
end
C
function f(t,x,y)
f=y
return
end
C
function g(t,x,y)
eps=10
g=-x-eps*(x**2-1)*y
return
end
RK4 for Chaos Theory: Lorenz Equations:
The Values of the parameters, Lorez used, and
We can solve this set of coupled 1st order equations by RK4 method.
FORTRAN Prog for RK4 Method:
C LORENZ EQUATIONS by RK4

implicit real*8(a-h,o-z)
open(1,file='rk4chaos.dat')
C
C Initial Values
C
t1=0
x1=0
y1=1.0
z1=0
C
h=0.01
h2=h/2
h6=h/6
C
write(1,*)y2,z2,x2,t2
C
do i=1,10000
a1=f1(x1,y1,z1)
b1=f2(x1,y1,z1)
c2=f3(x1,y1,z1)
C
x11=x1+h2*a1
y11=y1+h2*b1
z11=z1+h2*c1
C
a2=f1(x11,y11,z11)
b2=f2(x11,y11,z11)
c2=f3(x11,y11,z11)
C
x12=x1+h2*a2
y12=y1+h2*b2
z12=z1+h2*c2
C
a3=f1(x12,y12,z12)
b3=f2(x12,y12,z12)
c3=f3(x12,y12,z12)
C
x12=x1+h*a3
y12=y1+h*b3
z12=z1+h*c3
C
a4=f1(x13,y13,z13)
b4=f2(x13,y13,z13)
c4=f3(x13,y13,z13)
C
x2=x1+h6*(a1+2*2+2*a3+a4)
y2=y1+h6*(b1+2*b2+2*b3+b4)
z2=z1+h6*(c1+2*c2+2*c3+c4)
t2=t1+h
C
write(1,*)x2,y2,z2,t2
C
C Update…
C
x1=x2
y1=y2
z1=z2
t1=t2
enddo
C
stop
end
C
function f1(x,y,z)
sigma=10
f1=sigma*(y-x)
return
end
C
function f2(x,y,z)
rho=28
f2=x*(rho-z)-y
return
end
C
function f3(x,y,z)
beta=8/3.0
f3=x*y-beta*z
return
end
* Plotting x against y and z against z, we get the attractors.
Real Roots by Bisection Method
Theory:
Consider an equation, . If a real root exists between[ ], we search that root by
bisecting the interval; the midpoint being . If , is the root. If not
so, we choose either [ ] or [ ]. We check whether has opposite signs at the end
points (over any of the intervals). Then we take that new interval and bisect again. Thus we
keep searching the root according to some accuracy (we call it tolerance).
Algorithm:
1. Function defined in subroutine

2. Input:
1. If , No Root. Stop
2. Else set
3. If , is the root. Stop
4. If , set . Else set: .
5. Check if | | . Root = . Stop
6. Else go to step 3
Example:
Consider the equation: . Let us find a real root in the interval, [ ].
The following is the result of the iterations that is obtained by running the program.
FORTRAN Program:
C Bolzano Bisection Method
write(*,*)'Lower Limit, Upper Limit, Tolerance'

read(*,*)a,b,eps
C
C Testing if the root exists in this range
C
if(f(a)*f(b).gt.0.0)then
write(*,*)'No Root'
stop
endif
C
C Bisection Loop starts
C
do while(abs(b-a).ge.eps)
xm=0.5*(a+b)
if(f(xm).eq.0.0)go to 10
if(f(a)*f(xm).lt.0.0)then
b=xm
else
a=xm
endif
enddo
10 write(*,*)'Root = ',xm
stop
end
C
function f(x)
f=x**3-3*x-5
return
end
No. No.
Iterations Iterations
1 11.000000 12 2.2768555
2 6.5000000 13 2.2790527
3 4.2500000 14 2.2779541
4 3.1250000 15 2.2785034
5 2.5625000 16 2.2787781
6 2.2812500 17 2.2789154
7 2.1406250 18 2.2789841
8 2.2109375 19 2.2790184
9 2.2460938 20 2.2790356
10 2.2636719 21 2.2790270
11 2.2724609
Roots by Newton-Raphson Method:
Theory:
If we can find the derivative such that and also know that is monotonic in
[ ], we have:
Starting from an approximate value , we can generate the sequence
Algorithm:
1. Define and derivative in subroutine

2. Input: Guess value of , Tolerance, Max. no. of Iterations
3. Iteration in do loop:
4. Update:
5. Check if | | , then either or is the root. Stop.
6. Else continue until max iterations.
FORTRAN Program:
C Newton-Raphson Method
write(*,*)'Initial approx, Tolerance, Max iteration'

read(*,*)x,tol,maxit
do i=1,maxit
x0=x
x=x-f(x)/f1(x)
if(abs(x-x0).lt.tol)then
root=x
write(*,*)'Root =', root
stop
endif
if(i.eq.maxit)write(*,*)'No Solution'
enddo
stop
end
C
function f(x)
f=x**3-3*x-5
return
end
C
function f1(x)
f1=3*x**2-3
return
end
APPENDIX I
Least Square Fit:
Theory:
The following is a theory for data to be fitted with a straight line.
The equation of a straight line is y  mx  b

Consider the data points ( x1 , y1 ) , ( x2 , y 2 )…….etc.
n
Error is defined as  (m, b)   (mxi  b  yi ) 2 .
i 1
For the best fit, this error should be minimum.
 
Therefore,  0 and  0.
m b
  n 2
n

Now,  
m m  i 1
 ( mxi  b  y i )  = 
 i 1 m
(mxi  b  yi ) 2
n

=  2.(mx
i 1
i  b  yi )
m
(mxi  b  yi )
n n n n
= 2 (mxi  b  yi ) xi = 2m xi  2b xi   yi xi = 0 (1)
2
i 1 i 1 i 1 i 1
 n n
Similarly,  0  m xi  nb   yi (2)
b i 1 i 1
From (1) and (2),
 n
  n 
 n

 i 1
xi    yi 
    i 1
b 
 n n    n
2  m  
  xi  xi    y i xi 
 i 1 i 1   i 1 
n n n n n n n
n xi y i   y i  xi  y i  xi   xi  xi y i
2
Slope, m  i 1
n
i 1
n
i 1
and Intercept, b  i 1 i 1
n
i 1
n
i 1
.
n  xi  ( xi ) n xi  ( xi )
2 2 2 2
i 1 i 1 i 1 i 1
Example:
For the data points (1,2), (2,3), (3,4), (4,5)
n n
n  4,  xi  1  2  3  4  10 ,
i 1
y
i 1
i  2  3  4  5  14
n n
x y  1  2  2  3  3  4  4  5  40 , x  1  1  2  2  3  3  4  4  30
2
i i i
i 1 i 1
4  40  14  10 160  140 20 14  30  10  40 420  400 20
m    1, b     1.
4  30  10 2
120  100 20 4  30  10 2 120  100 20
Fig. Example of a linear least square fit
Simpson’s 1/3 Rule:
Theory:
Simpson’s 1/3 rule is where the integrand is approximated by a second order polynomial.
Fig. Integration of a function
b b
I   f ( x)dx   f 2 ( x)dx
a a
where f 2 ( x) is a second order polynomial given by

f 2 ( x)  a0  a1 x  a2 x 2
.
Choose
 a  b  a  b 
(a, f (a)),  , f  , and (b, f (b))
 2  2 
as the three points of the function to evaluate a0 , a1 and a 2 .
f (a)  f 2 (a)  a0  a1a  a2 a 2
ab ab ab ab
2
f   f2    a0  a1    a2  
 2   2   2   2 
f (b)  f 2 (b)  a0  a1b  a2 b 2
Solving the above three equations for unknowns, a0 , a1 and a2 give
ab
a 2 f (b)  abf (b)  4abf    abf (a)  b f (a)
2
a0   2 
a  2ab  b 2
2
ab ab
af (a)  4af    3af (b)  3bf (a)  4bf    bf (b)
a1    2   2 
a 2  2ab  b 2
 ab 
2 f (a)  2 f    f (b) 
 2 
a2   
a  2ab  b
2 2
Then
b
I   f 2 ( x)dx
a
 
b
  a0  a1 x  a 2 x 2 dx
a
b
 x2 x3 
 a0 x  a1  a2 
 2 3 a
b2  a2 b3  a 3
 a0 (b  a)  a1  a2
2 3
Substituting values of a 0 , a1 and a 2 gives
ba  ab 
b
f
a
2 ( x)dx 
6  f (a)  4 f    f (b)
 2  
Since for Simpson 1/3 rule, the interval a, b is broken into 2 segments, the segment width
ba
h
2
Hence the Simpson’s 1/3 rule is given by
h ab 
b
a f ( x)dx  3  f (a)  4 f  2   f (b)
Since the above form has 1/3 in its formula, it is called Simpson’s 1/3 rule.
Simpson’s 1/3 Rule (modified version)

(Multiple-segment Simpson’s 1/3 Rule)
Theory:
For better result, one can subdivide the interval a, b into n segments and apply Simpson’s
1/3 rule repeatedly over every two segments. Note that n needs to be even. Divide
interval a, b into n equal segments, so that the segment width is given by
ba
h .
n
Now
b xn
 f ( x)dx   f ( x)dx
a x0
where
x0  a
xn  b
b x2 x4 xn  2 xn
 f ( x)dx   f ( x)dx   f ( x)dx  ......   f ( x)dx   f ( x)dx

a x0 x2 xn  4 xn  2
Apply Simpson’s 1/3rd Rule over each interval,

 f ( x )  4 f ( x1 )  f ( x2 )   f ( x 2 )  4 f ( x3 )  f ( x 4 ) 
b
a f ( x)dx  ( x2  x0 ) 0 6 

 ( x4  x2 ) 
 6   ...
 f ( x n  4 )  4 f ( x n 3 )  f ( x n  2 )   f ( xn2 )  4 f ( xn1 )  f ( xn ) 
 ( xn2  xn4 )    ( xn  xn2 )  
 6   6 
Since
xi  xi 2  2h
i  2, 4, ..., n
then
 f ( x0 )  4 f ( x1 )  f ( x2 )   f ( x 2 )  4 f ( x3 )  f ( x 4 ) 
b
 f ( x)dx  2h 
a
6 

 2h 
 6   ...
 f ( x n  4 )  4 f ( x n 3 )  f ( x n  2 )   f ( xn2 )  4 f ( xn1 )  f ( xn ) 
 2h    2h  
 6   6 

h
 f ( x0 )  4 f ( x1 )  f ( x3 )  ...  f ( xn1 )  2 f ( x2 )  f ( x4 )  ...  f ( xn2 )  f ( xn )
3
APPENDIX II
Oscillations: a pictorial tour

(Results obtained from solving corresponding differential equations in computer)
Simple Harmonic Motion

(No damping)
𝑑 𝑋
𝑘𝑋
𝑑𝑡
Damped Harmonic Motion
(Critical Damping)
𝑑 𝑋 𝑑𝑋
𝜆 𝑘 𝑋
𝑑𝑡 𝑑𝑡

(Large damping)
(Small damping)
Relaxation Oscillation
Van der Pol Oscillator:

𝑑 𝑋 𝑑𝑋
𝜆 𝑋 𝑘𝑋
𝑑𝑡 𝑑𝑡
Exercises:
1. Consider the following matrix, ( ). Write a computer program to obtain the

corresponding transpose matrix . Also find the product . Display the result.
2. Write a computer program to obtain the first 100 numbers in the Fibonacci sequence and find
the golden ratio. [Golden ratio is the ratio of successive terms.] Is the Golden ratio tending
towards a fixed and finite number?
3. Write a program using the Newton-Raphson method to determine the roots of the equation:
.
4. Write a program to determine the roots of the equation:
that are close to .
5. Write a program to calculate the correlation coefficient for the 5 pairs of experimental
observations corresponding to some measurement: (2,5), (4,9), (5,11), (6,10), (8,12). Use the
∑ ∑ ∑
following formula:
√ ∑ ∑ √ ∑ ∑
6. Write a computer program to convert a binary number 10101 into its decimal
equivalent and then add the digits of the decimal number.
7. Compute numerically the following integral to verify the answer or obtain a value close
to that: ∫
8. You are given a set of height-weight data as follows: (170, 65), (172, 66), (181, 69), (157,
55), (150, 51). Write a program to find correlation coefficient, , where
∑ ̅ , ∑ ̅ and ∑ ̅ ̅ .
9. Write a program to convert the binary number 1000,000,000 into a decimal number.
10. = . Is that correct? Verify this through a computer program.
11. Write a program to find the transpose of the following matrix: ( ). Also
verify, .
12. Using Simpson’s modified 1/3 rd rule, write a program to calculate: ∫ .
13. Consider the number: 787. Is the number a prime? Test this through computer.
14. Compute this: ∫ , where we have
1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4

6.050 7.389 9.025 11.023 13.464 16.445 20.086 24.533 29.964
15. Write a program in computer to find out the through infinite series. Calculate up
to an accuracy level of . Find the value of .
16. Find the value of from the infinite series: up to 6th decimal
of accuracy.
17. Use Newton-Raphson method to find a non-zero solution of .
18. Write a computer program for a matrix A where the elements are
and another matrix B whose elements are and then compute the
matrix, .
19. For a cart sliding down an inclined plane, the relation between the distance and
time is given by, . A set of measurements are done as following.
(sec) 0.7 1.3 1.9 2.5 3.1 3.7 4.1 4.9 5.6 6.1 6.7 7.5 7.9 8.5 9.1
(cm) 1 4 9 16 25 36 49 64 81 100 121 144 169 196 225
Write a program to fit this data set. If you know the angle of inclination, , how do you
find the acceleration of gravity, from this measurements and your fitting exercise?
20. Given a number, 8956, write a program to compute the sum of the digits of this
number.
21. Given two matrices, ( ) and ( ), write a computer program to
compute .
compute, .
23. Write a program to verify approximately, . Calculate also
the percentage of error you incurred.
24. Write a suitable program to approximately verify, ∫ √ . Sate the accuracy
level that you are dealt with.
confirm that .
26. Write a program where the function is
defined in a subroutine and then obtain .
27. In a Nuclear detector, the counting rates vs. time are measured in the following table. If
you would draw a mean straight line through them, what would be the slope of that?
Time 0 15 30 45 60 75 90 105 120 135

(sec)
Counts 106 80 98 75 74 73 49 38 37 22
per sec
28. Convert to equivalent binary.
29. A set of 20 numbers are given: 1, 0.1, 5, 3, 10, -1, 4, 20, 1000, -9, 2, 14, 4.5, 0.9, 30, 9.8,
11, 22, 48, -10. Write a program to count how many numbers are between 0 to 10 and
how many are between 50 to 1000.
( )
30. While solving a cubic equation, we found a recursion relation: .
We are told that a root exists between 1 and 2. Starting from an initial guess, you have
to find out the root. Write a computer program for that.
31. Solve to 5 decimal places. You can start with an initial guess .
32. Verify this approximately up to an accuracy of order by
writing a computer program.
33. Write a program to show: .
This set of programs is intended to be distributed among physics students (U.G. and P.G.) for their
foundational coursework in Computation.
NOTE: In this basic course, I have avoided using formatted input and output. Format statements
can be used in the ‘read’ and ‘write’ statements…..
FORTRAN Programs for Physics Students

(With Theory and Algorithms)
Abhijit Kar Gupta
Department of Physics, Panskura Banamali College

Panskura R.S., East Midnapore, West Bengal, India, Pin Code: 721152
e-mail: kg.abhi@gmail.com

Fortran Programs For Physics Students

Uploaded by

Copyright:

Available Formats

Fortran Programs For Physics Students

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Fortran Programs For Physics Students

Uploaded by

Copyright:

Available Formats

FORTRAN PROGRAMS for Physics Students

(With THEORY and ALGORITHM)

Abhijit Kar Gupta

Area of a Triangle: Heron’s Formula

C Area of a Triangle: Heron’s Formula

C Check if a number is a perfect square

C To Check Leap Year

Maximum and Minimum from a list of Numbers:

C To find Max and Min numbers from a list.

C To find Max, Min and the Range.

C To find the digits and their sum

Theory: Factorial of a number, .

C Armstrong Numbers between 0 and 999

C Armstrong Numbers between two Numbers

Algorithm: (Numbers up to 610, for example)

C Fibonacci Numbers up to 610 (or any no. you like)

Golden Ratio = 1.6180…

1. Arrange the two numbers: bigger number = m, smaller number = n

C GCD & LCM

1. Input: , Tolerance value

C Exponential Series with tolerance value

Taylor series expansion:

1. Input: , tolerance value

C Sum of Cosine Series with tolerance value

Taylor series expansion of:

1. Input: , tolerance value

C Sum of Sine Series with tolerance value

From Decimal to Binary: (Integer)

2)156 0 Algorithm: (For Integer Numbers)

C From Decimal to Binary

C From Decimal to Binary

8)156 4 Algorithm: (For Integer Numbers)

C From Decimal to Octal

C From Decimal to Octal

C Addition of Two Matrices

Integration by Simpson’s 1/3rd Rule:

1. The function is defined in a subroutine.

Integration by Simpson’s 1/3rd Rule (Modified):

1. The function is defined in a subroutine.

C Simpson’s 1/3 Rule (Modified)

Least Square Method: Linear Fit

1. Input: x, y 4. Few data points generated with

C Least Square fit

Consider a particle projected at an angle from origin (

1. Input: initial velocity , angle of projection . The interval of time is set.

C Projectile Motion of a particle

Motion under Central Force

computational purpose, we take , . Thus we have, . We solve for and

C Motion under a Central Force

Consider the Oscillator equation in the form: . For computational purpose, we

Then the value of is given by .

RK4 for Van der Pole Oscillator:

RK4 for Chaos Theory: Lorenz Equations:

The Values of the parameters, Lorez used, and

FORTRAN Prog for RK4 Method:

C LORENZ EQUATIONS by RK4

Real Roots by Bisection Method

1. Function defined in subroutine

write(*,*)'Lower Limit, Upper Limit, Tolerance'

write(,)'Lower Limit, Upper Limit, Tolerance'

write(,)'Initial approx, Tolerance, Max iteration'