import numpy as np
import matplotlib.pyplot as plt
from math import pi
from math import cos, sin, exp
from math import log


def func(x):
    u = 0.5 -2*pi*x+sin(2*pi*x)
    return u


def B_k(k):# ДЛЯ ТОЧНОГО РЕШЕНИЯ ИЩУ КОЭФФИЦИЕНТЫ B_k
    A = 4 * sin(pi * k) / (2 * k - 1)
    B = 128 * (cos(pi * k)) / ((pi * pow((1 - 2 * k),2) * (4 * pow(k, 2) - 4 * k - 15)))

    return 2 * (A - B)


def exactsolution(X, t): # НАХОЖДЕНИЕ ТОЧНОГО РЕШЕНИЯ, РЯД ДО 1000, X - начальная сетка, t-время
    u = np.zeros(len(X))
    for j in range(len(u)):
        sum = 0
        for i in range(1, 40):
            sum = sum + B_k(i) * exp(-t * pow(pi * (2 * i - 1) / 2, 2)) * sin(pi * X[j] * (2 * i - 1) / 2)
        u[j] = sum + 0.5
    return u



def exactsolutioncheck(X, t): # НАХОЖДЕНИЕ ТОЧНОГО РЕШЕНИЯ, РЯД ДО 1000, X - начальная сетка, t-время
    sum = 0
    for i in range(1, 40):
        sum = sum + B_k(i) * exp(-t * pow(pi * (2 * i - 1) / 2, 2)) * cos(pi * X * (2 * i - 1) / 2)*(pi*(2*i-1)/2)
        u = sum + 0.5
    return u

def progonka(G, x):
    h = 1. / (len(x) - 1) # шаг по сетке
    tau = h/100
    A = (1. / (12. * tau)) - (1. / (2. * h * h))
    B = (5. / (6. * tau)) + (1. / (h * h))
    C = A
    P = np.zeros(len(x)-1)
    Q = np.zeros(len(x)-1)
    u = np.zeros(len(x))
    P[0] = 0.
    Q[0] = 1/2
    for j in range(1, len(x) - 1):
        P[j] = -C / (A * P[j - 1] + B) # пользуюсь учебником Михайлова, с.66
        Q[j] = (G[j] - A * Q[j - 1]) / (A * P[j - 1] + B)
    k = len(x) - 1
    u[k] = (4. * Q[k - 1] - P[k - 2] * Q[k - 1] - Q[k - 2]) / (3. - 4. * P[k - 1] + P[k - 2] * P[k - 1])
    while (k != 0):
        u[k - 1] = P[k - 1] * u[k] + Q[k - 1]
        k = k - 1
    return u

def Find_max_and_estab_time(e, t, x):#
    h = 1. / (len(x) - 1)
    tau = h/100
    firstcoeff = (1. / (12. * tau)) + (1. / (2. * h * h))
    secondcoeff = (5. / (6. * tau)) - (1. / (h * h))
    u_prev = np.zeros(len(x))
    G0 = np.zeros(len(x))
    for j in range(len(x)):
        u_prev[j] = func(x[j])
    u_next = u_prev
    L=0
    maX=0.
    if (e==0):
        while t > L*tau:
            L=L+1
            u_prev = u_next
            for j in range(1, len(x) - 1):
                G0[j] = u_prev[j + 1] * firstcoeff + u_prev[j] * secondcoeff + u_prev[j - 1] * firstcoeff
            u_next = progonka(G0, x)
            u_exact = exactsolution(x,L*tau)
            if max(np.abs(u_next - u_exact))>maX:
                maX = max(np.abs(u_next- u_exact))
        return maX
    if (e==1):
        for j in range(1, len(x) - 1):
            G0[j] = u_prev[j + 1] * firstcoeff + u_prev[j] * secondcoeff + u_prev[j - 1] * firstcoeff
        u_next = progonka(G0, x)
        while max(np.abs(u_next - u_prev)) > epsilon:
            L = L + 1
            u_prev = u_next
            for j in range(1, len(x) - 1):
                G0[j] = u_prev[j + 1] * firstcoeff + u_prev[j] * secondcoeff + u_prev[j - 1] * firstcoeff
            u_next = progonka(G0, x)
        return (L*tau)


