JAGANNATH UNIVERSITY

Department of Mathematics

ASSIGNMENT
Name : Saima Tasnim
ID : B180302053
Batch : 15th

Course : FORTRAN Programming Lab II

Course Topic : 8 Numerical Methods
Course Code : MTHL 3206
Course Teacher : Goutam Kumar Saha
(Assistant Professor)

Submission Date: 01/10/2023

 Bisection Method

Bisection is a root-finding method applies to any continuous functions with two known values of
opposite signs. The interval defined by these two values is bisected and a sub-interval in which
the function changes sign is selected. This sub-interval must contain the root. These two steps are
repeatedly executed until the root is in the form of the required precision level. This method is
also called as interval halving method, the binary method, or the dichotomy method.
Let 𝑓 be a continuous function defined on an interval [𝑎,𝑏] where 𝑓(𝑎) and 𝑓(𝑏) have opposite
signs. The bisection method is applicable for solving the equation 𝑓(𝑥)=0 for a real variable 𝑥. At
each step, the interval is divided into two parts/halves by computing the midpoint, 𝑐=(𝑎+𝑏)/2,
and the value of 𝑓(𝑐) at that point. Unless the root is 𝑐, there are two possibilities:𝑓(𝑎) and 𝑓(𝑐)
have opposite signs and bracket a root,𝑓(𝑐) and 𝑓(𝑏) have opposite signs and bracket a root. One
of the sub-intervals is chosen as the new interval to be used in the next step. This process is
carried out again and again until the interval is sufficiently small. If 𝑓(𝑎) and 𝑓(𝑐) have opposite
signs, then the value of 𝑏 is replaced by 𝑐. If 𝑓(𝑏) and 𝑓(𝑐) have opposite signs, then the value of
𝑎 is replaced by 𝑐.In the case that 𝑓(𝑐)=0,𝑐 will be taken as the solution and the process stops.

Example
Here, we have bisection method example problems with solution.
Question. Determine the root of the equation, f(x)=x3–x–2 for x∈[1,2].
Solution.
Given, f(x)=x3–x–2 for x∈[1,2]
Take, a=1,b=2
f(1)=1³–1–2=−2
f(2)=2³–2–2=+4
As the function is continuous, a root must lie within [1, 2].
1st Iteration: a1=1,b1=2,
c1=(2+1)/2=1.5
Hence, the function value at the midpoint is,
f(c1)=(1.5)3–(1.5)–2=−0.125
As f(c) gives a negative value, the value of a is replaced by c

The iterations are concisely summarized into a table below:

From the above table, it can be pointed out that, after 13 iterations, it becomes apparent that the
function converges to 1.521, which is concluded as the root of the polynomial.

Algorithm
To find a solution to f (x) = 0 given the continuous function f on the interval [a, b], where f (a)
and f (b) have opposite signs:
INPUT endpoints a, b; tolerance TOL; maximum number of iterations N0.
OUTPUT approximate solution p or message of failure.
Step 1 Set i = 1; FA = f (a).
Step 2 While i ≤ N0 do Steps 3–6.
Step 3 Set p = a + (b − a)/2; (Compute pi.) FP = f ( p).
Step 4 If FP = 0 or (b − a)/2 < TOL then OUTPUT (p); (Procedure completed successfully.)
STOP. Step 5 Set i = i + 1.
Step 6 If FA · FP > 0 then set a = p; (Compute ai, bi.) FA = FP else set b = p. (FA is
unchanged.) Step 7 OUTPUT (‘Method failed after N0 iterations, N0 =’, N0); (The procedure
was unsuccessful.) STOP.

Program
!sin x = 1 - x**2 correct upto five decimal places using bisection
f(x)=sin(x) - 1.0 + x**2
50 print*,'enter the trial solutions,tolerance and number of iteration'
read*,a,b,tol,n
if(f(a)*f(b).ge.0.0) go to 30
print*,' n a b p '
do i=1,n
P=(a+b)/2
Print*,i,a,b,p
if(abs(p-a).lt.tol)go to 20
if(f(a)*f(p).ge.0.0)then
a=p
else
b=p
end if
end do
30 print*,' The method fails '
go to 50
20 write(6,70) P
70 Format('The solution is=',f10.5)
stop
end

