MACM 316- D100 COMPUTING ASSIGNMENT 3

• Please follow the Guidelines for Computing Assignments found on the Canvas FAQ.

• Submit a one-page report via Crowdmark, grading scheme is also found on the FAQ.

Root Finding 2D Contour

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