Difference between revisions of "Numerical differentiation and integration"

From Process Model Formulation and Solution: 3E4
Jump to navigation Jump to search
m
 
(5 intermediate revisions by the same user not shown)
Line 1: Line 1:
== Part A: Numerical differentiation ==
== Part A: Numerical differentiation (01 to 04 November 2010) ==


<pdfreflow>   
Please [[Media:E-Numerical-differentiation-01-Nov-2010.pdf|download the lecture slides here]].
class_date        = 01 to 04 November <!-- (slides 1 to 11)<br>27 October (slides 10 to 19)<br>28 October (slides 20 to 31) -->
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 (08 to 15 November 2010) ==


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


== Part B: Numerical integration ==
=== Source code used in slides ===


<pdfreflow>   
'''Composite trapezoidal integration'''
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 ===


{| class="wikitable"
{| class="wikitable"
Line 81: Line 69:
     # Composite rule
     # Composite rule
     In = f(a)
     In = f(a)
     for k in range(2, n+1):
     for k in range(1, n):
         x  = x + h
         x  = x + h
         In += 2*f(x)
         In += 2*f(x)
Line 93: Line 81:
     print(trapezcomp(func, 0, np.pi/2, 1))  # 0.785398163397
     print(trapezcomp(func, 0, np.pi/2, 1))  # 0.785398163397
     print(trapezcomp(func, 0, np.pi/2, 2))  # 0.948059448969
     print(trapezcomp(func, 0, np.pi/2, 2))  # 0.948059448969
</syntaxhighlight>
|}
'''Romberg integration'''
{| class="wikitable"
|-
! MATLAB
! Python
|-
| width="50%" valign="top" |
In a file called '''<tt>trapezcomp.m</tt>'''
<syntaxhighlight lang="matlab">
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
</syntaxhighlight>
In a file called '''<tt>romberg.m</tt>'''
<syntaxhighlight lang="matlab">
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
</syntaxhighlight>
Then use this function as follows:
<syntaxhighlight lang="matlab">
>> 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
</syntaxhighlight>
| width="50%" valign="top" |
<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(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
</syntaxhighlight>
</syntaxhighlight>
|}
|}

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
Retrieved from "https://learnche.org/wiki_3E4/index.php?title=Numerical_differentiation_and_integration&oldid=1054"