Numerical differentiation and integration
Jump to navigation
Jump to search
The printable version is no longer supported and may have rendering errors. Please update your browser bookmarks and please use the default browser print function instead.
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
|