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

# # Trapezoidal rule when f'' is not finite

# In[17]:


import numpy as np


# The next function implements Trapezoid method.

# In[18]:


# Performs Trapezoid quadrature
def integrate(a,b,n,f):
    h = (b-a)/n
    x = np.linspace(a,b,n+1)
    y = f(x)
    res1 = 0.5*y[0] + np.sum(y[1:n]) + 0.5*y[n]
    return h*res1


# We will integrate
# $$
# f(x) = x \sqrt{x}, \qquad x \in [0,1]
# $$
# The exact integral is $2/5$ but $f''$ is not finite. Now we perform convergence test. The Peano kernel error formula gives the error estimate
# $$
# |E_n(f)| \le \frac{3 h^2}{16}
# $$

# In[19]:


f = lambda x: x * np.sqrt(x)
a, b = 0.0, 1.0
Ie = 2.0/5.0 # Exact integral
n, N = 2, 10

e = np.zeros(N)
for i in range(N):
    h = (b-a)/n
    I = integrate(a,b,n,f)
    e[i] = Ie - I
    if i > 0:
        print('%6d %18.8e %14.5g %18.8e'% (n,e[i],e[i-1]/e[i],3*h**2/16))
    else:
        print('%6d %18.8e %14.5g %18.8e'%(n,e[i],0,3*h**2/16))
    n = 2*n


# We still get $O(h^2)$ convergence and the error is smaller than the error estimate.