import math
from matplotlib import pyplot as plot


class CauchySolver:
    def __init__(self, u10, u20, l1, l2, p1, p2, x0, xn):
        def _phi(j, x):
            return 0.5 * (u10 + u20) * math.exp(l1 * x) + (-1) ** (j + 1) * 0.5 * (u10 - u20) * math.exp(l2 * x)

        self.phi = _phi
        self.u = {1: u10, 2: u20}

        self.p1 = p1
        self.p2 = p2

        self.x0 = x0
        self.xn = xn

    def get_setka(self, n):
        h = (self.xn - self.x0) / n
        x = [self.x0 + i * h for i in range(n + 1)]
        n = n + 1
        return x, n, h

    def y_euler(self, n):
        x, n, h = self.get_setka(n)
        u1 = []
        u2 = []
        for i in range(n):
            if i == 0:
                u1.append(self.u[1])
                u2.append(self.u[2])
            else:
                u1.append(u1[i - 1] + h * (self.p1 * u1[i - 1] + self.p2 * u2[i - 1]))
                u2.append(u2[i - 1] + h * (self.p2 * u1[i - 1] + self.p1 * u2[i - 1]))
        return u1, u2

    def y_euler_implicit(self, n):
        x, n, h = self.get_setka(n)
        u1 = []
        u2 = []
        for i in range(n):
            if i == 0:
                u1.append(self.u[1])
                u2.append(self.u[2])
            else:
                f1 = self.p1 * u1[i - 1] + self.p2 * u2[i - 1]
                f2 = self.p2 * u1[i - 1] + self.p1 * u2[i - 1]

                u1i = u1[i - 1] + h * f1
                u2i = u2[i - 1] + h * f2

                u1.append(u1[i - 1] + 0.5 * h * (f1 + self.p1 * u1i + self.p2 * u2i))
                u2.append(u2[i - 1] + 0.5 * h * (f2 + self.p2 * u1i + self.p1 * u2i))

        return u1, u2

    def y_runge_cutte_b(self, n):
        x, n, h = self.get_setka(n)
        u1 = []
        u2 = []
        for i in range(n):
            if i == 0:
                u1.append(self.u[1])
                u2.append(self.u[2])
            else:
                k11 = h * (self.p1 * u1[i - 1] + self.p2 * u2[i - 1])
                k21 = h * (self.p2 * u1[i - 1] + self.p1 * u2[i - 1])

                k12 = h * (self.p1 * (u1[i - 1] + k11 / 2) + self.p2 * (u2[i - 1] + k21 / 2))
                k22 = h * (self.p2 * (u1[i - 1] + k11 / 2) + self.p1 * (u2[i - 1] + k21 / 2))

                k13 = h * (self.p1 * (u1[i - 1] + k12 / 2) + self.p2 * (u2[i - 1] + k22 / 2))
                k23 = h * (self.p2 * (u1[i - 1] + k12 / 2) + self.p1 * (u2[i - 1] + k22 / 2))

                k14 = h * (self.p1 * (u1[i - 1] + k13) + self.p2 * (u2[i - 1] + k23))
                k24 = h * (self.p2 * (u1[i - 1] + k13) + self.p1 * (u2[i - 1] + k23))

                u1.append(u1[i - 1] + (k11 + 2 * k12 + 2 * k13 + k14) / 6)
                u2.append(u2[i - 1] + (k21 + 2 * k22 + 2 * k23 + k24) / 6)

        return u1, u2

    def eps(self, u1, u2):
        x, n, h = self.get_setka(len(u1) - 1)
        eps = []

        u11, u22 = [self.phi(1, xi) for xi in x], [self.phi(2, xi) for xi in x]

        m = max(abs(abs(u1[i]) - abs(u11[i])) + abs(abs(u2[i]) - abs(u22[i])) for i in range(len(u1)))

        return m

    def plot(self, j, y):
        x, n, h = self.get_setka(len(y) - 1)
        phi_x = [self.x0 + i * (self.xn - self.x0) / 1000 for i in range(1000)]
        phi_y = [self.phi(j, xi) for xi in phi_x]

        plot.xscale('log')
        plot.yscale('log')
        plot.plot(phi_x, phi_y, color='red', label='phi')
        plot.plot(x, y, color='blue', label='y')
        plot.legend()
        plot.show()

    def ppp(self, e):
    	x = [i for i in range(len(e))]
    	#plot.xscale('log')
    	plot.yscale('log')
    	plot.plot(x, e)
    	plot.show()

def main():
    n = 10

    solver = CauchySolver(u10=0, u20=1, l1=-201, l2=-1, p1=-101, p2=-100, x0=0, xn=5)

    e = []

    for i in range(15):
        u1, u2 = solver.y_runge_cutte_b(n)
        e.append(solver.eps(u1,u2))
        # u1_1, u2_2 = solver.y_euler_implicit(n * 2)
        # max_loc_eps = solver.eps(u1, u2)
        # max_loc_eps = solver.eps(u1_1, u2_2)
        # print('погрешность %.2E' % (max_loc_eps))
        # pz = math.log (abs((solver.eps(u1,u2)) /(solver.eps(u1_1,u2_2))), 2)
        # print('порядок %.2E' % (pz))
        # if i != 0:
            # pz = max_loc_eps_old / max_loc_eps
            # print('погрешность %.8E' % (pz,))

        # max_loc_eps_old = max_loc_eps
        n = n * 2

    solver.ppp(e)

    # u1, u2 = solver.y_runge_cutte_b(n)
    # solver.plot(1, u1)

main()