Unformatted Attachment Preview
Finding Roots and Extrema
CEE/MAE M20
Introduction to Computer Programming with MATLAB
Newton’s Method for Root Finding
How do we find the zero(s) of a
function using Newton’s method?
Goal
CEE/MAE M20
2
Newton’s Method for Root Finding
How do we find the zero(s) of a
function using Newton’s method?
1.) Make an initial guess x0
x0
CEE/MAE M20
3
Newton’s Method for Root Finding
How do we find the zero(s) of a
function using Newton’s method?
f (x0 )
1.) Make an initial guess x0
2.) Evaluate f(x0)
x0
CEE/MAE M20
4
Newton’s Method for Root Finding
How do we find the zero(s) of a
function using Newton’s method?
0
f (x0 )
f (x0 )
1.) Make an initial guess x0
2.) Evaluate f(x0)
3.) Calculate the derivative at x0
x0
CEE/MAE M20
5
Newton’s Method for Root Finding
How do we find the zero(s) of a
function using Newton’s method?
0
f (x0 )
f (x0 )
1.) Make an initial guess x0
2.) Evaluate f(x0)
3.) Calculate the derivative at x0
x0
CEE/MAE M20
6
Newton’s Method for Root Finding
How do we find the zero(s) of a
function using Newton’s method?
0
f (x0 )
f (x0 )
x1 = x0
1.) Make an initial guess x0
2.) Evaluate f(x0)
3.) Calculate the derivative at x0
4.) Calculate x1
f (x0 )
f 0 (x0 )
x0
CEE/MAE M20
7
Newton’s Method for Root Finding
How do we find the zero(s) of a
function using Newton’s method?
0
f (x0 )
f (x0 )
1.) Make an initial guess x0
2.) Evaluate f(x0)
3.) Calculate the derivative at x0
4.) Calculate x1
x1
x0
CEE/MAE M20
8
Newton’s Method for Root Finding
How do we find the zero(s) of a
function using Newton’s method?
0
f (x0 )
f (x0 )
1.) Make an initial guess x0
2.) Evaluate f(x0)
3.) Calculate the derivative at x0
4.) Calculate x1
5.) Evaluate f(x1)
6.) Calculate the derivative at x1
x1
x0
f (x1 )
0
f (x1 )
CEE/MAE M20
9
Newton’s Method for Root Finding
How do we find the zero(s) of a
function using Newton’s method?
0
f (x0 )
f (x0 )
x1
x0
1.) Make an initial guess x0
2.) Evaluate f(x0)
3.) Calculate the derivative at x0
4.) Calculate x1
5.) Evaluate f(x1)
6.) Calculate the derivative at x1
7.) Calculate x2…
Can we stop yet?
Iterate!
x2
f (x1 )
0
f (x1 )
CEE/MAE M20
10
Newton’s Method for Root Finding
What happens if we make a bad initial guess?
x0
f (x0 )
CEE/MAE M20
11
Newton’s Method for Root Finding
What happens if we make a bad initial guess?
x0
Root
f (x0 )
CEE/MAE M20
12
Newton’s Method for Root Finding
What happens if we make a bad initial guess?
x0
???
f (x0 )
CEE/MAE M20
13
Newton’s Method for Root Finding
How do we find the zero(s) of a
function using Newton’s method?
0
f (x0 )
f (x0 )
x1
x0
1.) Make an initial guess x0
2.) Evaluate f(x0)
3.) Calculate the derivative at x0
4.) Calculate x1
5.) Evaluate f(x1)
6.) Calculate the derivative at x1
7.) Calculate x2…
Can we stop yet?
Iterate!
x2
f (x1 )
0
f (x1 )
Let’s try to formulate this as pseudocode.
14
CEE/MAE M20
Function Handles
A way to pass the name of a function as an input to another function.
% A main script to call our
% function, MyFunction.m
function y = MyFunction(f, x, c)
% A simple function that evaluates
% f(x) and adds c.
x = pi;
c = 4;
y = f(x) + c;
y = MyFunction(@sin, x, c);
end
>> y = 4
MATLAB Review Functions as Parameters
What is the value of y when our main script is executed?
Suppose MyFunction is a function that expects a function as one of its input parameters. In a call
to MyFunction, the input function must be identified by its handle which is done by putting the
character “@” in front of its name.
15
CEE/MAE M20
Function Handles
We can pass our own functions by handle as well!
% A main script to call our
% function, MyFunction.m
function y = MyFunction(f, x, c)
% A simple function that evaluates
% f(x) and adds c.
x = pi;
c = 4;
y = f(x) + c;
y = MyFunction(@ComplicatedCubic, x, c);
end
This ComplicatedCubic
function is so simple, we don’t
need to store it in a separate
file…
function y = ComplicatedCubic(x)
% A cubic function we define ourselves
% (as opposed to simply using @sin)
y = x.^3 + x.^2 + x + 1;
end
CEE/MAE M20
16
Anonymous Function
We can store very simple functions without creating a
separate file using the idea of function handles.
% A main script to call our
% function, MyFunction.m
function y = MyFunction(f, x, c)
% A simple function that evaluates
% f(x) and adds c.
x = pi;
c = 4;
y = f(x) + c;
ComplicatedCubic = @(x) x.^3 + x.^2 + x + 1;
end
y = MyFunction(ComplicatedCubic, x, c);
Note: We don’t need the @ sign here,
because ComplicatedCubic is
already stored as a function handle.
17
CEE/MAE M20
CEE/MAE M20
Introduction to Computer Programming with MATLAB
H OMEWORK 8
Due: Tuesday, May 29, 2018, 11:55 pm
1. Spatial Hashing. In class, we discussed how to speed up the identification of particle-particle interactions like collision by discretizing the entire domain into a series of rectangular bins. In this problem,
you’ll construct your own approach to mapping particle data onto a 2D grid. The linear-indexing and
spatial conventions we’ll use are illustrated in Figure 1.
y max
5
9
13
2
6
10
14
3
7
11
15
4
8
12
16
...
1
2 y
y
0
0
x
2 x
...
x max
Figure 1: The location of each particle (blue dot) and bin-average location (red x) as calculated on a
4 × 4 grid, with xmax = ymax = 1 and h = 0.25. The following calculations assume bin numbering
starts in the upper-left corner, proceeds first by column, and that the lower-left corner of the grid
corresponds to (0, 0)
(a) Begin by defining two [1xN] vectors with uniformly distributed random values to represent the
x and y position of each particle. Set the size of your grid with variables xmax and ymax . You’re
free to scale and/or shift the particles to fall in any region of your domain, but you must choose
your limit variables such that 0 < xk < xmax and 0 < yk < ymax for k = 1, 2, . . . , N .
(b) Next, begin partitioning the grid. Specify the minimum cell dimension, h, that will be used along
with xmax and ymax to generate the integer number of bins in both the x- and y-directions, NX
and NY , respectively. Depending on your choice of h, the actual bin dimensions, ∆x and ∆y,
may be larger than requested (remember, the value h controls the minimum width and height).
(c) Construct a vector with the same size as x and y to hold the bin number of each particle. For
every particle k with x-position xk and y-position yk , the linearly-indexed bin number can be
calculated as
1
binNumk =
ymax − yk
xk
− 1 NY +
∆x
∆y
for k = 1, 2, ..., N
(1)
where d...e indicates rounding the bracketed quantity up to the nearest integer. Using the bin
numbers to classify points, calculate the average (x, y) position of the particles in each bin for
bins 1, 2, ..., NX × NY .
(d) Construct a single plot that displays both the location of each particle and the location of
each bin average as unconnected markers. You can plot the bin dividing lines as a visual
check of your algorithm by first calling grid on and then changing the interval spacing (e.g.,
xticks(0:dx:xMax)). Finally, display the particle ID’s contained within each bin to the
command window using exactly the formatting shown below (obtained using N = 50 and
rng('default')):
Bin
Bin
Bin
Bin
1:
2:
3:
4:
11 16 32 35 45
3 34 40
6 22 30
[]
..
.
Bin 16:
4 10 18 19 24
where we interpret the output as particle 11, described by x(11) and y(11), falls into Bin 1,
and so on. Test your code by generating results for both a relatively sparse data set, e.g., N = 20,
and dense data set, e.g., N = 200. In cases where a bin holds no particles, no average point
should appear on the plot and '[]' should be printed as part of the command window output.
(e) How does your code classify a point xk = 0, yk = ymax ? What bin should this particle logically
be assigned to, and (in words or pseudocode only) how could you modify your algorithm to
catch cases like these? Suppose each of these particles has a maximum interaction radius of
h (i.e., each particle only cares about other particles if their separation distance is ≤ h). How
many bins filled with particles do you need to check to identify all particles within distance h of
a particle located in Bin 1?
Note: Although you’re encouraged to adjust the partitioning values and data set to get a feel for
how spatial hashing works, please submit your code with N = 20, h = 0.25, and values uniformly
distributed on a xmax = ymax = 1 grid.
2. Newton’s Method. A zero-finding method of great interest is Newton’s method. It uses both the value
of the function and its derivative at each step. Suppose f and its derivative are defined at xc , where
“c” denotes current. The tangent line to f at xc has a zero at x+ = xc − f (xc )/f 0 (xc ) assuming that
f 0 (xc ) 6= 0. Newton’s method for finding a zero of the function f is based on the repeated use of this
formula. Complete the following function so that it performs as specified.
function [xc, fEvals] = Newton(f, x0, delta, fEvalMax)
% f is a HANDLE to a continuous function, f(x), of a single variable.
% x0 is an initial guess to a root of f.
% delta is a positive real number.
% fEvalsMax is a positive integer >= 2 that indicates the maximum
% number of f-evaluations allowed.
%
% Newton's method is repeatedly applied until the current iterate, xc,
% has the property that |f(xc)|