# Numerical differentiation and integration

## Part B: Numerical integration (08 to 15 November 2010)

### Source code used in slides

Composite trapezoidal integration

MATLAB Python
```function In = trapezcomp(f, a, b, n)
%    Composite trapezoidal function integration
%
%    INPUTS:
%    f:  the function to integrate
%    a:  lower bound of integration
%    b:  upper bound
%    n:  number of panels to create between ``a`` and ``b``

% Initialization
h = (b-a)/n;
x = a;

% Composite rule
In =f(a);
for k=2:n
x  = x + h;
In = In + 2. * f(x);
end
In = (In + f(b)) * h * 0.5;
% end of function
```

Then use this function as follows:

```>> f = @(x)sin(x);
>> I1 = trapezcomp( f, 0, pi/2, 1 )    % 0.7854
>> I2 = trapezcomp( f, 0, pi/2, 2 )    % 0.9481
```
```import numpy as np

def trapezcomp(f, a, b, n):
"""
Composite trapezoidal function integration

INPUTS:
f:  the function to integrate
a:  lower bound of integration
b:  upper bound
n:  number of panels to create between ``a`` and ``b``
"""

# Initialization
h = (b - a) / n
x = a

# Composite rule
In = f(a)
for k in range(1, n):
x  = x + h
In += 2*f(x)

return (In + f(b))*h*0.5

if __name__ == '__main__':
def func(x):
return np.sin(x)

print(trapezcomp(func, 0, np.pi/2, 1))  # 0.785398163397
print(trapezcomp(func, 0, np.pi/2, 2))  # 0.948059448969
```

Romberg integration

MATLAB Python

In a file called trapezcomp.m

```function In = trapezcomp(f, a, b, n)
%    Composite trapezoidal function integration
%
%    INPUTS:
%    f:  the function to integrate
%    a:  lower bound of integration
%    b:  upper bound
%    n:  number of panels to create between ``a`` and ``b``

% Initialization
h = (b-a)/n;
x = a;

% Composite rule
In =f(a);
for k=2:n
x  = x + h;
In = In + 2. * f(x);
end
In = (In + f(b)) * h * 0.5;
% end of function
```

In a file called romberg.m

```function I = romberg(f, a, b, p)
% Romberg integration
%
%    INPUTS:
%    f:  the function to integrate
%    a:  lower bound of integration
%    b:  upper bound
%    p:  number of rows in the Romberg table

I = zeros(p, p);
for k=1:p
% Composite trapezoidal rule for 2^k panels
I(k,1) = trapezcomp(f, a, b, 2^(k-1));

% Romberg recursive formula
for j=1:k-1
I(k,j+1) = (4^j * I(k,j) - I(k-1,j)) / (4^j - 1);
end

disp(I(k,1:k));   % display intermediate results
end
% end of function
```

Then use this function as follows:

```>> f = @(x)sin(x);
>> format long
>> p_rows = 4;
>> I = romberg(f, 0, pi/2, p_rows);
>> solution = I(p_rows, p_rows);
>> disp(solution)  # 1.000000008144020
```
```import numpy as np

def trapezcomp(f, a, b, n):
"""
Composite trapezoidal function integration

INPUTS:
f:  the function to integrate
a:  lower bound of integration
b:  upper bound
n:  number of panels to create between ``a`` and ``b``
"""

# Initialization
h = (b - a) / n
x = a

# Composite rule
In = f(a)
for k in range(1, n):
x  = x + h
In += 2*f(x)

return (In + f(b))*h*0.5

def romberg(f, a, b, p):
"""
Romberg integration

INPUTS:
f:  the function to integrate
a:  lower bound of integration
b:  upper bound
p:  number of rows in the Romberg table
"""

I = np.zeros((p, p))
for k in range(0, p):
# Composite trapezoidal rule for 2^k panels
I[k, 0] = trapezcomp(f, a, b, 2**k)

# Romberg recursive formula
for j in range(0, k):
I[k, j+1] = (4**(j+1) * I[k, j] - I[k-1, j]) / (4**(j+1) - 1)

print(I[k,0:k+1])   # display intermediate results

return I

if __name__ == '__main__':
def func(x):
return np.sin(x)

p_rows = 4
I = romberg(func, 0, np.pi/2, p_rows)
solution = I[p_rows-1, p_rows-1]
print(solution)                  # 1.00000000814
```