Завдання 2

import numpy as np
import pandas as pd
import matplotlib.gridspec as gridspec
import matplotlib.pyplot as plt

def dx_dt(x, y):
    return x * ( eps_x - gamma_x * y )

def dy_dt(x, y):
    return - y * ( eps_y - gamma_y * x )


def runge_kute_system(x, y, h, a, b):
    X = [x]
    Y = [y]
    X_Y = []

    # print(f'h = {h}')
    t = a
    # print(f't = {t}  |  x = {x}  |  y = {y}')
    while round(t, 3) <= b:
        t = t + h

        k1 = h * dx_dt(x, y)
        k2 = h * dx_dt(x + h / 2, y + k1 / 2)
        k3 = h * dx_dt(x + h / 2, y + k2 / 2)
        k4 = h * dx_dt(x + h, y + k3)
        dx = (k1 + 2 * k2 + 2 * k3 + k4) / 6
        x += dx
        X.append(x)

        r1 = h * dy_dt(x, y)
        r2 = h * dy_dt(x + h / 2, y + k1 / 2)
        r3 = h * dy_dt(x + h / 2, y + k2 / 2)
        r4 = h * dy_dt(x + h, y + k3)
        dy = (r1 + 2 * r2 + 2 * r3 + r4) / 6
        y += dy
        Y.append(y)

        # print(f't = {t:.3f}  |  x = {x:.3f}  |  y = {y:.3f}')

    X_Y.append(X)
    X_Y.append(Y)

    # print('\n------------------------------------------------------------\n')
    return X_Y


def grafik(k=3):
    plt.title(f"gamma_x={gamma_x}; gamma_y={gamma_y}; e_x={eps_x}; e_y={eps_y}")
    plt.xlabel("x(t)")
    plt.ylabel("y(t)")
    plt.grid()
    lines = []
    labels = []
    for i in dani:
        x_y = runge_kute_system(i[0], i[1], i[2], i[3], i[4])
        line, = plt.plot(x_y[0], x_y[1])
        lines.append(line)
        labels.append(f'x_0={i[0]}; y_0={i[1]}; h={i[2]}')
    plt.legend(lines, labels, loc='center left', bbox_to_anchor=(1, 0.5))
    plt.show()

    if k == 3:
        plt.title(f"x(t) (gamma_x={gamma_x}; gamma_y={gamma_y}; e_x={eps_x}; e_y={eps_y})")
        plt.xlabel("t")
        plt.ylabel("x")
        plt.grid()
        lines = []
        labels = []
        for i in dani:
            x_y = runge_kute_system(i[0], i[1], i[2], i[3], i[4])
            plt.plot(np.arange(a, b+2*h, h), x_y[0], label=f'x(t); x_0={i[0]}; y_0={i[1]}; h={i[2]}')
        plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))
        plt.show()

        plt.title(f"y(t) (gamma_x={gamma_x}; gamma_y={gamma_y}; e_x={eps_x}; e_y={eps_y})")
        plt.xlabel("t")
        plt.ylabel("y")
        plt.grid()
        lines = []
        labels = []
        for i in dani:
            x_y = runge_kute_system(i[0], i[1], i[2], i[3], i[4])
            plt.plot(np.arange(a, b + 2 * h, h), x_y[1], label=f'y(t); x_0={i[0]}; y_0={i[1]}; h={i[2]}')
        plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))
        plt.show()


def info(eps_x, eps_y, gamma_x, gamma_y):

    A = np.array([[0, - (eps_y * gamma_y)/gamma_x],
                  [(eps_x * gamma_x)/gamma_y, 0]])
    B = np.array([eps_y/gamma_x, eps_x/gamma_y])
    stat_point = np.array([eps_y/gamma_x, eps_x/gamma_y])
    print(f"Стаціонарна точка: ({stat_point[0]}, {stat_point[1]})")
    eigenvalues, eigenvectors = np.linalg.eig(A)
    print("lambda: ", eigenvalues)

    if eigenvalues[0] != eigenvalues[1] and np.iscomplexobj(eigenvalues[0]) == False:

        if (eigenvalues[0] < 0) and (eigenvalues[1] < 0):
            print(f'Харкатер точки спокою: стійкий вузол ; Стійкість точки спокою: асимптотично стійка.')

        if (eigenvalues[0] > 0) and (eigenvalues[1] > 0):
            print(f'Харкатер точки спокою: нестійкий вузол ; Стійкість точки спокою: нестійка.')

        if (eigenvalues[0] > 0 and eigenvalues[1] < 0) or (eigenvalues[0] < 0 and eigenvalues[1] > 0):
            print(f'Харкатер точки спокою: сідло ; Стійкість точки спокою: нестійка.')

    elif np.iscomplexobj(eigenvalues[0]) == True:

        if (eigenvalues[0] < 0) and (eigenvalues[1] < 0):
            print(f'Харкатер точки спокою: стійкий фокус ; Стійкість точки спокою: асимптотично стійка.')

        if (eigenvalues[0] > 0) and (eigenvalues[1] > 0):
            print(f'Харкатер точки спокою: нестійкий фокус ; Стійкість точки спокою: нестійка.')

        if (eigenvalues[0] > 0 and eigenvalues[1] < 0) or (eigenvalues[0] < 0 and eigenvalues[1] > 0):
            print(f'Харкатер точки спокою: центр ; Стійкість точки спокою: стійка.')

    else:
        if (eigenvalues[0] < 0) and (eigenvalues[1] < 0):
            print(f'Харкатер точки спокою: стійкий вузол ; Стійкість точки спокою: асимптотично стійка.')

        else:
            print(f'Харкатер точки спокою: нестійкий вузол ; Стійкість точки спокою: нестійка.')
    return stat_point

eps_x = 4
gamma_x = 0.3
eps_y = 0.4
gamma_y = 0.4


h = 0.001
a = 0
b = 60

point = info(eps_x, eps_y, gamma_x, gamma_y)
dani = [[0.01, 2.0, h, a, b],
        [0.01, 0.01, h, a, b],
        [1.2, 1.0, h, a, b],
        [0.01, 0.5, h, a, b],
        [0.01, 1.0, h, a, b],
        [1.2, 0.6, h, a, b],
        [0.01, 1.5, h, a, b],
        [0.01, 0.2, h, a, b],
        [0.0, 0.0, h, a, b]]



grafik()
eps_x = 1
point = info(eps_x, eps_y, gamma_x, gamma_y)

grafik()