# 4-ая лаба 2-ой вариант
def u_0yt(x,y,t):
    return 0

def u_1yt(x,y,t):
    return sin(y)*sin(t)

def u_x0t(x,y,t):
    return 0

def uy_x1t(x,y,t):
    return -sin(x)*sin(t)

def u_xy0(x,y,t):
    return 0

def f_xyt(x,y,t):
    return sin(x)*sin(y)*(cos(t)+2*sin(t))

def main3(u_0yt,u_1yt,u_x0t,u_x1t,u_xy0,x0,x1,y0,y1,h,t0,t1,tau):
    In = int((x1-x0)/h)
    Jn = int((y1-y0)/h)
    Kn = int((t1-t0)/tau)
    tau_2 = tau/2
    r=tau/(2*h*h)

    I = np.linspace(x0,x1,In) # массив с x
    J = np.linspace(y0,y1,Jn) # массив с y
    K = np.linspace(t0,t1,Kn) # массив с t

    w=np.zeros((Jn,In,Kn))

    # 0yt уже заполнен нулями
    # 1yt
    for k in range(0,Kn):
        for j in range(0,Jn):
            w[j,In-1,k] = u_1yt(x1,J[j],K[k])

    # x0t уже заполнен нулями
    # xy0 уже заполнен нулями

    for k in range (1,Kn): # идем по t

        P=np.zeros((Jn,In)) # промежуточный слой

        for j in range (1,Jn-1): # идем по y

            A=np.zeros((In-2,In-2))
            b=np.zeros(In-2)

            for n in range (0,In-2):
                A[n,n] = -(2*r+1)
            for n in range (0,In-3):
                A[n+1,n] = r
                A[n,n+1] = r
            for i in range (0,In-2):
                b[i] = -r*w[j+1,i,k-1]+(2*r-1)*w[j,i,k-1]-r*w[j-1,i,k-1]-r*h*h*f_xyt(I[i],J[j],K[k]-tau_2)

            #b[0] -= r*u_0yt(x0,J[j],K[k]-tau_2) вернёт 0, поэтому бесполезное действие
            b[In-3] -= r*u_1yt(x1,J[j],K[k]-tau_2)

            p = solve(A, b)
            #print(p[0])
            for i in range(0,In-2):
                P[j,i+1] = p[i]

            #P[j,0] = u_0yt(I[0],J[j],K[k]-tau_2)
            P[j,In-1] = u_1yt(I[In-1],J[j],K[k]-tau_2)

        #краевое первое не требуется тк оно 0
        '''
        for i in range(In):
            P[0,i] = u_x0t(i,Jn-1,K[k]-tau_2)

        '''
        for i in range(In):
            P[Jn-1,i] = P[Jn-2,i]+h*uy_x1t(i,Jn-1,K[k]-tau_2)
        for i in range (1,In-1):   # идем по x

            A=np.zeros((Jn-1,Jn-1))
            b=np.zeros(Jn-1)

            for n in range (0,Jn-1):
                A[n,n]=(1+2*r)
            for n in range (0,Jn-2):
                A[n+1,n]=-r
                A[n,n+1]=-r

            for n in range (0,Jn-1):
                A[n,Jn-3] += 1
                A[n,Jn-2] -=  -1

            for j in range (0,Jn-1):
                b[j] = r*P[j,i+1]+(1-2*r)*P[j,i]+r*P[j,i-1]+r*h*h*f_xyt(I[i],J[j],K[k]-tau_2)

            b[Jn-2] -= h*uy_x1t(I[i],J[j],K[k])

            p = solve(A, b)

            for j in range(0,Jn-1):
                w[j,i,k] = p[j]

    return I, J, K, w

def plot_3d_time(x,y,t,W,j):
    X, Y = np.meshgrid(x, y)

    fig = plt.figure(figsize=(8,6))
    ax = fig.add_subplot(111, projection='3d');
    # Plot a basic wireframe.
    #ax.plot_wireframe(X, Y, W[:,:,j],color='b');
    ls = LightSource(315, 90) #cm.gist_earth
    rgb = ls.shade(W[:,:,j], cmap=cm.coolwarm, vert_exag=0.01, blend_mode='soft')
    surf = ax.plot_surface(X, Y, W[:,:,j],facecolors=rgb);
    #fig.colorbar(surf, shrink=0.5, aspect=10)
    ax.set_title("thermal conductivity graph at time t="+str(t[j]))
    ax.set_xlabel("X")
    ax.set_ylabel("Y")
    ax.set_zlabel("Temp")
    plt.show();

X_4, Y_4, T_4, W_4 = main3(u_0yt,u_1yt,u_x0t,uy_x1t,u_xy0,0,math.pi/2,0,math.pi,0.1,0,10,0.1)

ip_w.interact(plot_3d_time,  x = ip_w.fixed(X_4),  y = ip_w.fixed(Y_4), t=ip_w.fixed(T_4), W=ip_w.fixed(W_4), j=ip_w.IntSlider(min=0, max=len(T_4)-1, step=1))

# 5-ая лаба 2-ой вариант

def u_0_y_t(x,y,t):
    return cosh(y)*exp(-3*t)

def u_1_y_t(x,y,t):
    return 0

def u_x_0_t(x,y,t):
    return cos(2*x)*exp(-3*t)

