<a href="https://colab.research.google.com/github/stephenbeckr/numerical-analysis-class/blob/master/Demos/Ch4_integration.ipynb" target="_parent"></a>

# Numerical Integration (quadrature)

- See also Prof. Brown's [integration notebook](https://github.com/cu-numcomp/numcomp-class/blob/master/Integration.ipynb) for CSCI-3656   [](https://colab.research.google.com/github/cu-numcomp/numcomp-class/blob/master/Integration.ipynb)
- Bengt Fornberg's talk [Gregory formulas and improving on the Trapezoidal rule](https://www.colorado.edu/amath/sites/default/files/attached-files/2019_unm_0.pdf)


```
import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import BarycentricInterpolator as interp

# From Table 9.2 in Quarteroni, Sacco and Saleri "Numerical Mathematics" (Springer, 2000)
ClosedNewtonCotesWeights = { 1:[1/2,1/2], 2:[1/3,4/3,1/3], 3:[3/8,9/8,9/8,3/8], 4:[14/45, 64/45, 24/45, 64/45, 14/45], 
    5:[95/288, 375/288,250/288, 250/288, 375/288, 95/288], 6:[41/140,216/140,27/140,272/140,27/140,216/140,41/140]}
ClosedNewtonCotesNames = {1:"n=1, Trapezoid", 2:"n=2, Simpson's", 3:"n=3, Simpson's 3/8", 4:"n=4, Boole's", 5:"n=5", 6:"n=6"}
```


```
f = lambda x : np.cos(x)
F = lambda x : np.sin(x)  # dF/dx = f
a,b = -1,2


# Other examples to try
# f = lambda x : x**(3/2)
# F = lambda x : 2/5*x**(5/2)
# a,b = 0,1

# f = lambda x : 1/(1+x**2)  # aka Runge's function
# F = lambda x : np.arctan(x)
# a,b = -5,5



I   = F(b) - F(a)
print("Integral I is {:.3f}".format(I))

x = np.linspace(a,b)
plt.fill_between( x, f(x), alpha=0.5);
plt.axvline(color='k');
plt.axhline(color='k');
```

### Try the Trapezoidal rule, n = 1


```
n = 1
print("Using the rule: ", ClosedNewtonCotesNames[n] )
weights   = ClosedNewtonCotesWeights[n]
(nodes,h) = np.linspace(a,b,n+1,retstep=True) # retstep tells it to return the spacing h

I_estimate = h*np.dot( weights, f(nodes) )
p          = interp(nodes,f(nodes))

x = np.linspace(a,b)
plt.fill_between( x, f(x), alpha=0.5);
plt.axvline(color='k');
plt.axhline(color='k');
plt.plot( x, p(x), 'r-', label="Interpolating polynomial" )
plt.legend()

print("True integral: {:.3f}, Estimate: {:.3f}, Abs. Error: {:.5f}".format(I,I_estimate,abs(I-I_estimate)))
nodes.tolist(),h,weights
```

### And Simpson's rule, n=2


```
n = 2
print("Using the rule: ", ClosedNewtonCotesNames[n] )
weights   = ClosedNewtonCotesWeights[n]
(nodes,h) = np.linspace(a,b,n+1,retstep=True) # retstep tells it to return the spacing h

I_estimate = h*np.dot( weights, f(nodes) )
p          = interp(nodes,f(nodes))

x = np.linspace(a,b)
plt.fill_between( x, f(x), alpha=0.5);
plt.axvline(color='k');
plt.axhline(color='k');
plt.plot( x, p(x), 'r-', label="Interpolating polynomial" )
plt.legend()

print("True integral: {:.3f}, Estimate: {:.3f}, Abs. Error: {:.5f}".format(I,I_estimate,abs(I-I_estimate)))
nodes.tolist(),h,weights
```

### n=3


```
n = 3
print("Using the rule: ", ClosedNewtonCotesNames[n] )
weights   = ClosedNewtonCotesWeights[n]
(nodes,h) = np.linspace(a,b,n+1,retstep=True) # retstep tells it to return the spacing h

I_estimate = h*np.dot( weights, f(nodes) )
p          = interp(nodes,f(nodes))

x = np.linspace(a,b)
plt.fill_between( x, f(x), alpha=0.5);
plt.axvline(color='k');
plt.axhline(color='k');
plt.plot( x, p(x), 'r-', label="Interpolating polynomial" )
plt.legend()

print("True integral: {:.3f}, Estimate: {:.3f}, Abs. Error: {:.5f}".format(I,I_estimate,abs(I-I_estimate)))
nodes.tolist(),h,weights
```

### n=4


```
n = 4
print("Using the rule: ", ClosedNewtonCotesNames[n] )
weights   = ClosedNewtonCotesWeights[n]
(nodes,h) = np.linspace(a,b,n+1,retstep=True) # retstep tells it to return the spacing h

I_estimate = h*np.dot( weights, f(nodes) )
p          = interp(nodes,f(nodes))

x = np.linspace(a,b)
plt.fill_between( x, f(x), alpha=0.5);
plt.axvline(color='k');
plt.axhline(color='k');
plt.plot( x, p(x), 'r-', label="Interpolating polynomial" )
plt.legend()

print("True integral: {:.3f}, Estimate: {:.3f}, Abs. Error: {:.5f}".format(I,I_estimate,abs(I-I_estimate)))
nodes.tolist(),h,weights
```

### n=5


