lab2

import numpy as np
import math
import pandas as pd

eps=0.1
h=0.001

def f(x,y):
    return x**2+8*x*y+0.001*x*y-x-y

def dFx(x,y):
    return (f(x+h,y)-f(x-h,y))/(2*h)

def dFy(x,y):
    return (f(x,y+h)-f(x,y-h))/(2*h)

def d2Fx(x,y):
    return (f(x+h,y)-2*f(x,y)+f(x-h,y))/(h**2)

def d2Fy(x,y):
    return (f(x,y+h)-2*f(x,y)+f(x,y-h))/(h**2)

def d2Fxy(x,y):
    return (f(x+h,y+h)-f(x+h, y-h)-f(x-h, y+h)-f(x-h,y-h))/(4*(h**2))

def gradF(x,y):
    return np.array([dFx(x,y), dFy(x,y)], float)

def sdevF(x,y):
    return  np.array([[d2Fx(x,y), d2Fxy(x,y)],[d2Fxy(x,y),d2Fy(x,y)]], float)

def Newton(x,y):
    xy = np.array([x,y], float)
    alpha = 1.0
    while True:
        hh = (-1) * np.dot(np.linalg.inv(sdevF(x,y)), gradF(x,y))
        xyk = xy + alpha*hh
        if f(x,y)-f(xyk[0],xyk[1]) <= eps * alpha * np.dot(gradF(xyk[0],xyk[1]),hh):
            break
        else:
            alpha = alpha * 0.5
            xy = np.copy(xyk)
    return  xyk

Newton(0,0)