Newton-Raphson method using Python

Aim-

Breaking Ice with Air cushion Vehicle - Finding Minimum Pressure with Newton-Raphson method

Theory- The Newton-Raphson method (also known as Newton's method) is a way to quickly find a good approximation for the root of a real-valued function f ( x ) = 0 f(x) = 0 f(x)=0. It uses the idea that a continuous and differentiable function can be approximated by a straight line tangent to it.

Provided Input-

For this challenge, we will be writing a python script to solve for the minimum value of pressure using the newton raphson method.

where p denotes the cushion pressure, h the thickness of the ice field, r the size of the air cushion, σ the tensile strength of the ice and β is related to the width of the ice wedge.

Take β = 0.5, r = 40 feet and σ = 150 pounds per square inch (psi).

Do not convert to the metric system, however, make sure the units are consistent before proceeding further with the problem.

i need to,

1. Use the newton raphson method to find out the value of pressure for h = 0.6 ft.

2. Find the optimal relaxation factor for this problem with the help of a suitable plot.

3. Tabulate the results of p for h = [0.6,1.2,1.8,2.4,3,3.6,4.2] assuming a suitable relaxation factor.

Program Code-

#Newton - Raphson method
import matplotlib.pyplot as plt

def f(p,h):
	return pow(p,3)*(1-pow(beta,2))+(0.4*h*pow(beta,2)-sigma*pow(h,2)/pow(r,2))*pow(p,2)+pow(sigma,2)*pow(h,4)*p/(3*pow(r,4))-pow((sigma*pow(h,2))/(3*pow(r,2)),3)


def fprime(p,h):
	return 3*pow(p,2)*(1-pow(beta,2))+(0.4*h*pow(beta,2)-sigma*pow(h,2)/pow(r,2))*2*p+pow(sigma,2)*pow(h,4)/(3*pow(r,4))


p_guess=10000
alpha=1
tol=1e-3
beta=0.5
r=40
sigma=150
pressure=[]
Thickness=[0.6,1.2,1.8,2.4,3.0,3.6,4.2]

iter=1

for h in Thickness:
	while(abs(f(p_guess,h))>tol):
		p_guess = p_guess  - alpha*(f(p_guess,h)/fprime(p_guess,h))
		iter = iter + 1
	pressure.append(p_guess)

for i in range(0,7):
	if i == 0:
		print(" : Thickness : pressure")
	print(" : ", Thickness[i]," : ", pressure[i])

print("Number of iterations = ",iter)
plt.plot(Thickness,pressure)
plt.title('Pressure Vs Thickness')
plt.xlabel('Thickness')
plt.ylabel('Pressure')
plt.show()

The above program will solve for the value of the cushion pressure for the different thickness of the ice field. It is using the Newton-Raphson method to find the roots of the equation which indicates the pressure of the ice upon the hovercraft. The Newton-Raphson method is an iterative method.

This method requires us to give an initial value, which is basically a guess value given by the solver (In this program it is the P_guess = 10000).

This program uses two functions(,i.e,f(p,h) and fprime(p,h)), which calculates the equation and its derivative for the inputed values,which are used for evaluation by the Newton-Raphson method.

Output-

The graph on the left shows the variation of the pressure relative to the thickness, and the right shows the output for the above-given program.

The graph delineates that the pressure upon the air cushion decreases as the thickness of the ice field decreases. At a certain point on the graph onwards, the value of the pressure is constant or with insignificant variations. From the graph, we can get that the minimum pressure applied to the air cushion is 0.066(psi). And this observed pressure is constant for a range of the variations in the thickness of the ice field. This variation ranges approximately from 2.4 feet to 0.6 feet.

The output at right shows a table that depicts the values of the pressure for the different values of the thickness of the ice field. The value of the minimum pressure approximated by us from the graph conforms with the value output at the table. A total of 35 iterations have taken place until a constant pressure value is reached for a unity relaxation factor.

So from these results, we can conclude that the minimum pressure will depend on the minimum thickness of the ice field. And for a unity relaxation factor, it seems that the deterministic values are varying gradually and its under-relaxing.

Program 2-

#Optimizing the relaxation factor
import matplotlib.pyplot as plt



def f(p,h,beta,sigma,r):
	t1=(1-pow(beta,2))*pow(p,3)
	t2=(0.4*h*pow(beta,2)-(sigma*pow(h,2))/(pow(r,2)))*pow(p,2)
	t3=(p*pow(sigma,2)*pow(h,4))/(3*pow(r,4))
	t4=pow((sigma*pow(h,2))/(3*pow(r,2)),3)
	t=t1+t2+t3-t4
	return t


def fprime(p,h,beta,sigma,r):

	tp1=3*pow(p,2)*(1-pow(beta,3))
	tp2=2*p*(0.4*h*pow(beta,2)-(sigma*pow(h,2))/(pow(r,2)))
	tp3=pow(sigma,2)*pow(h,4)/(3*pow(r,4))
	tp=tp1+tp2+tp3 
	return tp

# Input Values
h = 0.6				
beta = 0.5			
r = 40			
sigma = 150	
p_guess = 10000
Pressure =[]
iterations = []
tol = 1e-3
alpha_array =[]

for alpha in range(1,90):
	alpha = alpha*0.07
	p_guess = 10000
	alpha_array.append(alpha)

	iter = 1
	while (abs(f(p_guess,h,beta,sigma,r))>tol):		 
		p_guess = p_guess - alpha*(f(p_guess,h,beta,sigma,r)/fprime(p_guess,h,beta,sigma,r))	
		iter = iter + 1
	iterations.append(iter)
	Pressure.append(p_guess)


for i in range(0,len(iterations)):
	if i == 0:
		print ("Alpha Number of Iterations")
	print(alpha_array[i],'', iterations[i])

plt.plot(alpha_array, iterations)
plt.xlabel('Relaxation Factor')
plt.ylabel('No. of iterations')
plt.show()

The above program is used for optimizing the relaxation factor used.

For evaluating the optimum relaxation factor we are comparing the relaxation factor relative to the number of iterations for the corresponding relaxation factor.

In this program, there is a constant value usage of the thickness of the ice field.

Output-

The above graph describes the relationship between the relaxation factor relative to its corresponding number of iterations.

The curve depicted in the graph is an exponentially decreasing curve. From the graph we can understand that the no. of iterations is very high for lower values of the relaxation factor, therefore it is acting as an under-relaxation entity(,i.e, it is taking very small interval values for the iterations), but as the value of the relaxation factor increases there is a rapid drop in the no. of iterations(,i.e the values are taken with larger interval gap) and as the relaxation factor reaches an approximate value of 2.5 there is a slight increase again in the no. of iterations, this is indicating the overthrow nature of the interation depending on the alpha value.

Therefore for this case, the optimum value for the relaxation factor is 3.5,infereing from the graph. (,i.e, for a relaxation value of 3.5 the no.of iterations is two)