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 55
7
Blurred image window size 2121
8
Blurred image window size 101101
9
Speed of a blurring operation
(since our image size is 674 1200 and window is 101101)
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 Simpsons 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 Simpsons method (1)Using parabola P(x)to estimate area under the curve f(x) 38 Simpsons method (2)Simpsons 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.Simpsons 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:Simpsons methodAdding areas of all the N segments gives an approximation to the integral.40 Simpsons 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 methodsFor all the three methods we saw, as the number of subintervals N approaches infinity, approximation improves and approaches actual integralMidpoint method:However, the speed of improvement differs: Trapezoidal method: Simpsons method:Simpsons method approaches the integral faster43 Integration using scipy.integrateSee example at scipyintegrate.py.44 Gaussian integral See example at gaussianintegral.py.45 Finding rootsof 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/920348To find the roots in difficult cases, we must generate a sequence of approximationsuntil 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.htmlA 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 problemThere 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.