Output
 Fixed Point Iteration Method

The fixed point iteration method in numerical analysis is used to find an approximate solution to
algebraic and transcendental equations. Sometimes, it becomes very tedious to find solutions to
cubic, bi-quadratic, and transcendental equations; then, we can apply specific numerical methods
to find the solution; one among those methods is the fixed point iteration method.
The fixed point iteration method uses the concept of a fixed point in a repeated manner to
compute the solution of the given equation. A fixed point is a point in the domain of a function g
such that g(x) = x. In the fixed point iteration method, the given function is algebraically
converted in the form of g(x) = x.
Suppose we have an equation f(x) = 0, for which we have to find the solution. The equation can
be expressed as x = g(x). Choose g(x) such that |g’(x)| < 1 at x = x o where xo, is some initial guess
called fixed point iterative scheme. Then the iterative method is applied by successive
approximations given by xn = g(xn – 1), that is, x1 = g(xo), x2 = g(x1), and so on.

Example:
Find the first approximate root of the equation 2x3 – 2x – 5 = 0 up to 4 decimal places.
Solution:
Given f(x) = 2x3 – 2x – 5 = 0
As per the algorithm, we find the value of x o, for which we have to find a and b such that f(a) < 0
and f(b) > 0
Now, f(0) = – 5
f(1) = – 5;f(2) = 7
Thus, a = 1 and b = 2
Therefore, xo = (1 + 2)/2 = 1.5
Now, we shall find g(x) such that |g’(x)| < 1 at x = xo
2x3 – 2x – 5 = 0
⇒ x = [(2x + 5)/2]1/3
g(x) = [(2x + 5)/2]1/3 which satisfies |g’(x)| < 1 at x = 1.5
Now, applying the iterative method xn,= g(xn – 1) for n = 1, 2, 3, 4, 5, …
For n = 1; x1 = g(xo) = [{2(1.5) + 5}/2]1/3 = 1.5874
For n = 2; x2 = g(x1) = [{2(1.5874) + 5}/2]1/3 = 1.5989
For n = 3; x3 = g(x2) = [{2(1.5989) + 5}/2]1/3 = 1.60037
For n = 4; x4 = g(x3) = [{2(1.60037) + 5}/2]1/3 = 1.60057
For n = 5; x5 = g(x4) = [{2(1.60057) + 5}/2]1/3 = 1.60059
For n = 6; x6 = g(x5) = [{2(1.60059) + 5}/2]1/3 = 1.600597 ≈ 1.6006
The approximate root of 2x3 – 2x – 5 = 0 by the fixed point iteration method is 1.6006.

Algorithm
To find a solution to p = g( p) given an initial approximation p0:
INPUT initial approximation p0; tolerance TOL; maximum number of iterations N0.
OUTPUT approximate solution p or message of failure.
Step 1 Set i = 1.
Step 2 While i ≤ N0 do Steps 3–6.
Step 3 Set p = g( p0). (Compute pi.)
Step 4 If | p − p0| < TOL then OUTPUT ( p); (The procedure was successful.) STOP.
Step 5 Set i = i + 1.
Step 6 Set p0 = p. (Update p0.)
Step 7 OUTPUT (‘The method failed after N0 iterations, N0 =’, N0); (The procedure was
unsuccessful.) STOP.

Program
f(x)=sin(x)-x
g(x)=sin(x)
print*,'enter the trial solution , tolerance and number of iteration'
read*,a,tol,n
print*,' n a p'
do i=1,n
p=g(a)
Print*,i,a,p
if(abs(p-a).lt.tol)go to 20
a=p
end do
20 print*,'The solution is=',p
stop
end
Output
 Newton-Raphson Method

