'''The program produces accuracy-cost diagrams
for the Gaussian quadrature formulas of order 2, 4 and 6.

This script requires the file python13_subfunctions.py in order to work. 

Meaning of local variables:
b ... length of interval of integration
exact ... exact value of integral
N ... interval of integration is divided into 2**N subintervals
'''

import numpy as np
import matplotlib.pyplot as plt
from python13_subfunctions import gauss_ord2, gauss_ord4, gauss_ord6 

N = 10
error_list = np.zeros((3, N))


''' left image '''
b = 3
exact = np.sin(b)
for idx, k in enumerate(2**np.arange(N) + 1):
    ''' 
    idx = 0, 1, ..., N - 1
    k = 2**idx 
    '''
    
    x = np.linspace(0, b, k)
    
    error_list[0,idx] = np.abs(gauss_ord2(x,np.cos) - exact)
    error_list[1,idx] = np.abs(gauss_ord4(x,np.cos) - exact)
    error_list[2,idx] = np.abs(gauss_ord6(x,np.cos) - exact)


plt.subplot(121)
for idx in range(3):
     plt.loglog((idx+1)*2**np.arange(N), error_list[idx])
     plt.loglog((idx+1)*2**np.arange(N), error_list[idx], marker='o',
                           label = 'gauss '+str(2*(idx+1)))
plt.grid(True)
plt.title('cos(x)')
plt.xlabel('function evaluations')
plt.ylabel('error')
plt.legend(loc='lower left')
plt.axis([0.1, 1e4, 1e-18, 1])



''' right image '''
def f(x):
    return x**(5./2.)

b = 1.
exact = 2./7.*b**(7./2.)

for idx, k in enumerate(2**np.arange(N) + 1):
    ''' 
    idx = 0, 1, ..., N - 1
    k = 2**idx +1
    '''
    
    x = np.linspace(0, b, k)
    
    error_list[0,idx] = np.abs(gauss_ord2(x, f) - exact)
    error_list[1,idx] = np.abs(gauss_ord4(x, f) - exact)
    error_list[2,idx] = np.abs(gauss_ord6(x, f) - exact)


plt.subplot(122)
for idx in range(3):
    plt.loglog((idx+1)*2**np.arange(N), error_list[idx])
    plt.loglog((idx+1)*2**np.arange(N), error_list[idx], marker='o',
                           label = 'gauss '+str(2*(idx+1)))
plt.grid(True)
plt.title('$x^{5/2}$')
plt.xlabel('function evaluations')
plt.ylabel('error')
plt.legend(loc='lower left')
plt.axis([0.1, 1e4, 1e-18, 1])

plt.subplots_adjust(wspace=0.45) 
plt.show()