Demo entry 6634751

NewtonRaphsonR1

   

Submitted by anonymous on Aug 14, 2017 at 01:40
Language: Python console session. Code size: 2.0 kB.

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Sun Aug 13 15:11:19 2017

@author: Pamela
"""
import matplotlib
import numpy as np
import matplotlib.pyplot as plt

#Newton's Method
def Newton__Solver(r,H,dfdx,tol,nitmax,x_start):
    error = tol + 1
    nit = 1
    xn = np.zeros((2,2))
    xn[0] = x_start
    while (nit<nitmax) and (error>tol):
       gn = dfdx(xn[0]) #Evaluate (df/dx or grad(f))
       delta_xn = np.empty((1,2))
       h = H(xn[0])
       hh =([[h[0], h[1]],[h[2],h[3]]]) #Change size of hessian matrix
       det =  np.linalg.det(hh) #Calculate determinat of hessian matrix
       #if (det != 0):
       delta_xn = -np.linalg.solve(hh,gn) # dx = -inv(H) * grad
       xn[1] += xn[0]+delta_xn
       error = np.abs(delta_xn)
       print ("Iteration:",nit,"Approximate Root:",xn[1],"Error:",error)
       plt.plot(xn[:,0],xn[:,1],'r-o')
       plt.plot(x_start[0],x_start[1],'k-o')
       plt.text(x_start[0],x_start[1], r'$(2, 2)$')
       nit += 1
       #else:
          #print ("Does not converge")
       return xn
   
# Define Function
def f(x):
    x1 = x[0]
    x2 = x[1]
    obj = np.sin(x1)*np.sin(x2)
    return obj
# Define Gradient (first derivatives)
def dfdx(x):
    x1 = x[0]
    x2 = x[1]
    grad = []
    grad.append(np.cos(x1)*np.sin(x2))
    grad.append(np.sin(x1)*np.cos(x2))
    return grad
# Define Hessian (second derivatives)
def H(x):
    x1 = x[0]
    x2 = x[1]
    hessian = []
    hessian.append(-np.sin(x1)*np.sin(x2))
    hessian.append(np.cos(x1)*np.cos(x2))
    hessian.append(-np.cos(x1)*np.cos(x2))
    hessian.append(-np.sin(x1)*np.sin(x2))
    return hessian
# Initial parameters
x_start = [2.0, 2.0]
tol = 1e-12
nitmax = 12

# Create contour plot
i1 = np.arange(-4.0, 4.0, 0.1)
i2 = np.arange(-4.0, 4.0, 0.1)
x1_mesh, x2_mesh = np.meshgrid(i1, i2)
f_mesh = np.sin(x1_mesh)*np.sin(x2_mesh)
plt.contour(x1_mesh, x2_mesh, f_mesh)
#plt.title('r(x,y) = sin(x)*sin(y)')
plt.xlabel('x')
plt.ylabel('y')

xn = Newton__Solver(f,H,dfdx,tol,nitmax,x_start)

This snippet took 0.00 seconds to highlight.

Back to the Entry List or Home.

Delete this entry (admin only).