#!/usr/bin/python
# -*- coding: utf-8 -*-

"""  Ejemplo del Metodo de Diferencias Finitas para resolver la ecuacion diferencial parcial:
    Uxx=-Pi*Pi*cos(Pix)
        xi<=U<=xf
        U(xi)=vi y U(xf)=vf
 con solucion analitica:
     cos(pix)

 Para conocer mas del metodo ver:
     http://132.248.182.159/acl/MDF

 Autor: Antonio Carrillo Ledesma
     http://mmc.geofisica.unam.mx/acl
"""


# 2.X compatible
from __future__ import print_function
import sys

if sys.version[0] == "2":
    input = raw_input
import math


def LadoDerecho(x):
    "Lado derecho de la ecuacion diferencial parcial"
    return -math.pi * math.pi * math.cos(math.pi * x)


def Jacobi(A, b, N, iter):
    "Resuelve Ax=b usando el metodo Jacobi"
    sum = 0
    x = [0] * N
    xt = [0] * N
    for m in range(iter):
        sum = A[0][2] * x[1]
        xt[0] = (1.0 / A[0][1]) * (b[0] - sum)
        for i in range(1, N - 1):
            sum = A[i][0] * x[i - 1] + A[i][2] * x[i + 1]
            xt[i] = (1.0 / A[i][1]) * (b[i] - sum)
        sum = A[N - 1][0] * x[N - 2]
        xt[N - 1] = (1.0 / A[N - 1][1]) * (b[N - 1] - sum)
        for i in range(N):
            x[i] = xt[i]
    return x


def GaussSeidel(A, b, N, iter):
    "Resuelve Ax=b usando el metodo GaussSeidel"
    sum = 0
    x = [0] * N
    for m in range(iter):
        sum = A[0][2] * x[1]
        x[0] = (1.0 / A[0][1]) * (b[0] - sum)
        for i in range(1, N - 1):
            sum = A[i][0] * x[i - 1] + A[i][2] * x[i + 1]
            x[i] = (1.0 / A[i][1]) * (b[i] - sum)
            sum = A[N - 1][0] * x[N - 2]
            x[N - 1] = (1.0 / A[N - 1][1]) * (b[N - 1] - sum)

    return x


# MDF1DDBAND
if __name__ == "__main__":
    xi = 0.0  # Inicio del diminio
    xf = 1.0  # Fin del dominio
    vi = 1.0  # Valor en la frontera xi
    vf = -1.0  # Valor en la frontera xf
    n = 11  # Particion
    N = n - 2  # Incognitas
    h = (xf - xi) / (n - 1)
    # Incremento en la malla

    # Declaracion de la matriz A y los vectores b y x
    A = []  # Matriz A
    for i in range(N):
        A.append([0] * 3)
    b = [0] * N  # Vector b
    x = [0] * N  # Vector x

    # Stencil
    R = 1 / (h * h)
    P = -2 / (h * h)
    Q = 1 / (h * h)

    # Primer renglon de la matriz A y vector b
    A[0][1] = P
    A[0][2] = Q
    b[0] = LadoDerecho(xi) - vi * R
    print(b)
    # Renglones intermedios de la matriz A y vector b
    for i in range(1, N - 1):
        A[i][0] = R
        A[i][1] = P
        A[i][2] = Q
        b[i] = LadoDerecho(xi + h * (i + 1))
    # Renglon final de la matriz A y vector b
    A[N - 1][0] = R
    A[N - 1][1] = P
    b[N - 1] = LadoDerecho(xi + h * N - 2) - vf * Q

    # Resuelve por el metodo Jacobi
    for i in range(N):
        for j in range(3):
            print(" ", A[i][j], end="")
        print("")

    print(b)
    x = Jacobi(A, b, N, N * 14)
    print(x)
    # Solucion
    print(xi, vi)
    for i in range(N):
        print(xi + h * (i + 1), x[i])
        print(xf, vf)

    # Resuelve por el metodo Gauss-Seidel
    for i in range(N):
        for j in range(3):
            print(" ", A[i][j], end="")
        print("")

    print(b)
    x = GaussSeidel(A, b, N, N * 7)
    print(x)
    # Solucion
    print(xi, vi)
    for i in range(N):
        print(xi + h * (i + 1), x[i])
    print(xf, vf)