Numerical differentiation and integration
Revision as of 14:27, 15 November 2010 by Kevindunn (talk | contribs) (→Source code used in slides)
Part A: Numerical differentiation
<pdfreflow> class_date = 01 to 04 November button_label = Create my course notes! show_page_layout = 1 show_frame_option = 1 pdf_file = E-Numerical-differentiation-01-Nov-2010.pdf </pdfreflow>
Part B: Numerical integration
<pdfreflow> class_date = 08 to 15 November button_label = Create my course notes! show_page_layout = 1 show_frame_option = 1 pdf_file = E-Numerical-integration-10-Nov-2010.pdf </pdfreflow>
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
|