Difference between revisions of "Numerical differentiation and integration"

From Process Model Formulation and Solution: 3E4
Jump to navigation Jump to search
m
m
 
Line 2: Line 2:


Please [[Media:E-Numerical-differentiation-01-Nov-2010.pdf|download the lecture slides here]].
Please [[Media:E-Numerical-differentiation-01-Nov-2010.pdf|download the lecture slides here]].


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

Latest revision as of 21:30, 7 February 2017

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