<a href="https://colab.research.google.com/github/maxigaarp/teaching/blob/master/optimization/Cap11-Metodo%20de%20Newton-Raphson.ipynb" target="_parent"></a>


```python
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from sympy import diff, Symbol
import numpy as np

"""
Basic version of Newton-Raphson Algorithm with symbolic math.

Just define the number of variables, the function f(z1,z2,...,zn) and the algorithm 
does the calculation. 

Issues: Convergence is not always granted. With more than 2 variables the plot of the 
third cell is not available.
"""

n = 3 # numero de variables
z = [Symbol('z'+str(i+1)) for i in range(n)] # definicion de variables z1,z2,...,zn

f = 0.6*z[0]**4+2*z[1]**4+3*z[0]**3- 2*z[1]**2 +30*z[0] +3 #definicion de funcion f(z)
f = (z[0]-2)**2 + (z[1]-1)**2 -z[2]*(z[0]+z[1]-1)#definicion de funcion f(z)

#f = z[0]**2 + 2*z[1]**2 #definicion de funcion f(z)

```


```python
np.set_printoptions(formatter={'float': lambda x: "{0:0.3f}".format(x)})



gradient=[f.diff(zi) for zi in z] # calculo del gradiente f'(z)
hessian=[[f.diff(zi).diff(zj) for zi in z] for zj in z]



x=[np.array([0.5,0.5,-1])] #iteraciones del metodo del gradiente, valor inicial x[0]=(2,1)

#t=0.1 #valor del paso fijo
maxiter=1000 #numero maximo de iteraciones

print("Initial Point: "+str(x[0]))

for k in range(maxiter):
    nablaf_xk=np.array([nf.evalf(subs = dict(zip(z,x[k]))) for nf in gradient],dtype="float32") #evaluar x[k] en gradiente
    hessianf_xk=np.matrix([[hij.evalf(subs = dict(zip(z,x[k]))) for hij in hi] for hi in hessian],dtype="float32") #evaluar x[k] en matriz hessiana

    x.append(x[k]-np.tensordot(np.linalg.inv(hessianf_xk),nablaf_xk,1))
    
    if np.linalg.norm(nablaf_xk)<0.00001: # termino del algoritmo
        print("Optimo: "+str(x[k+1]))
        break
    else:
        print("Iteration "+str(k+1)+", gradient: "+str(nablaf_xk)+", x: "+str(x[k+1]))
```

    Initial Point: [0.500 0.500 -1.000]
    Iteration 1, gradient: [-2.000 0.000 0.000], x: [1.000 0.000 -2.000]
    Optimo: [1.000 0.000 -2.000]



```python
import plotly.graph_objects as go

xd=np.linspace(-10,5,101)
yd=np.linspace(-5,5,101)
X,Y=np.meshgrid(xd,yd)
Z=np.array([[f.evalf(subs = dict(zip(z,[x,y]))) for x in xd] for y in yd]).astype(float)
fig = go.Figure()
fig.add_trace(go.Surface(z=Z, x=X, y=Y, opacity=0.3))

xp=np.round([p[0] for p in x])
yp=np.round([p[1] for p in x])
zp=np.array([f.evalf(subs = dict(zip(z,[p[0],p[1]]))) for p in x]).astype(float)

fig.add_trace(go.Scatter3d(x=xp, y=yp, z=zp,marker=dict(size=4,color="red"),
                                line=dict(color='darkblue',width=2)))

fig.update_layout(title='Newton-Raphson Optimization', autosize=True,
                  width=1000, height=1000,
                  margin=dict(l=65, r=50, b=65, t=90))

fig.show()

```


<html>
<head><meta charset="utf-8" /></head>
<body>
    <div>
            
                
            
            <div id="33fe5040-c915-4a33-82b7-3d51101b9959" class="plotly-graph-div" style="height:1000px; width:1000px;"></div>
            
        </div>
</body>
</html>


ipynb
e858271f8931a6fe14c635a466dfd6c9dad3863d
Jupyter Notebook
2f6f7aa448bab44553e88f26a42984b487d1d8a4
MIT
maxigaarp/teaching
optimization/Cap11-Metodo de Newton-Raphson.ipynb
2f6f7aa448bab44553e88f26a42984b487d1d8a4
MIT
maxigaarp/teaching
optimization/Cap11-Metodo de Newton-Raphson.ipynb
2f6f7aa448bab44553e88f26a42984b487d1d8a4
MIT
maxigaarp/teaching
optimization/Cap11-Metodo de Newton-Raphson.ipynb