Difference between revisions of "Numerical differentiation and integration"
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 | ||
h = (b-a)/n; | h = (b-a)/n; | ||
x = a; | x = a; | ||
% | % 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
|