Difference between revisions of "Numerical differentiation and integration"

From Process Model Formulation and Solution: 3E4
Jump to navigation Jump to search
Line 31: Line 31:
<syntaxhighlight lang="matlab">
<syntaxhighlight lang="matlab">
function In = trapezcomp( f, a, b, n )
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
% Initialization
h = (b-a)/n;
h = (b-a)/n;
x = a;
x = a;


% COMPOSITE RULE
% Composite rule
In =f(a);
In =f(a);
for k=2:n
for k=2:n
Line 55: Line 62:
| width="50%" valign="top" |
| width="50%" valign="top" |
<syntaxhighlight lang="python">
<syntaxhighlight lang="python">
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(2, n+1):
        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
</syntaxhighlight>
</syntaxhighlight>
|}
|}

Revision as of 13:20, 10 November 2010

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-08-Nov-2010.pdf </pdfreflow>

Source code used in slides

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(2, n+1):
        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