Newton-Rhapson method

The Newton-Rhapson method is an iterative technique that is used to solve for the roots of a function. In this case, we use the Newton-Rhapson method to solve for the cushion pressure required to cut a certain thickness of ice using an air cushion vehicle.

The thickess, denoted by h in the equation is varies from 0.6 to 4.2 and the corresponding pressures are found using the Newton-Rhapson method.

A second code was written to figure out the iterations required for convergence for relaxation factors ranging from 0 to 2.

Code that measures pressure required for different ice thickness:

import matplotlib.pyplot as plt

#Defining the function and the derivative of the function
def f(p):
	term1 = (1-pow(beta,2))
	term2 = (0.4*h*(pow(beta, 2))) - (sig*(pow(h,2))/(pow(r,2)))
	term3 = (pow(sig,2))*(pow(h,4))/(3*(pow(r,4)))
	term4 = pow((sig*pow(h,2))/(3*(pow(r,2))),3)

	return pow(p,3)*term1 + pow(p,2)*term2 + p*term3 - term4


def fp(p):
	term1 = (1-pow(beta,2))
	term2 = (0.4*h*(pow(beta, 2))) - (sig*(pow(h,2))/(pow(r,2)))
	term3 = (pow(sig,2))*(pow(h,4))/(3*(pow(r,4)))
	term4 = pow((sig*pow(h,2))/(3*(pow(r,2))),3)

	return 3*pow(p,2)*term1 + 2*p*term2 + term3


#Defining the constants
beta = 0.5
r = 40 #size of air cushion
sig = 150*144 #tensile strength of ice
a = 1 #relaxation factor
p = 40
tol = 1e-4

iter = 1

h1 = []
p1 = []
for h in range(1,8):
	h = h*0.6
	h1.append(h)
	
	while (abs(f(p)) > tol):
		p = (p - a*(f(p)/fp(p)))
		iter = iter + 1

	print(p, iter)
	p1.append(p)

	p = 40
	iter = 1	

plt.figure()
plt.plot(h1,p1)	
plt.xlabel(\'Thickness of the ice\')
plt.ylabel(\'Cushion pressure required\')
plt.savefig(\'phplot.png\')
plt.show()

Results:

The relationship between p and h is tabulated and shown in the form of a graph as well. It is clear that as thickness of ice increases, the pressure required increases as well. The relationship is not linear but parabolic with the pressure required increasing exponentially.

Code that finds the iterations required for different relaxation factors:

import matplotlib.pyplot as plt

#Defining the function and the derivative of the function
def f(p):
	term1 = (1-pow(beta,2))
	term2 = (0.4*h*(pow(beta, 2))) - (sig*(pow(h,2))/(pow(r,2)))
	term3 = (pow(sig,2))*(pow(h,4))/(3*(pow(r,4)))
	term4 = pow((sig*pow(h,2))/(3*(pow(r,2))),3)

	return pow(p,3)*term1 + pow(p,2)*term2 + p*term3 - term4


def fp(p):
	term1 = (1-pow(beta,2))
	term2 = (0.4*h*(pow(beta, 2))) - (sig*(pow(h,2))/(pow(r,2)))
	term3 = (pow(sig,2))*(pow(h,4))/(3*(pow(r,4)))
	term4 = pow((sig*pow(h,2))/(3*(pow(r,2))),3)

	return 3*pow(p,2)*term1 + 2*p*term2 + term3

#Defining the constants
beta = 0.5
r = 40
sig = 150*144
h = 0.6
p = 40
tol = 1e-4

iter = 1

a1 = []
it = []
for a in range(1,200):
	a = a*0.01
	a1.append(a)
	
	while (abs(f(p)) > tol):
		p = (p - a*(f(p)/fp(p)))
		iter = iter + 1

	print(p, iter, a)
	it.append(iter)
	p = 40
	iter = 1	

plt.figure()
plt.plot(a1,it)	
plt.xlabel(\'alpha\')
plt.ylabel(\'No. of iterations\')
plt.savefig(\'aviter.png\')
plt.show()

Relaxation factor v alpha plot:

The relationship between the relaxation factor and the no. of iterations is a curious one. It seemingly takes up the same number of iterations to converge to the value of p for a relaxation factor from about 0.8 to 1.5. Since this is an under relaxation study, the optimal relaxation factor can be taken as 0.99 or 1. However, the range of relaxation that gave the least number of iterations was from 0.99 to 1.12 with 12 iteration each except for at 1.1 which converged at the 11th iteration.