def Graphics(t, x):# попробую вернуть максимум
    #вывод значения на k-м слое для tau = h*h/4
    h = 1. / (len(x) - 1.)
    tau=h/100
    firstcoeff = (1. / (12. * tau)) + (1. / (2. * h * h))
    secondcoeff = (5. / (6. * tau)) - (1. / (h * h)) # коэффициенты перед u[j-1],u[j],u[j+1] для нижнего слоя
    u_prev = np.zeros(len(x))
    for j in range(len(x)):
        u_prev[j] = 0.5 - 2. * pi * x[j] + sin(2. * pi * x[j])
    u_next=u_prev
    G0 = np.zeros(len(x))
    L=0
    while ((t/tau)-L) >0 :# получаю значения на нужном слое
        #print(t/tau)
        u_prev = u_next
        for j in range(1, len(x) - 1):
            G0[j] = u_prev[j + 1] * firstcoeff + u_prev[j] * secondcoeff + u_prev[j - 1] * firstcoeff
        u_next = progonka(G0, x)
        L=L+1
    return u_next

A = np.zeros
epsilon = 0.0000001
x10 = np.linspace(0., 1., 11)
x20 = np.linspace(0., 1., 21)
x15 = np.linspace(0.,1.,16)
x30 = np.linspace(0.,1.,31)
x25 = np.linspace(0., 1., 26)
x40 = np.linspace(0.,1.,41)
x80 = np.linspace(0.,1.,81)
e_t2 =Find_max_and_estab_time(1, 10, x10)
e_t3 = Find_max_and_estab_time(1, 10, x15)
e_t4 = Find_max_and_estab_time(1, 10, x20)
e_t5 = Find_max_and_estab_time(1, 10, x30)
e_t6 = Find_max_and_estab_time(1, 10, x40)
e_t7 = Find_max_and_estab_time(1, 10, x80)
#print(e_t2, "Время установления для 10 узлов")
#print(e_t3, "Время установления для 15 узлов")
#print(e_t4, "Время установления для 20 узлов")
#print(e_t5, "Время установления для 30 узлов")
#print(e_t6, "Время установления для 40 узлов")
#print(e_t7, "Время установления для 80 узлов")
n10=Find_max_and_estab_time(0, e_t2, x10)
n15=Find_max_and_estab_time(0, e_t3, x15)
n20=Find_max_and_estab_time(0, e_t4, x20)
n30=Find_max_and_estab_time(0, e_t5, x30)
n40=Find_max_and_estab_time(0, e_t6, x40)
n80 = Find_max_and_estab_time(0, e_t7, x80)
print(log(n10/n20,2.),"Порядок сходимости при N=10")
print(log(n15/n30,2),"Порядок сходимости при N=15")
print(log(n20/n40,2.),"Порядок сходимости при N=20")
print(log(n40/n80,2),"Порядок сходимости при N=40\n")
print(n10,"Погрешность на сетке из 10 узлов")
print(n15,"Погрешность на сетке из 15 узлов")
print(n20,"Погрешность на сетке из 20 узлов")
print(n30,"Погрешность на сетке из 30 узлов")
print(n40,"Погрешность на сетке из 40 узлов")
plt.ioff()
fig, (ax1, ax2,ax3) = plt.subplots(nrows=1, ncols=3, figsize=(19,5))
t1=0
t2=1
t3=2
ax1.set_title("t=0")
ax2.set_title("t=1")
ax3.set_title("t=2")
ax1.plot(x80,exactsolution(x80,t1))
ax1.plot(x10, Graphics(t1,x10))
ax1.plot(x15, Graphics(t1,x15))
ax1.plot(x20, Graphics(t1,x20))
ax1.plot(x30, Graphics(t1,x30))
ax1.plot(x40, Graphics(t1,x40))
ax1.plot(x80, Graphics(t1,x80))
ax2.plot(x80,exactsolution(x80,t2))
ax2.plot(x10, Graphics(t2,x10))
ax2.plot(x15, Graphics(t2,x15))
ax2.plot(x20, Graphics(t2,x20))
ax2.plot(x30, Graphics(t2,x30))
ax2.plot(x40, Graphics(t2,x40))
ax2.plot(x80, Graphics(t2,x80))
ax3.plot(x80,exactsolution(x80,t3))
ax3.plot(x10, Graphics(t3,x10))
ax3.plot(x15, Graphics(t3,x15))
ax3.plot(x20, Graphics(t3,x20))
ax3.plot(x30, Graphics(t3,x30))
ax3.plot(x40, Graphics(t3,x40))
ax3.plot(x80, Graphics(t3,x80))
ax1.grid()
ax2.grid()
ax3.grid()
lgnd = ax1.legend(['Точное решение','N=10','N=15','N=20','N=30','N=40','N=80'],)
lgnd = ax2.legend(['Точное решение','N=10','N=15','N=20','N=30','N=40','N=80'],)
lgnd = ax2.legend(['Точное решение','N=10','N=15','N=20','N=30','N=40','N=80'],)

plt.show()