Breaking Ice with Air cushion Vehicle - Find minimum pressure with Newton-Raphson method
Objective:
Write a code to solve given equation by using Newton Raphson method.
1.Use relaxation factor of 1 and explain your results.
2.How will you choose the optimum relaxation factor for this case? and what is the optimum relaxation factor for this case? Show with the help of a graph, assume h = 0.6.
3.Explain your results with respect to the appropriate relaxation factor.
Introduction:
Newton-Raphson Method
The Newton-Raphson method is a way to quickly find a good approximation for the root of a real-valued function f(x)=0f(x) = 0f(x)=0. It uses the idea that a continuous and differentiable function can be approximated by a straight line tangent to it.
Suppose you need to find the root of a continuous, differentiable function f(x)f(x)f(x), and you know the root you are looking for is near the point x=x0x = x_0x=x0. Then Newton's method tells us that a better approximation for the root is
x1 = x0 -(f(x0))/(f'(x0))
This process may be repeated as many times as necessary to get the desired accuracy. In general, for any x-value xn the next value is given by
xn+1 = xn -(f(xn))/(f'(xn))
Given:
Determine the minimum cushion pressure needed to break a given thickness of ice
using an air cushion vehicle, Muller (“Ice Breaking with an Air Cushion Vehicle”) in
Mathematical Modeling: Classroom Notes in Applied Mathematics, SIAM 1987) derived
the equation
where,
p = the cushion pressure
β = related to the width of the ice wedge = 0.5
h = thickness of the ice field
r = the size of the air cushion = 40
σ = the tensile strength of the ice = 150 psi= (150*144) psf
code:
Case-1. With relaxation factor 1
import matplotlib.pyplot as plt
import math
from tabulate import tabulate
import numpy as np
# Input
β = 0.5 # related to the width of the ice wedge
r = 40 # the size of the air cushion
σ = 150*144 # the tensile strength of the ice
height = [0.6,1.2,1.8,2.4,3.0,3.6,4.2]
alpha = 1 # Relaxation factor
# Define function
def f(p):
term1 = pow(p,3)*(1-pow(β,2))
term2 = pow(p,2)*((0.4*h*pow(β,2))-((σ*pow(h,2))/pow(r,2)))
term3 = p*(pow(σ,2)*pow(h,4)/(3*pow(r,4)))
term4 = (σ*pow(h,2)/(3*pow(r,2)))
term5 = pow(term4,3)
f = term1+ term2 + term3 -term5
return f
# define first derivative of f(p)
def fprime(p):
term1 = 3*pow(p,2)*(1-pow(β,2))
term2 = 2*p*((0.4*h*pow(β,2))-(σ*pow(h,2)/pow(r,2)))
term3 = (σ**2*pow(h,4)/(3*pow(r,4)))
fp = term1 + term2+ term3
return fp
tol = 1e-4 # tolerance
iter = 1
p_guess = 20 # initial gues value of p
α = 1 # Relaxation factor
Preasure = []
iteration = []
for h in height:
while(abs(f(p_guess)) > tol):
p_guess = p_guess -α*(f(p_guess)/fprime(p_guess))
iter = iter +1
Preasure.append(p_guess)
iteration.append(iter)
# TO tabulate data in table
for i in range(len(height)):
if i == 0:
print("height : Preasure")
headers = ['iteration', 'Thichness' ,'Preasure']
table =tuple(zip(iteration , height, Preasure))
print(tabulate(table, headers= headers , tablefmt ='grid'))
# To plot the graph
plt.plot(height,Preasure)
plt.xlabel('Thichness of ice field')
plt.ylabel('Preasure in psf')
plt.show()
Output:
1.Thickness of ice field vs preasure plot.
fig: Preasure vs Thichness
2. Tabulate the data of thickness and Preasure.
Result:
From above graph and table ,we can see that as thickness of ice field increeses cushoin preasure increases non-linearly(Quadraticaly).
Case-2: Find optimum relaxation Factor
""" case 2
for h = 0.6 & find the optimum Relaxation factor """
import matplotlib.pyplot as plt
import math
from tabulate import tabulate
import numpy as np
#Define function
def f(p):
term1 = pow(p,3)*(1-pow(β,2))
term2 = pow(p,2)*((0.4*h*pow(β,2))-(σ*pow(h,2)/pow(r,2)))
term3 = p*(σ**2*pow(h,4)/(3*pow(r,4)))
term4 = (σ*pow(h,2)/(3*pow(r,2)))
term5 = pow(term4,3)
f = term1+ term2 + term3 -term5
return f
# define first derivative of f(p)
def fprime(p):
term1 = 3*pow(p,2)*(1-pow(β,2))
term2 = 2*p*((0.4*h*pow(β,2))-(σ*pow(h,2)/pow(r,2)))
term3 = (pow(σ,2)*pow(h,4)/(3*pow(r,4)))
fp = term1 + term2+ term3
return fp
# Input
β = 0.5 # width of ice wedge
r = 40 # size of air cushion
σ = 150*144 # tensile strength of ice
h = 0.6 # thichness of ice field
a =[]
iteration=[]
tol = 1e-4 # tolerance
alpha = np.arange(0.2,2,0.005)
for α in alpha:
iter = 1
p_guess = 10 # initial gues value of p
while(abs(f(p_guess)) > tol):
p_guess = p_guess -α*(f(p_guess)/fprime(p_guess))
iter = iter +1
a.append(α)
iteration.append(iter)
print('minimum iteration:',min(iteration))
headers = ['alpha','iteration']
table = tuple(zip(a,iteration))
print(tabulate(table,headers =headers,tablefmt = 'grid'))
# Ploting
plt.plot(a,iteration)
plt.xlabel('Relaxation Factor(α)')
plt.ylabel('No.of iterations')
plt.savefig('alpha vs iteration')
plt.show()
Output:
1.Relaxation factor vs No. of iteration plot.
fig: Relaxation factor vs No. of iteration plot
2.Table showing optimum relaxation factor for minimum no. of iteration
Result:
Optimum value of Relaxation factor is 1.055 and minimum no.of iteration are 7.
Case-3: Explain your results with respect to the appropriate relaxation factor.
""" Case3 : with optimum value of relaxation factor """
import matplotlib.pyplot as plt
import math
from tabulate import tabulate
import numpy as np
# Input
β = 0.5 # related to the width of the ice wedge
r = 40 # the size of the air cushion
σ = 150*144 # the tensile strength of the ice
height = [0.6,1.2,1.8,2.4,3.0,3.6,4.2]
# Define function
def f(p):
term1 = pow(p,3)*(1-pow(β,2))
term2 = pow(p,2)*((0.4*h*pow(β,2))-((σ*pow(h,2))/pow(r,2)))
term3 = p*(pow(σ,2)*pow(h,4)/(3*pow(r,4)))
term4 = (σ*pow(h,2)/(3*pow(r,2)))
term5 = pow(term4,3)
f = term1+ term2 + term3 -term5
return f
# define first derivative of f(p)
def fprime(p):
term1 = 3*pow(p,2)*(1-pow(β,2))
term2 = 2*p*((0.4*h*pow(β,2))-(σ*pow(h,2)/pow(r,2)))
term3 = (σ**2*pow(h,4)/(3*pow(r,4)))
fp = term1 + term2+ term3
return fp
tol = 1e-4 # tolerance
iter = 1
p_guess = 20 # initial gues value of p
α = 1.005 # Relaxation factor
Preasure = []
iteration = []
for h in height:
while(abs(f(p_guess)) > tol):
p_guess = p_guess -α*(f(p_guess)/fprime(p_guess))
iter = iter +1
Preasure.append(p_guess)
iteration.append(iter)
# TO tabulate data in table
for i in range(len(height)):
if i == 0:
print("height : Preasure")
headers = ['iteration', 'Thichness' ,'Preasure']
table =tuple(zip(iteration , height, Preasure))
print(tabulate(table, headers= headers , tablefmt ='grid'))
# To plot the graph
plt.plot(height,Preasure)
plt.xlabel('Thichness of ice field')
plt.ylabel('Preasure in psf')
plt.show()
Output:
1.Thickness of ice field vs preasure plot
fig: Preasure vs Thichness
2. Tabulate the data of thickness and Preasure.
Result:
From above graph and table ,we can see that as thickness of ice field increeses cushoin preasure increases non-linearly(Quadraticaly).
Explanation of code:
1.Define function(f(p_guess)) and first derivative of function(fprime(p_guess))
use 'return' to recall them outside the function.
2.To get accuracy upto four decimal point use 'tol = 1e-4', use 'np.arange(0.2,2,0.005)'
to increase the value of alpha from 0.2 to 2 by interval of 0.005
3.In order to consider only positive values of f(p_guess) use
'while(abs(f(p_guess)) > tol):'.
4.Use 'append' command to store the values of p_guess and iter valus into empty preasure and iteration matrix.
5. To tabulate the data call tabulate module from tabulate define heading of each coloumn using 'Headers' and to sort the data into coloumn use
'table = tuple(zip(a,iteration))'
6.To print the table use 'print(tabulate(table,headers =headers,tablefmt = 'grid'))'
Error occured:
1.Getting horizontal line in alpha vs no. of iteration plot.
fig: Error occured
Steps taken to overcome the error:
Enter the values of 'p_guess' and 'iter' inside the for loop.
References:
1.Newton-Raphson Method [https://brilliant.org/wiki/newton-raphson-method]
2.To create table in python [https://pypi.org/project/tabulate]