To determine minimum pressure while breaking ice with air cushion vehicle using Newton-Raphson method in Python

AIM

1) To use the Newton-Raphson method to find value of pressure at 0.6 feet

2) Find the optimal relaxtion factor with the help of suitable plot.

3) To determine pressure for given range of ice thickness

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. It uses the idea that a continuous and differentiable function can be approximated by a straight line tangent to it.

It is an iterative technique. It is principally used in the situation where the equation is complex to solve and we need to find the root for the equation for one of the variables involved in the equatioin. In order to do so, we guess a value for which the equation would be 0. Then, an iterative process starts and this guess value is constantly modified in order to reach the actual answer.

At the guess value, the function is assumed to be linear irrespective of its nature at that point. The point at which the straight line intersects the line of f(x) = 0 is projected on the curve perpendicularly and this point is taken as the new guess value. This process is repeated till the guess value becomes constant to the desired tolerance.

Now, from the equation of the straight line assumed at the guess value, we can say that

`(df)/(dx) = (f(x_(guess)) - f(x_(n ew)))/x_(guess - x_(n ew))`

Now, we assume that `f(x_(n ew guess))` = 0

`therefore` the equation simplifies to`fprime(p,h)`

`x_(n ew)` = `x_(guess)` - `fprime(x_(guess))`/`f(x_(guess))`

Now, in order to visualise the change better, we use a constant known as relaxation factor (`alpha`) to vary the rate of change of guess value. When `alpha <1`, this implies that we are being cautios and is known as under relaxation. Here, we want to vary the guess value at the slower rate than it is intended to be. It is helpful when the function is complex and a number of affecting parameters are involved.

`alpha` can also be greater than 1. This is known as over relaxation. This can help in getting to the result in lesser time and number of iterations. However, too much over relaxation can cause overshooting, meaning that the change in guess value is so rapid that the solver reaches the good value but skips it and reaches a farther vaue in the next iteration. This way, it oscillates about the good value before reaching to the result and takes a lot of time. Hence, relaxation factor should be suitably provided.

OBJECTIVE

The objective of this program is to use the Newton-raphson method to determine the cushion pressure for differet ice thickness and also to determine suitable relaxation factor for a given ice thickness. The function relating various parameters involved is given below:

In order to use the Newton-Raphson, first we define a function f(p,h) and then we differentiate the function with respect to p in oreder to find the minimum value of pressure and define `fprime(p,h)`.

For varying ice thicknes, we define an array containing the given thickness values and provided the rest of the inputs. We set the guess value for pressure = 1 and tolerance upto 4 decimal places.

using for loop, we provide the different values of thickness to the function and store the corresponding values of minimum pressure in an array using append command.

Here, while loop is used to make sure that the loop runs till the value of the function becomes constant upto the given tolerance.

Once the loop is exit, the required values are printed and the plot between pressure and thickness of ice is made.

Determining optimum value of `alpha`

To determine optimum relaxation factor, similar procedure is used but this time, thickness is kept constant to 0.6 feet and an array for relaxation factor values is created.

Again, for loop is used using different values in Rf array but this time, number of iterations required to get the result is obtained and stored in an array. While loop is used for the similar purpose as above.

Once it is completed, the value of `alpha` corresponding to minimum number of iterations is determined and the pressure value is found using that value.

PROGRAM

# To determine minimum pressure using Newton-Raphson method
import math
import matplotlib.pyplot as plt
import numpy as np

# Function
def f(p,h):
	a = pow(p,3)*(1-pow(beta,2))
	b = (0.4*h*pow(beta,2) - sigma*pow(h,2)/pow(r,2))*pow(p,2)
	c = (pow(sigma,2)*pow(h,4)/(3*pow(r,4)))*p 
	d = sigma*pow(h,2)/(3*pow(r,2))
	return  a + b + c - pow(d,3)

# Derivative w.r.t p
def fprime(p,h):
	a = 3*pow(p,2)*(1-pow(beta,2)) 
	b = 2*(0.4*h*pow(beta,2) - sigma*pow(h,2)/pow(r,2))*p
	c =  pow(sigma,2)*pow(h,4)/(3*pow(r,4))
	return a + b + c

beta = 0.5
#size of air cushion
r = 40 #feet
# thickness of ice
h = np.linspace(0.6,4.2,7)
#tensile strength of ice
sigma = 150 #psi
sigma = sigma*144


p = 1
pressure = []
#alpha = 0.4
tol = 1e-04
iter = 1

# Perforing Newton-Raphson
for t in h:
	while(abs(f(p,t))>tol):
		p = p - (f(p,t)/fprime(p,t))
		iter = iter + 1
	pressure.append(p)
pressure = np.array(pressure)/144 #converting pressure to psi

print(' height     pressure')
for i in range(0,7):
	print(' ', h[i], '   :   ', pressure[i])

# Ploting
plt.plot(h,pressure)
plt.xlabel('height')
plt.ylabel('minimum pressure')
plt.title('Minimum pressure vs height')
plt.grid('on')
plt.show()

To determine the optimum value for relaxation factor

The optimum relaxation factor is that for which minimum iterations are required to get to the final answer.

# To determine optimum value of relaxation factor
import math
import matplotlib.pyplot as plt
import numpy as np

# Function
def f(p,h):
	a = pow(p,3)*(1-pow(beta,2))
	b= (0.4*h*pow(beta,2) - sigma*pow(h,2)/pow(r,2))*pow(p,2)
	c = (pow(sigma,2)*pow(h,4)/(3*pow(r,4)))*p 
	d = sigma*pow(h,2)/(3*pow(r,2))
	return  a + b + c - pow(d,3)

# Derivative of function w.r.t p
def fprime(p,h):
	a = 3*pow(p,2)*(1-pow(beta,2)) 
	b= 2*(0.4*h*pow(beta,2) - sigma*pow(h,2)/pow(r,2))*p
	c =  pow(sigma,2)*pow(h,4)/(3*pow(r,4))
	return a + b + c

# inputs
beta = 0.5
# thickness of ice
h = 0.6
# tensile strength
sigma = 150*144 
# Size of air cushion
r = 40
# relaxation factor
alpha = np.linspace(0.1,1.9,400)

iterations = []
tol = 1e-04

# performing Newton-raphson for different relaxation fctors
for j in alpha:
	i = 1
	p_guess = 100
	while(abs(f(p_guess,h))>tol):
		p_guess = p_guess - j*(f(p_guess,h)/fprime(p_guess,h))
		i = i + 1
	iterations.append(i)
print(min(iterations))

# determining optimum relaxation factor
for i in range(len(iterations)):
	if iterations[i] == min(iterations): # alpha corresponding to min. iterations
		print(alpha[i])

# Plotting
plt.plot(alpha,iterations)
plt.xlabel('relaxation factor')
plt.ylabel('iterations')
plt.title('iterations vs relaxation factor')
plt.grid('on')
plt.show()

RESULTS

OBSERVATION AND CONCLUSION

From the plots, we can see that the value of pressure increases with increase in the thickness of ice. Now, considering the plot of relaxation factor with number of iterations, we see that the number of iterations decreases with increase in `alpha` but only up to a certain value, beyond which it starts increasing. Thus, we get an optimum value of alpha corresponding to which, the number of iterations are minimum to reach the results. If the value is increased any further, number of iterations required starts increasing as the guess value oscillates about the good value.

The minimum value of pressure corresponding to given thickness of ice is determined and displayed using the Newton-Raphson method and the optimum value of relaxation factor is determined.