#!/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)