Description
One dimensional root finding has many applications, and a common situation is called continuation,
where a sequence of (related) roots are to be solved for. An elementary example of a continuation is the
numerical construction of a contour (or level curve) of a 2D function. A contour plot is a 2D representation
of a function of two variables, f(x, y) that plots all curves of the form f(x, y) = c, for various constants c.
Figure 1 shows a colour plot of the HI(x, y) function, including the zero contour (HI(x, y) = 0 is the white
dashed curve) . While it may seem that finding solutions to HI(x, y) = 0 is a 2D problem, there is a way
to apply the 1D methods presented in class.
-2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5
x
-2.5
-2
-1.5
-1
-0.5
0
0.5
1
1.5
2
2.5
y
-1 -0.8 -0.6 -0.4 -0.2 0 0.2
x
0
0.2
0.4
0.6
0.8
1
1.2
y
P0
P1
P2
P3
0
1
2
Figure 1: and Figure 2
Suppose we have found the points {(x1, y1), … (xn, yn)} on our zero contour. To find our next point,
(xn+1, yn+1), we will root search only on a circle, centred at Pn = (xn, yn), with a radius ∆s. The points
on this circle are a function of angle θn, and can write Pn+1 = (xn+1, yn+1) as
xn+1(θn) = xn + ∆s cos(θn) ; yn+1(θn) = yn + ∆s sin(θn) . (1)
We can find θn by solving the equation HI(xn+1(θn), yn+1(θn)) = 0, which is now a 1D root-finding problem
in θn. The radius ∆s is a numerical parameter for this procedure, and as you will discover, controls the
graphical quality of the contour produced. Figure 1 shows the zero contour of HI(x, y) computed by this
method (in red). The computation uses only ∆s = 0.6 and finds 24 points on the contour, so the computed
contour does not appear smooth.
A cartoon for our procedure is Figure 2 where the point P3 on the contour (thick blue curve) has just been
located at the angle θ2. Note that there is actually a second point on the (green) circle around P2, but it
is the previous point P1! Hence, special care needs to be taken that our root-finding algorithm doesn’t
DJM 1
MACM 316- D100 COMPUTING ASSIGNMENT 3
reverse itself. We can now repeat the procedure on the (magenta) circle around P3 to iterate the search
procedure to find P4.
But how can we find the first point? We can use the same angle solve procedure provided that we begin
from a point P0 that only needs to be within a distance ∆s from our contour. Note in the cartoon, P0 is
not on the contour, but the general procedure still works. For this assignment, you can manually set P0
as an ’eyeballed’ point that sits near your contour.
The goal of this assignment is to compare the bisection and secant methods, in terms of efficiency and robustness, when computing the zero contour of 2D functions. Begin by downloading
CA3 demo.m, a code that does the root-finding using Matlab’s fzero command — you will need to replace
fzero calls inside the contour finding loop. Make sure you download the most recent versions of BiSsqrt.m
and SMsqrt.m from Canvas to help you modify CA3 demo.m. We will benchmark our code on the HI(x, y)
function.
Both the bisection and secant methods need an initial interval — for this you will need to think about
the geometry. It turns out that the “right” interval for finding the angle θn is [θn−1 − δ, θn−1 + δ], where
δ gives HI(x(θn), y(θn)) opposite signs at the endpoints. You should start by setting δ = π/2 for both
methods, but you will be required to experiment with this parameter. Consider all of the following when
completing your report:
• Start by reproducing Figure 1 using your bisection code, with ∆s = 0.6 and δ = 3π/4.
• If you try to redo the above with secant method, the method will fail for the same choice of ∆s (try it
and see). Instead, set ∆s = 0.1 and δ = π/2 to trace the contour using the secant method.
• Test both root finders for a variety of ∆s values smaller than 0.1. What value of ∆s is needed for a
visually smooth curve, say for a plot as printed on a full letter-size page?
• Using half the ∆s value you found above, experiment next by changing δ. Your discussion should
report on any δ values that produce convergence failures, as well as how computational cost of
tracing the contour changes with δ. Cost should be considered as the average number of function
evaluations required per point traced on the contour.
• What have you learned about the trade-offs between using bisection or the secant method for tracing
contours?
Your report should include a plot of the HI(x, y) function, with your computed zero contour traced in red.
Make sure you include which root finding method you used to create this plot, as well as the values of ∆s
and δ.
Optional: If you wish, instead of tracing contours of the HI function, you may instead trace contours of
the cellular function:
f(x, y) = sin(πx) sin(πy) + a1 cos
π(πx − y)
2
+ a2 cos
π(x − πy)
2
where you should choose your own personal a1, a2 coefficients. The file CA3 newfunc.m produces contour
plots of these functions, with the zero contours plotted as white dashes. You may use a function of this
form to complete the δ testing of each method; if you do this, then the plot you include with your report
should be of f(x, y) instead of HI(x, y).
DJM 2