Untitled

# -*- coding: utf-8 -*-
"""
Created on Thu Oct 19 11:25:12 2023

Newton-Rhapson 2D
+ 26.10.2023



@author: 48609
"""

import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import solve, norm
from matplotlib import cm

plt.close('all')

#A = 5
#mx = 1e-1
#my = 5e-2
#def f1(x,y):
 #   return A*np.exp((x*mx)**2+(y*my)**2)


def f1(x,y):
    return 2*np.sin(x) + 3*np.cos(y)**2

def f2(x,y):
    return x**2 + 2*y**2 - 5

# vector function:
# vector xy contains x (as xy[0]) and y (as xy[1])
def f(xy):
    return np.array([f1(xy[0],xy[1]), f2(xy[0],xy[1])])

def FI(xy):
    return np.array([[2*np.cos(xy[0]) , -6*np.cos(xy[1])*np.sin(xy[1])], [2*xy[0], 4*xy[1]]])

xy0 = np.array([0.5,2])
xy = [xy0]

more = True
cond_suc = 1e-9    # condition on success (finding the zero point)
max_count = 1000    # maximum number of loops
count = 0

while more:
    count += 1
    correction = solve(FI(xy[-1]), f(xy[-1]))
    xy.append(xy[-1] - correction)
    result = norm(f(xy[-1]))
    if result < cond_suc:
        print("Zero point (f1={}, f2={}) found at point x={}, y={} after {} loops".format(f(xy[-1])[0],f(xy[-1])[1],xy[-1][0],xy[-1][-1], count))
        more= False
    if cond_suc > max_count:
        print("Too many loops (={}), last point: (x={}, y={}) with value {}".format(count, xy[-1][0], xy[-1][1],f(xy[-1])))
    
xy1 = xy0 - correction

correction = solve(FI(xy1), f(xy1))
xy2 = xy1 - correction

#plots

xs = -3
xe = 3
Nx = 500
xp = np.linspace(xs,xe,Nx)
ys = -2
ye = 4
Ny = 500
yp = np.linspace(ys,ye,Ny)
[Xp,Yp] = np.meshgrid(xp,yp)

Z1 = f1(Xp,Yp)
Z2 = f2(Xp,Yp)

#fig = plt.figure()
#ax = fig.add_subplot(111,projection='3d')
#ax.set(xlabel='x',ylabel='y',zlabel='Z1',title='function f1')
#ax.plot_surface(Xp,Yp,Z1,cmap=cm.jet)


#ax.plot_surface(Xp,Yp,Z2,cmap=cm.hot)

fig = plt.figure()
ax = fig.add_subplot(111)
ax.contour(Xp,Yp,Z1,[0],colors='b')
ax.contour(Xp,Yp,Z2,[0],colors='g')
ax.set(xlabel='x',ylabel='y',title='function f1 (blue) and f2 (green)')
ax.axis('equal')

xy = np.array(xy).T

ax.plot(*xy0, 'or')     # *  - means use two components
ax.plot(xy[0],xy[1], 'x-r')
#ax.plot(*xy2, 'xr')
ax.plot(*xy[:,-1], 'xb')

ax.set_title('function f1 (blue) and f2 (green)\nZero point (f1={:.3e}, f2={:.3e})\nfound at point x={:.3f}, y={:.3f} after {} loops'.format(f(xy[:,-1])[0],f(xy[:,-1])[1],xy[0,-1],xy[1,-1], count))