```
n = 5
print("Using the rule: ", ClosedNewtonCotesNames[n] )
weights   = ClosedNewtonCotesWeights[n]
(nodes,h) = np.linspace(a,b,n+1,retstep=True) # retstep tells it to return the spacing h

I_estimate = h*np.dot( weights, f(nodes) )
p          = interp(nodes,f(nodes))

x = np.linspace(a,b)
plt.fill_between( x, f(x), alpha=0.5);
plt.axvline(color='k');
plt.axhline(color='k');
plt.plot( x, p(x), 'r-', label="Interpolating polynomial" )
plt.legend()

print("True integral: {:.3f}, Estimate: {:.3f}, Abs. Error: {:.5f}".format(I,I_estimate,abs(I-I_estimate)))
nodes.tolist(),h,weights
```

### n=6


```
n = 6
print("Using the rule: ", ClosedNewtonCotesNames[n] )
weights   = ClosedNewtonCotesWeights[n]
(nodes,h) = np.linspace(a,b,n+1,retstep=True) # retstep tells it to return the spacing h

I_estimate = h*np.dot( weights, f(nodes) )
p          = interp(nodes,f(nodes))

x = np.linspace(a,b)
plt.fill_between( x, f(x), alpha=0.5);
plt.axvline(color='k');
plt.axhline(color='k');
plt.plot( x, p(x), 'r-', label="Interpolating polynomial" )
plt.legend()

print("True integral: {:.3f}, Estimate: {:.3f}, Abs. Error: {:.5f}".format(I,I_estimate,abs(I-I_estimate)))
nodes.tolist(),h,weights
```

## Let's try different kinds of functions


```
def tryAllRules( f, F, a, b):
  err = []
  for n in range(1,6+1):
    weights   = ClosedNewtonCotesWeights[n]
    (nodes,h) = np.linspace(a,b,n+1,retstep=True)
    I_estimate = h*np.dot( weights, f(nodes) )
    I   = F(b) - F(a) # True answer
    err.append( abs(I_estimate - I))
  return np.array( err )


f = lambda x : np.cos(x)
F = lambda x : np.sin(x)  # dF/dx = f
a,b = -1,2
err1 = tryAllRules( f, F, a, b)

# Other examples to try
f = lambda x : x**(3/2)
F = lambda x : 2/5*x**(5/2)
a,b = 0,1
err2 = tryAllRules( f, F, a, b)

f = lambda x : x**(11/2)
F = lambda x : 2/13*x**(5/13)
a,b = 0,1
err3 = tryAllRules( f, F, a, b)

# Runge's function
f = lambda x : 1/(1+x**2)
F = lambda x : np.arctan(x)
a,b = -5,5
err4 = tryAllRules( f, F, a, b)

print("Rows are different n, columns are different functions")
print(np.array2string( np.array([err1,err2,err3,err4]).T, precision=2))
```

    Rows are different n, columns are different functions
    [[1.56e+00 1.00e-01 3.46e-01 2.36e+00]
     [6.65e-02 2.37e-03 2.76e-02 4.05e+00]
     [2.86e-02 1.29e-03 1.24e-02 6.65e-01]
     [9.15e-04 3.03e-04 1.22e-04 3.73e-01]
     [5.11e-04 2.11e-04 6.92e-05 4.39e-01]
     [1.04e-05 8.73e-05 5.90e-07 1.12e+00]]


### Let's examine Runge's function more closely
$$f(x) = \frac{1}{1+x^2}$$

Our error wasn't going down, but the function is $C^\infty(\mathbb{R})$. Did we make a mistake?

No, our formula was correct, the issue is that the $f'(\xi)$ term (and $f''(\xi)$, etc.) are very large.  One way to think of this issue is that the function has a **singularity** (though it is on the imaginary axis, at $\pm i$).

(Btw, how do you prounce Runge? It's German, and you can listen to native speakers say it [at Forvo](https://forvo.com/search/Runge/))


```
import sympy
from sympy.abc import x
from sympy import init_printing
from sympy.utilities.lambdify import lambdify
init_printing()
import matplotlib as mpl
mpl.rcParams['mathtext.fontset'] = 'cm'
mpl.rcParams.update({'font.size': 20})


f = lambda x : 1/(1+x**2)
F = lambda x : np.arctan(x)
a,b = -5,5

g = 1/(1+x**2) # symbolic version

gNumerical = lambdify(x,g) # avoid sympy plotting
xGrid   = np.linspace(a,b,150)
plt.figure(figsize=(10,8)) 
plt.plot( xGrid, gNumerical(xGrid),label='$f(x)$' )


#k = 3  # order of derivative
for k in range(1,6):
  dg = lambdify(x,sympy.diff(g,x,k))
  plt.plot( xGrid, dg(xGrid), label="$f^{("+str(k)+")}(x)$");
plt.axvline(color='k');
plt.axhline(color='k');
#plt.legend(prop={'size': 20});
plt.legend()
plt.title("Runge's function");

#sympy.plot(g);  # sympy plots are not so nice

# sympy.plot(sympy.diff(g,x,k));
```

4a1d80690b0f592ad469ed5004482cf20c638569
ipynb
Jupyter Notebook
Demos/Ch4_integration.ipynb
skhadem/numerical-analysis-class
a022fc2e73254800a1e193f94280446223ef01ae
BSD-3-Clause
2020-08-25T19:11:46.000Z
2021-11-12T21:49:44.000Z
Demos/Ch4_integration.ipynb
sabrinazhengliu/numerical-analysis-class
1585b68d3f3c375259507bf18f03ba0332851edb
BSD-3-Clause
2020-09-01T21:44:12.000Z
2020-09-01T21:44:12.000Z
Demos/Ch4_integration.ipynb
sabrinazhengliu/numerical-analysis-class
1585b68d3f3c375259507bf18f03ba0332851edb
BSD-3-Clause
2020-08-25T21:25:17.000Z
2021-03-03T04:14:26.000Z
Qwen/Qwen-72B
1. YES
2. YES