Optimization of Ice Pressure using Python

import numpy
import scipy
import matplotlib.pyplot as pyplot



def f(p,h):
	t1 = pow(p,3)*(1 - pow(b,2))
	t2 = pow(p,2)*((0.4*h*pow(b,2) - sigma*pow(h,2)/pow(r,2)))
	t3 = (pow(sigma,2)*pow(h,4)/(3*pow(r,4)))*p
	t4 = sigma*pow(h,2)/(3*pow(r,2))
	return t1+t2+t3+t4

def df(p,h):
	dt1 = 3*pow(p,2)*(1 - pow(b,2))
	dt2 = 2*p*((0.4*h*pow(b,2) - sigma*pow(h,2)/pow(r,2)))
	dt3 = (pow(sigma,2)*pow(h,4)/(3*pow(r,4)))
	return dt1+dt2+dt3

b = 0.5   
r = 40
sigma = 150*12*12
height = [0.6,1.2,1.8,2.4,3.0,3.6,4.2]
t = 0.6
alpha = 0.5
error = 1e-3
p_guess = 500
pressure = []
iter = 1

for h in height:
	p_guess = 500
	while abs(f(p_guess,h)>error):
		p_guess = p_guess - (f(p_guess,h)/df(p_guess,h))
		iter = iter + 1

	pressure.append(p_guess)


print(pressure)
print('No of Iteration',iter)
pyplot.plot(height,pressure)
pyplot.xlabel('Height')
pyplot.ylabel('pressure')
pyplot.show()

Output

[10.623201749071786, 9.018651268466796, 4.36265759918318, 2.608408632174104, 255.79162125322364, 34.63972943108442, 277.142174983423]
No of Iteration 415