The Newton-Raphson Method is referred to as one of the most commonly used techniques for
finding the roots of given equations. It can be efficiently generalized to find solutions to a system
of equations. Moreover, we can show that when we approach the root, the method is
quadratically convergent. Newton Raphson's method is an efficient technique to solve the
equations numerically. It gives us better approximations in terms of solutions.

Let x0 be the approximate root of f(x) = 0 and let x1 = x0 + h be the correct root.
Then f(x1) = 0
⇒ f(x0 + h) = 0….(1)

By expanding the above equation using Taylor’s theorem, we get:

f(x0) + hf1(x0) + … = 0
⇒ h = -f(x0) /f’(x0)

Therefore, x1 = x0 – f(x0)/ f’(x0)

Now, x1 is a better approximation than x0.

Similarly, the successive approximations x2, x3, …., xn+1 are given by

This is called Newton Raphson's formula.

Example:
Find a root of an equation f(x)=x3-x-1 using Newton Raphson method
Solution:
Here x3-x-1=0
Let f(x)=x3-x-1
∴f′(x)=3x2-1

Here

x 0 1 2

-
f(x) -1 5
1

Here f(1)=-1<0andf(2)=5>0
∴ The root lies between 1 and 2
Algorithm
To find a solution to f (x) = 0 given an initial approximation p0:

INPUT initial approximation p0; tolerance TOL; maximum number of iterations N0.

OUTPUT approximate solution p or message of failure.

Step 1 Set i = 1.

Step 2 While i ≤ N0 do Steps 3–6.

Step 3 Set p = p0 − f ( p0)/f ( p0). (Compute pi.)

Step 4 If | p − p0| < TOL then OUTPUT (p); (The procedure was successful.) STOP.

Step 5 Set i = i + 1.

Step 6 Set p0 = p. (Update p0.)

Step 7 OUTPUT (‘The method failed after N0 iterations, N0 =’, N0); (The procedure was
unsuccessful.) STOP.

Program
!find root of x^3-x-1=0 up to 4 decimal places using newton raphson
!D[x^3-x-1, x]

f(x)= x**3 - x - 1
g(x)= 3*x**2 - 1
print*,'enter the trial solution , tolerance and number of iteration'
read*,a,tol,n
print*,' n a p'
do i=1,n
p=a-(f(a)/g(a))
write(6,12) i,a,p
12 format(i3,f10.4,f10.4)
if(abs(p-a).lt.tol)go to 20
a=p
end do
20 write(6,70) p
70 Format('The solution is =', f10.4)
stop
end
Output
 False Position Method

An ancient method of solving an equation in one variable is the false position method (method of
false position) or regula falsi method. In simple words, the method is described as the trial and
error approach of using “false” or “test” values for the variable and then altering the test value
according to the result.
Consider an equation f(x) = 0, which contains only one variable, i.e. x.
To find the real root of the equation f(x) = 0, we consider a sufficiently small interval (a, b)
where a < b such that f(a) and f(b) will have opposite signs. According to the intermediate value
theorem, this implies a root lies between a and b.
Also, the curve y = f(x) will meet the x-axis at a certain point between A[a, f(a)] and B[b, f(b)].
Now, the equation of the chord joining A[a, f(a)] and B[b, f(b)] is given by:

Let y = 0 be the point of intersection of the chord equation (given above) with the x-axis. Then,

This can be simplified as:

Thus, the first approximation is x1 = [af(b) – bf(a)]/ [f(b) – f(a)]

