# Untitled

unknown
plain_text
7 months ago
2.4 kB
2
Indexable
Never
```# -*- 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.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()