from math import sqrt
import numpy as np
import matplotlib.pyplot as plt

# Точность
Epsilon = 1e-12

# Пусть c = 7.0, в этом случае, согласно теории, есть 2 корня
c = 7.0  # Параметр c, при котором есть 2 корня
X1 = (1379.0 - sqrt(19159.0))/196.0  # Первый корень (расчитан теоретически на листочке)
X2 = (1379.0 + sqrt(19159.0))/196.0  # Второй корень (расчитан теоретически на листочке)

# Наша функция, зависящая от X, после того, как мы подставили параметр c в систему, и подстановкой выразили Y через X
def Foo1(X):
    return 392.0*X**2 - 5516.0*X + 19209.0

# Первая производная от вышеуказанной функции
def Foo2(X):
    return 784.0*X - 5516.0

# Метод Ньютона
def NyutonMethod(X0, f, f1, Eps, ResRateArr, needed_X, solutions_exist):
    Xk = X0  # Искомое решние
    Xk_1 = Xk + 1.0  # Xk+1
    Count = 0  # Счётчик итераций
    first_iteration = True  # Флажок для первой итерации
    print("Расчёты итераций для Xk...")
    while abs(Xk_1 - Xk) >= Eps:
        if first_iteration:
            first_iteration = False
        else:
            Xk = Xk_1
        Xk_1 = Xk - f(Xk)/f1(Xk)  # Xk+1 = Xk - f(Xk) / f'(Xk)
        Count += 1  # Делаем инкремент итерации на счётчике
        nevyazka = abs(needed_X - Xk_1)
        ResRateArr.append(nevyazka)  # Запоминаем невязку на каждой итерации
        print(Count, ") Xk =", Xk_1, ", |X* - Xk| =", nevyazka)
        if Count == 30 and not solutions_exist:
            print("За 30 итераций метод Ньютона так и не смог найти никаких корней...")
            return Count
    print("Корень функции, полученный с помощью метода Ньютона с начальным приближением X0 =", X0, ", равен", Xk)
    return Count

print("\nПараметр c =", c, "(один из множества случаев, когда система имеет 2 корня)")
print("Вид функции при c =", c, ": f(x) = 392x^2 - 5516x + 19209, f'(x) = 784x - 5516")
print("Корни этой функции, расчитанные теоретически:", X1, ",", X2)

X_0 = 4.0
print("\nПоиск первого корня методом Ньютона с начальным приближением X0 =", X_0)
NM_ResRateArr_1 = []  # Массив ординат (норм невязки на разных итерациях) для графика Метода Ньютона для 1-ого корня
NM_IterAmount_1 = NyutonMethod(X_0, Foo1, Foo2, Epsilon, NM_ResRateArr_1, X1, True)  # Делаем расчёты и запоминаем количество итераций
NM_IterArr_1 = np.arange(1, NM_IterAmount_1 + 1)  # Массив абсцисс (количества итераций) для графика Метода Ньютона для 1-ого корня

X_0 = 10.0
print("\nПоиск второго корня методом Ньютона с начальным приближением X0 =", X_0)
NM_ResRateArr_2 = []  # Массив ординат (норм невязки на разных итерациях) для графика Метода Ньютона для 2-ого корня
NM_IterAmount_2 = NyutonMethod(X_0, Foo1, Foo2, Epsilon, NM_ResRateArr_2, X2, True)  # Делаем расчёты и запоминаем количество итераций
NM_IterArr_2 = np.arange(1, NM_IterAmount_2 + 1)  # Массив абсцисс (количества итераций) для графика Метода Ньютона для 2-ого корня

# Пусть c = 11sqrt(2)/2, в этом случае, согласно теории, есть толькой один корень (это один из четырёх таких случаев)
c = sqrt(121.0 / 2.0)  # Параметр c, при котором есть только один корень
X1 = 5*sqrt(2)  # Первый и единственый корень (расчитан теоретически на листочке)

# Наша функция, зависящая от X, после того, как мы подставили параметр c в систему, и подстановкой выразили Y через X
def Foo3(X):
    return X**2 - 10*sqrt(2)*X + 50

# Первая производная от вышеуказанной функции
def Foo4(X):
    return 2*X - 10*sqrt(2)

print("\nПараметр c =", c, "(один из четырёх случаев, когда система имеет только один корень)")
print("Вид функции при c =", c, ": f(x) = x^2 - 10sqrt(2) + 50, f'(x) = 2x - 10sqrt(2)")
print("Корень этой функции, расчитанный теоретически:", X1)

X_0 = 5.0
print("\nПоиск первого и единственного корня методом Ньютона с начальным приближением X0 =", X_0)
NM_ResRateArr_3 = []  # Массив ординат (норм невязки на разных итерациях) для графика Метода Ньютона для 1-ого корня
NM_IterAmount_3 = NyutonMethod(X_0, Foo3, Foo4, Epsilon, NM_ResRateArr_3, X1, True)  # Делаем расчёты и запоминаем количество итераций
NM_IterArr_3 = np.arange(1, NM_IterAmount_3 + 1)  # Массив абсцисс (количества итераций) для графика Метода Ньютона для 1-ого корня

# Пусть c = 10, в этом случае, согласно теории, корней вообщей нету (дискриминант отрицательный, пересечений нету)
c = 10.0  # Параметр c, при котором нет корней

# Наша функция, зависящая от X, после того, как мы подставили параметр c в систему, и подстановкой выразили Y через X
def Foo5(X):
    return 800*X**2 - 11960*X + 49401

# Первая производная от вышеуказанной функции
def Foo6(X):
    return 1600*X - 11960

print("\nПараметр c =", c, "(один из множества случаев, когда система вообще не имеет корней)")
print("Вид функции при c =", c, ": f(x) = 800x^2 - 11960x + 49401, f'(x) = 1600x - 11960")

X_0 = 5.0
print("\nПопытки найти какие-то корни Методом Ньютона с начальным приближением X0 =", X_0)
NM_ResRateArr_4 = []  # Массив ординат (норм невязки на разных итерациях) для графика Метода Ньютона
NM_IterAmount_4 = NyutonMethod(X_0, Foo5, Foo6, Epsilon, NM_ResRateArr_4, X1, False)  # Делаем расчёты и запоминаем количество итераций
NM_IterArr_4 = np.arange(1, NM_IterAmount_4 + 1)  # Массив абсцисс (количества итераций) для графика Метода Ньютона

# Строим графики сходимостей
plt.semilogy(NM_IterArr_1, NM_ResRateArr_1, label="NM1, c = 7.0 (1-ый корень из двух)")
plt.semilogy(NM_IterArr_2, NM_ResRateArr_2, label="NM2, c = 7.0 (2-ой корень из двух)")
plt.semilogy(NM_IterArr_3, NM_ResRateArr_3, label="NM3, c = 11sqrt(2)/2 (единственный корень)")
plt.semilogy(NM_IterArr_4, NM_ResRateArr_4, label="NM4, c = 10 (нет корней)")
plt.title("Диаграммы сходимости")
plt.xlabel("Номер итерации")
plt.ylabel("Норма невязки на этой итерации")
plt.legend()
plt.grid(True)
plt.show()