Also, x1 is the root of f(x) if f(x1) = 0.
If f(x1) ≠ 0 and if f(x1) and f(a) have opposite signs, then we can write the second approximation
as:
x2 = [af(x1) – x1f(a)]/ [f(x1) – f(a)]
Similarly, we can estimate x3, x4, x5, and so on.
Example
Find a root for the equation 2ex sin x = 3 using the false position method and correct it to three
decimal places with three iterations.
Solution:
Given equation: 2ex sin x = 3
This can be written as: 2ex sin x – 3 = 0
Let f(x) = 2ex sin x – 3
So, f(0) = 2e0 sin 0 – 3
=0–3
= -3 < 0
And
f(1) = 2e1 sin 1 – 3
= 2e sin 1 – 3
= 1.5747 > 0
That means the root of f(x) lies between 0 and 1.
Let a = 0 and b = 1.
The first approximation = x1 = [af(b) – bf(a)]/ [f(b) – f(a)]
= [0(1.5747) – 1(-3)]/ [1.5747 – (-3)]
= 0.6557
Thus, x1 = 0.6557
Now, substitute x = 0.6557 in f(x).
So, f(0.6557) = 2e0.6557 sin(0.6557) – 3
= 2.3493 – 3
= -0.6507 < 0
As we know, f(1) > 0
That means a root lies between 0.6557 and 1.
Let a = 0.6557
The second approximation = x2 = [af(x1) – x1f(a)]/ [f(x1) – f(a)]
= [0.6557(1.5747) – 1(-0.6507)]/ [1.5747 – (-0.6507)]
= 1.6832/2.2254
= 0.7563
Therefore, x2 = 0.7563
Let us substitute 0.7563 in f(x).
So, f(0.7563) = 2e0.7563 sin(0.7563) – 3
= 2.9239 – 3
= -0.0761 < 0
We know that f(1) > 0
Thus, a root lies between 0.7563 and 1.
Hence, the third approximation = x3 = [af(x2) – x2f(a)]/ [f(x2) – f(a)]
= [(0.7563)(1.5747) – 1(-0.0761)]/ [1.5747 – (-0.0761)]
= (1.1909 + 0.0761)/1.6508
= 1.2670/1.6508
= 0.7675
So, x3 = 0.7675
Therefore, the best approximation of the root up to three decimal places is 0.768 (up to three
decimal places).

Algorithm
To find a solution to f (x) = 0 given the continuous function f on the interval [ p0, p1] where f
( p0) and f ( p1) have opposite signs:
INPUT initial approximations p0, p1; tolerance TOL; maximum number of iterations N0.
OUTPUT approximate solution p or message of failure.
Step 1 Set i = 2; q0 = f ( p0); q1 = f ( p1).
Step 2 While i ≤ N0 do Steps 3–7.
Step 3 Set p = p1 − q1( p1 − p0)/(q1 − q0). (Compute pi.)
Step 4 If | p − p1| < TOL then OUTPUT (p); (The procedure was successful.) STOP.
Step 5 Set i = i + 1; q = f ( p).
Step 6 If q · q1 < 0 then set p0 = p1; q0 = q1.
Step 7 Set p1 = p; q1 = q.
Step 8 OUTPUT (‘Method failed after N0 iterations, N0 =’, N0); (The procedure unsuccessful.)
STOP.
Program
f(x)=cos(x) - x*exp(x)
50 print*,'enter the trial solutions,tolerance and number of iteration'
read*,a,b,tol,n
if(f(a)*f(b).ge.0.0) go to 30
print*,' n a b p '
do i=1,n
p=(a*f(b)-b*f(a))/(f(b)-f(a))
Print*,i,a,b,p
if(abs(p-a).lt.tol)go to 20
if(f(a)*f(p).ge.0.0)then
a=p
else
b=p
end if
end do
30 print*,' The method fails '
go to 50
20 Print*,'The solution is=',p
stop
end

Output
 Newton’s Forward Interpolation Method

Interpolation is the method of finding the value of a function for any intermediate value for an
independent variable, while the process of calculating the value of the function outside the given
range is called extrapolation.
Forward Differences: The differences y1 – y0, y2 – y1, y3 – y2… yn – yn–1 when denoted by
dy0, dy1, dy2… dyn–1 are respectively, called the first forward differences.

Formula

Example
Algorithm:

Step 1: Start the program

Step 2: Read n (No. of arguments)

Step 3: For i = 0 to n – 1

Read x i &y i [0]

End i

Step 4: Construct the Forward Difference Table

For j = 1 to n – 1

For i = 0 to n − 1 – j

y i [j] = y[i + 1][j − 1] − y[i][j − 1]

End i

End j

Step 5: Print the Forward Difference Table

For i = 0 to n – 1

