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

Objective

> To find out the value of pressure for h = 0.6 ft by the below given equation using Newton-Raphson method.

   `p^3(-beta^2)+(0.4hbeta^2-(sigmah ^2)/r^2)p^2+((sigma^2h^4)/(3r^4))p-((sigmah^2)/(3r^2))^3= 0`

   `p` = Cushion pressure

   `beta` = Width of the ice wedge 

   `r` = Size of air cushion

   `sigma` = Tensile strength of ice

   `h` = Thickness of the ice wedge

> Find the optimal relaxation factor for this problem using a suitable plot.

> To tabulate the results of pressure(p) for h = [0.6,1.2,1.8,2.4,3,3.6,4.2] with a suitable relaxation factor.

Newton-Raphson method

        In numerical analysis, Newton's method, also known as the Newton–Raphson method, named after Isaac Newton and Joseph Raphson, is a root-finding algorithm which produces successively better approximations to the roots (or zeroes) of a real-valued function.

The most basic version starts with a single-variable function f defined for a real variable x, the function's derivative f ′, and an initial guess x₀ for a root of f. If the function satisfies sufficient assumptions and the initial guess is close, then

is a better approximation of the root than x₀. Geometrically, (x₁, 0) is the intersection of the x-axis and the tangent of the graph of f at (x₀, f (x₀)): that is, the improved guess is the unique root of the linear approximation at the initial point. The process is repeated as

until a sufficiently precise value is reached. This algorithm is first in the class of Householder's methods, succeeded by Halley's method. The method can also be extended to complex functions and to systems of equations.

Program

#Import required modules
import math
import matplotlib.pyplot as plt
import numpy as np
from   tabulate import tabulate

#Function definition
def func (p):
	t1 = pow(p,3)*(1-pow(beta,2))
	t2 = (0.4*h*pow(beta,2)-(sig*pow(h,2))/pow(r,2))*pow(p,2)
	t3 = ((pow(sig,2)*pow(h,4))/(3*pow(r,4)))*p
	t4 = pow(((sig*pow(h,2))/(3*pow(r,2))),3)
	return(t1+t2+t3-t4)
def funcp (p):
	t1 = (3*pow(p,2))*(1-pow(beta,2))
	t2 = 2*(0.4*h*pow(beta,2)-(sig*pow(h,2))/pow(r,2))*p
	t3 = ((pow(sig,2)*pow(h,4))/(3*pow(r,4)))
	return(t1+t2+t3)



#Case1 - To find minimum pressure(p) required when h = 0.6
#Inputs
p    = 30
h    = 0.6
beta = 0.5
r    = 40
sig  = 21600
tol  = 1e-4
itr  = 1
while(abs(func(p))>tol):
	p = p - (func(p)/funcp(p))
	itr = itr+1
print(itr)
print(p)



#Case2 - To find optimal relaxation factor(alpha) 
alpha= np.linspace(0.01,1.99,200)
itr  = []
for a in alpha:
	p = 30
	i = 1
	while(abs(func(p))>tol):
		p = p - (a*(func(p)/funcp(p)))
		i = i+1
	itr.append(i)
headers_1 = ['Alpha','Iterations']
data_1    = tuple(zip(alpha,itr))
print(tabulate(data_1,headers=headers_1,tablefmt='grid'))



#Case3 - To find minimum pressure for different heights
height   = [0.6,1.2,1.8,2.4,3,3.6,4.2]
alp      = 1.124
pressure = []
for h in height:
	p = 30
	while(abs(func(p))>tol):
		p = p - (alp*(func(p)/funcp(p)))
	pressure.append(p)
headers_2 = ['Height','Pressure']
data_2    = tuple(zip(height,pressure))
print(tabulate(data_2,headers=headers_2,tablefmt='grid'))



#Plotting Iteration vs Relaxation factor 
plt.figure()
plt.plot(alpha,itr)
plt.xlabel('Relaxation factor')
plt.ylabel('Iteration')
plt.title('Iteration vs Relaxation factor')

#Plotting Pressure vs Height
plt.figure()
plt.plot(height,pressure)
plt.xlabel('Height')
plt.ylabel('Pressure')
plt.title('Pressure vs Height')
plt.show()

Program explanation

> Required modules are imported

> Two functions are defined

> The first function is used to solve the given equation(`f`) by splitting it into four terms

    `t1 = p^3(-beta^2)`

    `t2 = (0.4hbeta^2-(sigmah ^2)/r^2)p^2`

    `t3 = ((sigma^2h^4)/(3r^4))p`

    `t4 = ((sigmah^2)/(3r^2))^3`

> Then the equation is differentiated with respect to p

> The second function is defined accordingly to solve the derivative of the equation(`f'`)

> The defined functions are reused in the same program for all the below cases

Case-1 (To find the minimum pressure)

> The pressure value is assumed to be.

    `p` = 30

> The other inputs given are as follows.

    `beta` = 0.5

    `r` = 40 ft

    `sigma` = 150 * 144 = 21600 

    `h` = 0.6

    tol = 1e-4

> An interative loop is created with the given tolerance and function.

> The loop calculates the minimum pressure value by Newton Raphson method using the defined functions.

> And so, for case 1 the minimum pressure value required to break the given thickness of ice is acquired.

Case-2 (To find optimal relaxation factor)

> The relaxation factor(`alpha`) is used to slow down or pase up the convergence of the value.

> Without an optimal relaxation factor the iterations will be increased.

> The inputs for this case is same as first case but, an additional variable is included which is `alpha` as relaxation factor.

> Alpha is given from (0.01 - 1.99).

> Using linspace command the alpha value is split with number of samples as 200.

> To get the iteration value for each alpha sample an open array is assigned to variable itr.

> A for loop is made with all the samples in alpha.

> The first iteration is given as i = 1

>  Using while loop the minimum value of pressure is calculated by Newton Raphson method with the relaxation factor `alpha`.

    `p = p - alpha((f(p))/(f'(p)))`

> As the loop runs the interation value changes accordingly since i = i+1

> By tabulating the alpha sample to the interation array the optimal relaxation factor is found.

> The Iteration vs alpha is plotted.

Case-3 (To find minimum pressure for different thickness of ice)

> The different thickness of ice is given in an array.

> The opltimal relaxation factor is given as `alpha = 1.124`

> The other inputs are considered to be same as the previous cases.

> To store the pressure values an open array for pressure is created.

> A for loop is made with the values given in thickness array.

> Inside the for loop a while loop is created by which the minimum pressure value is calculated for all the thickness values by Newton Raphson method using the defined functions.

> The calculated pressure values are stored in the pressure array using the command, append.

> Then the pressure array and thickness array are tabulated and plotted. 

Results of the program

Case-1

> From the above image it is seen that the minimum pressure p = 4.2390 when thickness of ice h = 0.6ft 

> The minimum pressure is found using Newton Raphson method in 11 iterations

Case-2

> The above table is a part of the tabulation from the output window in which left row denotes alpha values, and the right row denotes the number of iterations

> It is seen that for alpha value of 1.1243 the interation value is least

> The above given plot show the number of iterations for all the alpha values

> From the table and the plot it is seen that the optimal relaxation factor is 1.1243 

Case-3

> The above table shows the minimum pressure required to break the ice according to their thickness

> It is seen that as the thickness of the ice increases the pressure also increases 

> As the thickness increases the cross sectional area of the ice increases therefore, the pressure required to break the ice also increases

> The same is also seen in the plot shown above where pressure required to break the ice increases as the thickness increases.