def u_x_1_t(x,y,t):
    return (5/4)*cos(2*x)*exp(-3*t)

def u_x_y_0(x,y,t):
    return cos(2*x)*cosh(y)

def main4(u_0_y_t, u_1_y_t, u_x_0_t, u_x_1_t, u_x_y_0, x0, x1, y0, y1, h, t0, t1, tau):
    In = int((x1-x0)/h)
    Jn = int((y1-y0)/h)
    Kn = int((t1-t0)/tau)
    tau_2 = tau/2
    r=tau/(2*h*h)

    I = np.linspace(x0,x1,In) # массив с x
    J = np.linspace(y0,y1,Jn) # массив с y
    K = np.linspace(t0,t1,Kn) # массив с t

    w=np.zeros((Jn,In,Kn))

    # 0yt
    for k in range(Kn):
        for j in range(Jn):
            w[j, 0, k] = u_0_y_t(0, J[j], K[k])

    # 1yt  уже заполнен нулями

    # x0t
    for k in range(Kn):
        for i in range(In):
            w[0, i, k] = u_x_0_t(I[i], 0, K[k])

    # xy0
    for j in range(Jn):
        for i in range(In):
            w[j, i, 0] = u_x_y_0(I[i], J[j], 0)

    #x1t
    for k in range(Kn):
        for i in range(In):
            w[Jn-1, i, k] = u_x_1_t(I[i], Jn-1, K[k])



    for k in range (1,Kn): # идем по t

        P=np.zeros((Jn,In)) # промежуточный слой

        for j in range (1,Jn-1): # идем по y

            A=np.zeros((In-2,In-2))
            b=np.zeros(In-2)

            for n in range (0,In-2):
                A[n,n]= 2*r+1
            for n in range (0,In-3):
                A[n+1,n]= -r
                A[n,n+1]= -r
            for i in range (0,In-2):
                b[i] = w[j,i+1,k-1]

            b[0] += r*u_0_y_t(I[0],J[j],K[k]-tau_2)
            b[In-3] += r*u_1_y_t(I[In-1],J[j],K[k]-tau_2)

            p = solve(A, b)

            for i in range(0,In-2):
                P[j,i+1] = p[i]

            #P[j,0] = u_0_y_t(I[0],J[j],K[k]-tau_2)
            #P[j,In-1] = u_1_y_t(I[In-1],J[j],K[k]-tau_2)

        for j in range (1,Jn-1):   # идем по y 2 раз

            A=np.zeros((In-2,In-2))
            b=np.zeros(In-2)

            for n in range (0,In-2):
                A[n,n]=(2*r+1)
            for n in range (0,In-3):
                A[n+1,n]=-r
                A[n,n+1]=-r

            for i in range (0,In-2):
                b[i] = P[j,i+1]

            b[0] -= -r*u_0_y_t(I[0],J[j],K[k])

            b[In-3] -= -r*u_1_y_t(I[In-1],J[j],K[k])

            p = solve(A, b)

            for i in range(0,In-3):
                w[j,i+1,k] = p[i]

        """
        for i in range (1,In-1):   # идем по x

            A=np.zeros((Jn-2,Jn-2))
            b=np.zeros(Jn-2)

            for n in range (1,Jn-2):
                A[n,n] = 1+2*r
            for n in range (0,Jn-3):
                A[n+1,n] = -r
                A[n,n+1] = -r

            for j in range (0,Jn-2):
                b[j] = P[j+1,i]

            b[0] += r*w[j,0,k]
            b[Jn-3] += r*w[j,In-1,k]


            p = solve(A, b)

            for j in range(0,Jn-2):
                w[j+1,i,k] = p[j]
           """

    return I, J, K, w

X_5, Y_5, T_5, W_5 = main4(u_0_y_t, u_1_y_t, u_x_0_t, u_x_1_t, u_x_y_0, 0, pi/4, 0, log(2), 0.01, 0, 3, 0.1)

ip_w.interact(plot_3d_time,  x = ip_w.fixed(X_5),  y = ip_w.fixed(Y_5), t=ip_w.fixed(T_5), W=ip_w.fixed(W_5), j=ip_w.IntSlider(min=0, max=len(T_5)-1, step=1))

# 6-ая лаба 2-ой вариант

def K(x,s):
    return (x**2)*exp(x*s)

def func(x):
    return 1-x*(exp(x)-exp(-x))

def Fredholm_2_Rect(K, func, a, b, h):
    n = int((b-a)/h)
    x = np.linspace(a,b,n)
    wt = 0.5
    wj = 1
    A = np.zeros((n,n))
    b = np.zeros(n)
    for i in range(n):
        A[i,0] = -h*wt*K(x[i],x[0])
        for j in range(1,n-1):
            A[i,j] = -h*wj*K(x[i],x[j])
        A[i,n-1] = -h*wt*K(x[i],x[n-1])
        A[i,i] += 1
        b[i] = func(x[i])
    y = solve(A, b)
    return x , y

def main5(K, func, a, b, h):
    x, y = Fredholm_2_Rect(K, func, a, b, h)
    fig = plt.figure(figsize=(16,8))
    plt.plot(x,y,label='')
    plt.title('Approximate solution',fontsize=24)
    plt.xlabel('x')
    plt.ylabel('u(x)')
    plt.show()

main5(K, func, -1, 1, 0.001)