For j = 0 to n − 1 – i

Print y i [j]

End j

Next Line

End i

Step 6: Read a (Point of Interpolation)

Step 7: Assign h = x[1] − x[0] (Step Length)

Step 8: Assign u = (a − x[0])/h

Step 9: Assign sum = y 0 [0] & p = 1.0

Step 10: For j = 1 to n – 1

p = p ∗ (u − j + 1)/j

sum = sum + p ∗ y 0 [j]

End j

Step 11: Display a & sum

Step 12: Stop

Program

real x(10),y(10),a(10),b(10),c(10),d(10)

print*,'Enter number of points'

read*,n

print*,'Enter the values of x and y'

read*,(x(i),i=1,n),(y(i),i=1,n)

do i=1,n-1
a(i)=y(i+1)-y(i)
enddo

do i=1,n-2
b(i)=a(i+1)-a(i)
enddo

do i=1,n-3
c(i)=b(i+1)-b(i)
enddo

do i=1,n-4
d(i)=c(i+1)-c(i)
enddo

print*,'enter x to calculate y'

read*,x0

u=(x0-x(1))/(x(2)-x(1))

yk=y(1)+u*a(1)+(u*(u-1.0)*b(1))/2.0+(u*(u-1.0)*(u-2.0)*c(1))/6.0+
(u*(u-1.0)*(u-2.0)*(u-3.0)*d(1))/24.0
 write(6,20)yk
20 format('The value of y is = ',f10.4)
stop
end

Output
 Newton’s Backward Interpolation Method
Backward Differences: The differences y1 – y0, y2 – y1, ……, yn – yn–1 when denoted by
dy1, dy2, ……, dyn, respectively, are called first backward difference. Thus, the first backward
differences are:

Formula

Example
Algorithm

Step 1: Start the program

Step 2: Read n (No. of arguments)

Step 3: For i = 0 to n – 1

Read x i &y i [0]

End i

Step 4: Construct the Backward Difference Table

For j = 1 to n – 1

For i = j to n – 1

y i [j] = y[i][j − 1] − y[i − 1][j − 1]

End i

End j

Step 5: Print the Backward Difference Table

For i = 0 to n – 1

For j = 0 to i

Print y i [j]

End j

End i

Step 6: Read a (Point of Interpolation)

Step 7: Assign h = x[1] − x[0] (Step Length)

Step 8: Assign u = (a − x[n − 1])/h

Step 9: Assign sum = y n − 1 [0] & p = 1.0

Step 10: For j = 1 to n – 1

p = p ∗ (u + j − 1)/j

sum = sum + p ∗ y n − 1 [j]

End j

Step 11: Display a & sum

Step 12: Stop

Program

real x(10),y(10),a(10),b(10),c(10),d(10)

print*,'Enter number of points '

read*,n

print*,'Enter the values of x and y'

read*,(x(i),i=1,n),(y(i),i=1,n)

do i=n,2,-1
a(i)=y(i)-y(i-1)
enddo

do i=n,3,-1
b(i)=a(i)-a(i-1)
enddo

do i=n,4,-1
c(i)=b(i)-b(i-1)
enddo

do i=n,5,-1
d(i)=c(i)-c(i-1)
enddo

print*,'enter x to calculate y'

read*,x0

u=(x0-x(n))/(x(2)-x(1))
yk=y(n)+u*a(n)+(u*(u+1.0)*b(n))/2.0+(u*(u+1.0)*(u+2.0)*c(n))/6.0+
(u*(u+1.0)*(u+2.0)*(u+3.0)*d(n))/24.0

write(6,20)yk
20 format('The value of y is = ',f10.4)

stop
end
Output
 Lagrange Interpolation Method
Newton’s forward and backward interpolation formulae can be used only when the values
of x are at equidistant. If the values of x are at equidistant or not at equidistant, we use
Lagrange’s interpolation formula.

