#!/usr/bin/env python3
# -*- coding: utf-8 -*-


#Constructs matrix and outputs 2nd largest eig 

import numpy as np
import matplotlib.pyplot as plt


def build_matrix(a,N): 
    Nt=4*N**2+1
    P=np.zeros((Nt,Nt))

    for i in range(Nt):
        xcurr=i-2*N**2
        P[i,i]=1-2*a
        if (i<Nt-1):
            P[i,i+1]=(1-xcurr/2/N**2)*a*(np.abs(xcurr)<=2*N**2)
        if (i>0):
            P[i,i-1]=(1+xcurr/2/N**2)*a*(np.abs(xcurr)<=2*N**2)
    # obtains second largest eig
    e,v=np.linalg.eig(P)
    e.sort()
    return P,e[-2]

a=1/4
N=[1,2,4,8,10,20]
for n in N:
    _,e=build_matrix(a,n)
    e=np.real(e) # to remove the 0j imaginary part 
    print('N = '+str(n) + ': expected rate of convergence (in TV norm): C*|lambda_2|^n = C*'+str(round(e,6))+'^n')