COMP 208
Lecture 21: Numerical integration
November 18, 2019 Instructor: Jonathan Campbell
• • •
Course updates
Midterm/A2 grades posted.
A4 coming out tomorrow.
Tutorials on Assignment 4 this week:
• •
Wednesday, Nov. 20, 11:30 a.m. – 12:30 p.m., TR 3120 TBD
2
• • •
Image analysis pt. 2
Numerical integration
Finding roots of non-linear equations, pt. 1
Plan for today
3
Image analysis (2)
Recall: blurring an image
•
Blurring is achieved by replacing each pixel by the mean (average) value of all the pixels surrounding it,
within a window of size NxN.
5
Original image
6
Blurred image — window size 5×5
7
Blurred image — window size 21×21
8
Blurred image — window size 101×101
9
Speed of a blurring operation
•
(since our image size is 674 × 1200 and window is 101×101)
•
•
•
•
Our blurring function is very slow
— around 25 billion operations (674 × 1200 × 101 × 101)
Compare the running time with the following code that uses the scipy package (make sure scipy is installed first).
meanfilter_scipy.py
The scipy code is much faster, as they have made many optimizations and numerical tricks in their code.
10
•
•
•
•
Edge detection is a method to identify a set of points (edges or lines) at which image brightness changes sharply.
Edge Detection
There are many functions in skimage for edge detection, but we will use the Canny algorithm implemented by skimage.feature.canny()
See canny_edge_detection.py.
Edge detection on sharper images tends to produce
too many lines. We can decrease these extraneous details by blurring the image before doing edge detection.
11
Original image
12
Image after edge detection
13
Blurred image with edge detection
14
Functions as arguments
Functions are also objects
def func1(x):
print(“This is func1”, x)
func1(10)
# This is func1 10
func2 = func1 func2(123)
#This is func1 123
print(type(func1), type(func2))
#
16
Functions can be passed as arguments
def func(g, x):
print(“Calling g(x): “, g(x))
def g1(x): return x * x
def g2(x):
return x * x * x
func(g1, 10) # Calling g(x): 100 func(g2, 10) # Calling g(x): 1000
17
Numerical integration
•
The definite integral of a function of a single variable, f(x), between two limits a and b
can be viewed as the area S under the curve of the function (and above the x-axis).
Definite integral
•
Numerical integration algorithms try to estimate this area.
Image source
• • •
We then estimate the area under the curve in each segment. Finally, sum these areas.
Numerical integration
•
Divide the region between a and b into N segments.
Our approach will be as follows:
•
Midpoint method
Trapezoidal method Simpson’s method
We will see three algorithms for estimating segment area:
• • •
20
Running example
We will use the following function in the interval [0, 5]
to compare the three algorithms.
21
• • •
We will divide the region from a to b into N equal rectangles. The width of each segment is
The first segment starts at a and the last segment ends at b.
Each segment starting point is calculated as follows:
Midpoint method
22
Width of each segment is
23
Midpoint method (2)
•
•
•
We estimate the area under the curve of f(x) by computing the areas of each rectangle (width * height at midpoint).
Then, to compute an approximation to the integral,
we just have to sum these areas.
The approximation will improve
as the number of segments is increased.
24
Midpoint method, n=10
25
Midpoint method, n=20
26
Midpoint method implementation
def midpoint_integrate(f, a, b, N): area = 0.0
dx = (b – a) / N # width of segment x = a # start segments at a
for i in range(N):
area += dx * f(x + dx / 2)
x += dx # go to next segment
return area
27
•
•
Midpoint method error
We know the integral of the example function is given by:
We can compare this true area value with what we found
to find the error of our midpoint method.
28
Midpoint method error (2)
def g(x):
return 1 / (1 + x ** 2)
actual_integral = 1.37340076695
estimated_integral = midpoint_integrate(g, 0, 5, 10)
print(estimated_integral) # 1.373543428316664 error = abs(estimated_integral – actual_integral) print(error) # 0.00014266136666418738
29
Midpoint method error (3)
N estimated_integral
—– ——————–
10 1.3735434283e+00
20 1.3734392602e+00
40 1.3734103959e+00
80 1.3734031745e+00
160 1.3734013689e+00
320 1.3734009174e+00
640 1.3734008046e+00
1280 1.3734007764e+00
2560 1.3734007693e+00
5120 1.3734007675e+00
10240 1.3734007671e+00
error
—————-
1.4266136666e-04
3.8493263529e-05
9.6289197895e-06
2.4075766312e-06
6.0191232953e-07
1.5047571300e-07
3.7615273785e-08
9.4000809359e-09
2.3462840559e-09
5.8283222693e-10
1.4196355203e-10
30
Area of a trapezoid
31
Trapezoidal method
•
•
• •
The trapezoidal method divides segments into trapezoids instead of rectangles.
Recall that the area of a trapezoid is given by
Thus, the area of each segment under the curve is given by:
We then sum the segment areas to get the approximation.
32
Width of each segment is
Trapezoidal method, n=10
34
Trapezoidal method implementation
def trapezoidal_integrate(f, a, b, N): area = 0.0
dx = (b – a) / N x=a
for i in range(N):
area += dx * (f(x) + f(x + dx)) / 2 x += dx
return area
35
Trapezoidal method error
N estimated_integral
—– ——————–
10 1.3731040812e+00
20 1.3733237548e+00
40 1.3733815075e+00
80 1.3733959517e+00
160 1.3733995631e+00
320 1.3734004660e+00
640 1.3734006917e+00
1280 1.3734007481e+00
2560 1.3734007622e+00
5120 1.3734007658e+00
10240 1.3734007667e+00
error
—————-
2.9668571989e-04
7.7012176612e-05
1.9259456542e-05
4.8152683754e-06
1.2038458712e-06
3.0096677106e-07
7.5245528253e-08
1.8815126790e-08
4.7075219278e-09
1.1806169375e-09
2.9889868358e-10
36
Smallest N for a given error tolerance
for N in range(1, 10000):
estimate = trapezoidal_integrate(g, 0, 5, N)
# Here we choose error tolerance to be 10e-8
if abs(estimate – actual_integral) < 10e-8: print(“Accuracy achieved at N =”, N) break37 Simpson’s method (1)Using parabola P(x)
to estimate area under the curve f(x) 38 Simpson’s method (2)•••Simpson’s method generally finds a better approximation to the area under the curve in each segment than the trapezoidal method (which uses a line instead of a parabola).Given any three points, there is a unique polynomial (parabola) called the interpolating polynomial, that passes through these points.Simpson’s method fits a parabola through the curve at three points — the value of the function at the two endpoints (a and b) and midpoint of the interval: 39 •The area under the parabola in the interval
is given by:Simpson’s method•Adding areas of all the N segments
gives an approximation to the integral.40 Simpson’s method implementationdef simpsons_integrate(f, a, b, N): area = 0.0dx = (b – a) / N x=afor i in range(N):area += dx * (f(x) + 4*f(x + dx / 2)
x += dx return area41+ f(x + dx)) / 6 Simpson’s method errorN estimated_integral—– ——————– 10 1.3733969793e+00 20 1.3734007584e+00 40 1.3734007664e+00 80 1.3734007669e+00160 1.3734007669e+00320 1.3734007669e+00640 1.3734007669e+00 1280 1.3734007669e+00 2560 1.3734007669e+00 5120 1.3734007669e+0010240 1.3734007669e+00error—————-3.7876621872e-068.5498519375e-095.3898707719e-103.8371528177e-117.0714545330e-125.1121329392e-124.9893422727e-124.9857895590e-124.9855675144e-124.9795723100e-124.9829029791e-1242 Accuracy of the three methods•For all the three methods we saw, as the number of subintervals N approaches infinity, approximation improves and approaches actual integral•Midpoint method:However, the speed of improvement differs: • • ••Trapezoidal method: Simpson’s method:Simpson’s method approaches the integral faster43 Integration using scipy.integrateSee example at scipyintegrate.py.44 Gaussian integral See example at gaussianintegral.py.45 Finding roots
of non-linear equations ••Often in engineering analysis, for a function f(x),
we need to find the values of x for which f(x) = 0.Roots of a functionThese values are known as the roots of the equation f(x),
also known as the zeroes or solutions of the function.47 Roots of a function (2)•••Some functions can be solved in closed form
(i.e., their solutions can be represented in terms of functions that we know)*.E.g.,: can be factored to (x-1)(x+1),
with the solution of f(x) = 0 being x = -1, 1.But some functions cannot be factored or represented in other ways.* ref: https://math.stackexchange.com/a/920348••To find the roots in difficult cases, we must generate a sequence of approximations
until we (hopefully) obtain a value close to the actual root.Roots of a function (3) Note that we will concentrate on finding the real roots to functions, not the complex ones.49 ref: http://www.xaktly.com/GraphingPolynomials.html•••A root where the value of y
does not change on either side.Double root E.g., for f(x) = x2,
f(x) is always positive.We will not be able
to detect these roots.51 ••In general, an equation may have any number of (real) roots.Number of rootsFor example, has only one root,
whereas has an infinite number of roots.52 Methods of finding roots••• •its physical context may suggest the approximate location.The methods we will see here all require a starting point,
i.e., an estimate of the root.•If the equation is associated with a physical problem
There are a few different ways of estimating a root value:a search can be carried out for estimated root values (we will see this).the function could be plotted (a visual procedure).53 •incremental search method• • • •bisection methodsecant method Newton-Raphson methodOutlineWe will see a few different methods next class.54 ••• •Numerical Methods in Engineering with Python (Kiusalaas)
Chapter 4 & 6; PDF copy available for McGill students
https://mcgill.on.worldcat.org/oclc/830085795Resources Introduction to Python for Science (Pine)
https://github.com/a301-teaching/pyman/tree/python3 Look over example code.Post on discussion boards if you have questions.55
Reviews
There are no reviews yet.