Let y = f(x) be a function such that f (x) takes the values y 0 , y1 , y2 ,......., yn corresponding to x=
x0 , x1, x2 ..., xn That is yi = f(xi),i = 0,1,2,...,n . Now, there are (n + 1) paired values (x i, yi),i = 0,
1, 2, ..., n and hence f ( x) can be represented by a polynomial function of degree n in x.

Then Lagrange’s formula is

Example

Using Lagrange’s interpolation formula find y(10) from the following table:

Solution:
Here the intervals are unequal. By Lagrange’s interpolation formula we have
Algorithm

Step 1. Start

Step 2. Read number of data (n)

Step 3. Read data Xi and Yi for i=1 ton n

Step 4. Read the value of independent variables say xp whose corresponding value of dependent
say yp is to be determined.

Step 5. Initialize: yp = 0

Step 6. For i = 1 to n

Set p = 1

For j =1 to n

If i ≠ j then

Calculate p = p * (xp - Xj)/(Xi - Xj)

End If

Next j

Calculate yp = yp + p * Yi

Next i

Step 7. Display the value of yp as interpolated value.

Step 8. Stop

Program

real x(10),y(10)
print*,'enter the number of points'
read*,n
print*,'Enter the values of x and y'

read*,(x(i),i=1,n),(y(i),i=1,n)
print*,'enter x to calculate y'
read*,x0
yp=0.0
do j=1,n
yd=1.0
do i=1,n
if(i.eq.j)go to 7
yd=yd*((x0-x(i))/(x(j)-x(i)))
7 enddo

yp=yp+yd*y(j)
end do

write(6,12)yp
12 format('The value is =',f10.4)
stop
end

Output
NEWTON'S DIVIDED DIFFERENCE INTERPOLATION
METHOD
Examples
1. Find Solution using Newton's Divided Difference Interpolation formula

x f(x)

300 2.4771

304 2.4829

305 2.4843

307 2.4871
x = 301

Solution:
The value of table for x and y

x 300 304 305 307

2.482
y 2.4771 2.4843 2.4871
9

Numerical divided differences method to find the solution

Algorithm

Input:

- Data points (x_i, y_i) for i = 0, 1, ..., n

Initialize: - Initialize an empty divided difference table D.

Step 1: Fill in the first column of the divided difference table with y-values.

For i = 0 to n: D[i][0] = y_i

Step 2: Calculate divided differences.

For j = 1 to n: For i = 0 to n - j: D[i][j] = (D[i+1][j-1] - D[i][j-1]) / (x[i+j] - x[i])

Step 3: Construct the interpolating polynomial.

The interpolating polynomial P(x) can be written as:

P(x) = D[0][0] + (x - x_0)D[0][1] + (x - x_0)(x - x_1)D[0][2] + ...

Step 4: Use the polynomial for interpolation. We can evaluate P(x) at any desired x-value to
approximate the corresponding y-value.

Output: - Interpolating polynomial P(x) for the given data points.

Program

program NewtonInterpolation

implicit none

integer, parameter :: max_points = 100

real(8) :: x(max_points), f(max_points),

divided_diff(max_points)

real(8) :: interp_x, interp_result

integer :: n, i, j
 real(8) :: prod

! Read the number of data points

write(*,*) "Enter the number of data points (n):"

read(*,*) n

! Read the data points (x, f(x))

write(*,*) "Enter the data points (x(i), f(i)):"

do i = 1, n

read(*,*) x(i), f(i)

end do

! Calculate the divided differences

divided_diff = f

do j = 2, n

do i = n, j, -1

divided_diff(i) = (divided_diff(i) - divided_diff(i-

1)) / (x(i) - x(i-j+1))

end do

end do

! Enter the value at which you want to interpolate

write(*,*) "Enter the value at which you want to interpolate

(interp_x):"

read(*,*) interp_x

! Perform interpolation

interp_result = divided_diff(1)
 prod = 1.0

do i = 2, n

prod = prod * (interp_x - x(i-1))

interp_result = interp_result + divided_diff(i) * prod

end do

! Display the interpolation result

write(*,*) "Interpolated result at x =", interp_x, "is",

interp_result

end program NewtonInterpolation

Output