Numerical differentiation and integration

From Process Model Formulation and Solution

Jump to: navigation, search

[edit] Part A: Numerical differentiation

Class date: 01 to 04 November
I want my notes with:  

  pages per physical page

Use page frames?


[edit] Part B: Numerical integration

Class date: 08 to 15 November
I want my notes with:  

  pages per physical page

Use page frames?

[edit] 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
Personal tools
Namespaces
Variants
Actions
Navigation
Toolbox