I am trying to write the basic Newton's method without pre-built solvers.
This is the function:
## definition of variables
x_1, x_2 = sym.symbols("x_1 x_2")
a_T=np.array([[0.3],[0.6],[0.2]])
b_T=np.array([5,26,3])
c_T=np.array([40,1,10])
u= x_1-0.8
v= x_2-(a_T[0]+a_T[1]*u**2*(1-u)**(1/2)-a_T[2]*u)
alpha= -b_T[0]+(b_T[1]*u**2)*(1+u)**(1/2)+(b_T[2])*u
beta= c_T[0]*v**2*(1-c_T[1]*v)/(1+c_T[2]*u**2)
## function
f = alpha**(-beta)
I calculated the gradient and the Hessian and defined the other parameters:
## gradient
gradient_cal = sym.Matrix(1,2,sym.derive_by_array(f, (x_1, x_2)))
## hessian
hessian_cal = sym.Matrix(2, 2, sym.derive_by_array(gradient_cal, (x_1, x_2)))
# initial guess
x_A= Matrix([[1],[0.5]])
xk = x_A
#tolerance
epsilon= 1e-10
#maximum iterations
max_iter=100
And the function itself:
def newton(gradient_cal,hessian_cal,xk,epsilon,max_iter):
for k in range(0,max_iter):
fxk = gradient_cal.evalf(subs={x_1:xk[0], x_2:xk[1]})
if fxk.norm() < epsilon:
print('Found solution after',k,'iterations.')
return xk
Dfxk = hessian_cal.evalf(subs={x_1: xk[0], x_2: xk[1]})
if Dfxk == 0:
print('Zero derivative. No solution found.')
return None
A=hessian_cal.evalf(subs={x_1: xk[0], x_2: xk[1]})
B=gradient_cal.evalf(subs={x_1: xk[0], x_2: xk[1]})
pk= (A.inv().dot(B))
xk = np.subtract(xk, pk)
print('Exceeded maximum iterations. No solution found.')
return None
approx = newton(gradient_cal,hessian_cal,x_A,epsilon,max_iter)
print(approx)
The following error shows up:
TypeError: Shape should contain integers only.
I checked it and saw that the Hessian contains "I" values. Therefore I am not sure if the calculations of the gradient and the Hessian are correct.
Does anyone have a better solution to calculate the gradient and the Hessian for such a complex function?
