Numerical differentiation and integration

Part A: Numerical differentiation (01 to 04 November 2010)

Please download the lecture slides here.

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

Please download the lecture slides here.

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