Ice Breaker : An application of Newton raphson method of finding roots

INTRODUCTION

Every system and process are generally mathematically represented by system of equations. These equations could be very simple ex: x+3 = 0 and finding root of x is straight forward. The value at which function as zero value or when function crosses x axis is called as the roots of the equation. As the order of the equation increases, number of times the equation crosses over the x axis increases and hence the roots increase.

Sometimes, the equation that has to be solved could be amalgm of Exponential, Logarithmic terms in addition to algebraic equations and are termed as Transcendental equations. For these equations, finding roots analytically could be stressing and next to impossible. In these scenario, a numerical iteration is used, where the guess value is fed into x method and the result obtained is used as the new guess value for next iteration until the value of the function is within the tolerance band.

There are various schemes such as Regular Falsi method, Secant method, Newton raphson method. In this work, we would try to solve Ice Breaker equation by applying Newton Raphson method

Newton raphson method: Xnew = Xold - f(Xold)/f'(Xold), In each successive iteration, previous X value is replaced with newly found one. However, the speed of convergence could be controlled by introducing alpha term in the equation and modified NR equation is as shown:

Xnew = Xold - alpha*f(Xold)/f'(Xold)

ICE BREAKER EQUATION:

Ice breaker, a slow moving air vehicle which induces circumferential crack in ice by bending moments
Pressure of air cushion is sufficient to depress water and force air below icefield and breaks ice. However the phenomenon is function of thickness of ice(h), Width of air cushion(r), Breaking stress(sigma), Pressure(p).

Python Code:

This problem is divided into two parts. In first part, since the value of NR method can change or reduce or increase iteration based on alpha value chosen, at constant value of ice thickness, optimal alpha value is found. This value is then used in the second part. In second part, ice thickness is varied and pressure values are obtained.

The ice breaker equation and its derivative are defined in the function at the very start. NR method is intialised, Ice breaker conditions are passed and NR method is executed until value of function is within tolerance specified. The results are then plotted

# -*- coding: utf-8 -*-
"""
Created on Mon Jul 13 11:41:56 2020

@author: sweekar
This is coded by Sweekar, please reach out at [email protected]
This code is to find Minimum pressure at the ice thickness specified
Ice breaker, a slow moving air vehicle which induces circumferential crack in ice by bending moments
Pressure of air cushion is sufficient to depress water and force air below icefield
"""
import matplotlib.pyplot as plt
import numpy as np

def f(h,r,beta,sigma,p_guess):
    term1 = 1-pow(beta,2)
    term2 = ((0.4*h*pow(beta,2))-(sigma*pow(h/r,2)))
    term3 = (pow(sigma,2)*pow(h/r,4))/3
    term4 = pow(((sigma*pow(h/r,2))/3),3)
    return term1*pow(p_guess,3)+term2*pow(p_guess,2)+term3*p_guess-term4

def fprime(h,r,beta,sigma,p_guess):
    term1 = 1-pow(beta,2)
    term2 = ((0.4*h*pow(beta,2))-(sigma*pow(h/r,2)))
    term3 = (pow(sigma,2)*pow(h/r,4))/3
    return 3*term1*pow(p_guess,2)+2*term2*p_guess+term3

#Inputs
hcheck = 0.6
r = 40
beta = 0.5
sigma = 150

#Initialising Newton raphson solver
p_guess = 80
alpha = np.linspace(0.05,2,40)
count = 0
tol = 1e-6
countmat = []
p = []

for j in range(0,len(alpha)):
    h = hcheck
    while (abs(f(h,r,beta,sigma,p_guess))>tol):
        p_guess = p_guess - alpha[j]*(f(h,r,beta,sigma,p_guess)/fprime(h,r,beta,sigma,p_guess))
        count = count + 1
    countmat.append(count)
    p.append(p_guess)
    count = 0
    p_guess = 80

#Plotting for Alpha dependency
plt.figure(1)
plt.plot(alpha,countmat,color='red')
plt.xlabel('Alpha')
plt.ylabel('Iterations')
plt.title('Alpha dependency')
plt.savefig('Alpha Dependency')
plt.figure(2)
plt.plot(alpha,p,color='blue')
plt.xlabel('Alpha')
plt.ylabel('Pressure')
plt.title('Pressure dependency on alpha')
plt.savefig('Pressure dependency on alpha')

#Checking for Ice thickness and minimum pressure point
#Choosing optimal alpha either as 0.85 or 1.25
#Choosing under relaxation, aplha optimal is 0.85
alphaopt = 0.85
pressure = []
iteramat = []
itera = 0
h = [0.6,1.2,1.8,2.4,3,3.6,4.2]
for i in range(0,len(h)):
    while(abs(f(h[i],r,beta,sigma,p_guess))>tol):
        p_guess = p_guess-alphaopt*(f(h[i],r,beta,sigma,p_guess)/fprime(h[i],r,beta,sigma,p_guess))
        itera = itera+1
    iteramat.append(itera)
    pressure.append(p_guess)
    itera = 0
    p_guess = 80

#Plotting for Height dependency
plt.figure(3)
plt.plot(h,pressure,color='red')
plt.xlabel('Thickness')
plt.ylabel('pressure')
plt.title('pressure vs height')
plt.savefig('pressure vs height')
plt.figure(4)
plt.plot(h,iteramat,color='blue')
plt.xlabel('Thickness')
plt.ylabel('Iterations')
plt.title('Check on aplha chosen')
plt.savefig('Thickness vs iterations')

#Printing values
print('  Thickness  |  Pressure   | Iterations  ')
for k in range(0,7):
    print('  %f   |  %f   |   %d  '%(h[k],pressure[k],iteramat[k]))
print('minimum value of pressure is %f psi and occurs at height of %f feet'%(pressure[0],h[0]))
    

RESULTS AND PLOTS:

Part 1: Alpha value is varied between 0.05 and 2.0 in steps of 0.05, Convergence plot is obtained and variation of pressure values is obtained as shown

From these graphs, it could be noted that as alpha increases, number of iterations required to reach the value exponentially reduces and after certain alpha value, no significant improvement in iterations is seen whereas pressure starts to osciallte and become unstable.

Hence the alpha value has to be between 0.75 and 1.3. Choosing underrelaxation, alpha between 0.75 and 1, least value of pressure could be seen at 0.85 with relatively stable in range of 0.8 to 0.9.

Hence 0.85 is chosen as the optimal alpha

Part 2:

In this step, Thickness is varied and pressure value is noted

From this it could be noted that least pressure occurs at the minimum thickness which in this case is 0.6feet.

Or least value of pressure is 0.5129psf

DISCUSSION ON THE RESULTS

1) The correcteness of value of alpha chosen can be seen from second part based on number of iterations taken to converge

The iterations reduces smoothly and no significant change in total iterations taken and hence the alpha chosen is right

2) As the thickness increases, the pressure requried increases which is logical

3) since h and r are in same units, no need to change unit of sigma and p would be in same units as that of sigma

4) Constraint in the equation is that thickness of ice has to be greater than zero. If ice thickness is zero, pressure is zero and that becomes minimum. Mathematically correct, however, it defeats the purpose of having ice breaker in first place. Hence constraint is h>0