# %% Import of packages:

import matplotlib.pyplot as plt
import numpy as np


# %% Function f(x):

def fun_f(x):
    y = 3*(x**2)   
    return y


# %% Generating x_array and k_array:

x_a = 0
x_b = 1
dx = 0.1
N = int((x_b - x_a)/dx + 1)
x_array = np.linspace(x_a, x_b, N)
k_array = np.arange(0, N)
 
# %% Trapezoid integration on sum form:

y_array = fun_f(x_array)
S_sum = dx*np.sum((y_array[0:-1] + y_array[1:]))/2

# %% Trapezoid integration on recursive form:

S_km1 = 0  # S_km1 means S_(k-1).

for k in k_array[1:]:
    y_km1 = fun_f(x_array[k-1])
    y_k = fun_f(x_array[k])
    S_k = S_km1 + dx*(y_km1 + y_k)/2
    S_km1 = S_k  # Shift to prepare for next iteration
        
S_recursive = S_k

# %% Presenting results:

print('S_sum =', f'{S_sum:.3f}')
print('S_recursive =', f'{S_recursive:.3f}')

# %% Plotting:

x_a = 0
x_b = 1
dx_cont = 0.001
N_cont = int((x_b - x_a)/dx_cont + 1)
x_cont_array = np.linspace(x_a, x_b, N_cont)
y_cont_array = fun_f(x_cont_array) 

plt.close('all')  # Closes all figures before plotting
plt.figure(num=1, figsize=(16/2.54, 16/2.54))

plt.plot(x_cont_array, y_cont_array, 'b', label='y_cont')
plt.plot(x_array, y_array, 'or', label='y')
plt.plot(x_array, y_array, '-g')
plt.legend()
plt.fill_between(x_array, 0, y_array, facecolor='y')
plt.grid(which='both', color='grey')
plt.xlabel('x')

# plt.savefig('trapezoid_integration_0_1.pdf')  # pdf
plt.show()