Breaking Ice with Air cushion Vehicle - Find minimum pressure with Newton-Raphson method

Here we are going to sove the air cushion problem in python using the newton raphson method.

The data and formula for this are given below

• if we solve for the different height and using the relaxation factor 1 we get the following result show in the figure.

the code for the different hieight is also give where we get this result after runnning this code.

# importing the modules
import numpy as np
import matplotlib.pyplot as plt

# difining the function
def f(p,r,sigma,beta,h):
	term1 = pow(p,3)*(1-pow(beta,2))
	term2 = ((0.4*h*pow(beta,2))-((sigma*pow(h,2))/(pow(r,2))))*pow(p,2)
	term3 = (pow(sigma,2)*pow(h,4)*p)/(3*pow(r,4))
	term4 = pow(((sigma*pow(h,2))/(3*pow(r,2))),3)
	term = term1 + term2 + term3 - term4
	return term


# defining the derivative of the function
def fdesh(p,r,sigma,beta,h):
	
	term1 = 3*pow(p,2)*(1-pow(beta,2))
	term2 = ((0.4*h*pow(beta,2))-((sigma*pow(h,2))/(pow(r,2))))*2*p
	term3 = (pow(sigma,2)*pow(h,4))/(3*pow(r,4))
	term = term1 + term2 + term3
	return term


#variable for calculation
# gussed value
xg = 500
# relaxation factor
alpha = 1
# tolerance
tol = 1e-3

# variable of the model
r=40
sigma=150
beta=0.5
height=[0.6,1.2,1.8,2.4,3.0,3.6,4.2]

# creating the table heading
print('Height   pressure')
# initialing the pressure array
pressure=[]
#for the different height we have to make for loop to calculate the pressure
for h in height:
	while(abs(f(xg,r,sigma,beta,h))>tol):
		# formula of the newton raphson method
		xg = xg - alpha*(f(xg,r,sigma,beta,h)/fdesh(xg,r,sigma,beta,h))
	pressure.append(xg)
	print(h,' -- ',xg)

# making the plot
plt.plot(height,pressure)
plt.show()

Below is the graph of the pressure vs height

Explaination for the code---

-first we import the module numpy and matplotlib for the calculation purpose and making the plot

-now we difine the original function for which we wnat to sove the variable it will require some inlput like gusse value , radius , sigma , beta and height. to call the function we will require this defination

-we will also require the derivative of the functon in the newton raphson method so we also define this next

-now we will declare the different variable. first is guess value which we to give initiallly to start the newton raphson method iterations we declared as 500.relacation fator alpha as 1.tolerance for the answer is declared as 1e-3.

-Then we store the pressure value which we get through the newton raphson method now if the value is more than the tolerance value then it will loop agin into the while loop. it will execute until the we get the value less then th tolerance value

-after that we solve for the different height given in the height array

-next we print the value of the height and the pressure and then plot the pressure vs height

-The image of the plot and the table is given above

• Now we will try to find the optimum relaxation factor in this problem

In here we will try to vary the relaxation factor from 0.1 to 5 devided into 200 parts in between and we will keep track on the iteration for the calculation

Now optimum relaxation factor is one which have a minimum iteration so will make program to calculate the iteration and make plot of this the program is given below.

# importing the modules
import numpy as np
import matplotlib.pyplot as plt

# defining the function
def f(p,r,sigma,beta,h):
	
	term1 = pow(p,3)*(1-pow(beta,2))
	term2 = ((0.4*h*pow(beta,2))-((sigma*pow(h,2))/(pow(r,2))))*pow(p,2)
	term3 = (pow(sigma,2)*pow(h,4)*p)/(3*pow(r,4))
	term4 = pow(((sigma*pow(h,2))/(3*pow(r,2))),3)
	term = term1 + term2 + term3 - term4
	return term

#defining the derivative of the function
def fdesh(p,r,sigma,beta,h):
	term1 = 3*pow(p,2)*(1-pow(beta,2))
	term2 = ((0.4*h*pow(beta,2))-((sigma*pow(h,2))/(pow(r,2))))*2*p
	term3 = (pow(sigma,2)*pow(h,4))/(3*pow(r,4))
	term = term1 + term2 + term3
	return term


#variable for calculation
xg = 500
alpha = np.linspace(0.1,5,200)
tol = 1e-3

# variable of the model
r=40
sigma=150
beta=0.5
h=0.6

#counter variable and itteration array
ct = 1
itter=[]

#creating the 'for' loop for iteration calculation
for a in alpha:
	ct=1
	while(abs(f(xg,r,sigma,beta,h))>tol):
		xg = xg - a*(f(xg,r,sigma,beta,h)/fdesh(xg,r,sigma,beta,h))
		ct = ct + 1
	xg=500
	itter.append(ct)

#plotting part	
plt.plot(alpha,itter)
plt.xlabel('Relaxation factor')
plt.ylabel('Iteration')
plt.show()

Explaination for the code

-first of all we will imprt the module required for the program these are the numpy and matplotlib

-Then we will define the function and its derivative in starting

-Then we declare the variable for the calculation requirement like guess value,relaxation factor and toleraance value for the answer

-Declaring the variable of the model like thickness, sigma , beta and height.

-Now to keep track of the iteration we have to declare the counter variable and iteration araay, itetation is array because each value of the relaxation factor we have to calculate the iteration.

-now we create for loop for alpha wich is the relaxation factor and declare counter as 1 at starting the for the value grater then the tolerance we calculate the guess value as explained in the above program and counter value is increased by one and so on goes the while loop. When the value is less than the tolerance we go outside the while loop. and we again declare the guess value as 500 and save the iteration in the array.

-Now we plot the graph of the iteratio vs relaxation factor.

Resuts of this program discussed below

So the optimum relaxation factor is 3, where the iteration is minimum as show in the above graph