ENGM3356 Numerical Methods and Partial Differential EquationsWinter 2018Assignment #2: Solutions

ENGM3356 Numerical Methods and Partial Differential EquationsWinter 2018: pg 24)a) We use the logarithms because it ensures a function where the estimated radon levels neverbecome negative (which would have no physical meaning). Thus, although the logarithmsthemselves may become negative in a polynomial fit, the actual radon levels are obtainedby the transformationef(x), wheref(x) is the value of the fitted polynomial at the pointx.Now since we are attempting to map out the excursions above the 20 ppm level using rootfinding techniques, we want to convert the fitted polynomial to one in which the excursionscorrespond to roots. This is done by translating the curve by subtracting ln(2010−6) sothat our final function becomesf(x) =2.18 + 5.0768x2.6870x2+ 0.5129x30.0318x4This function has roots where the radon level crosses the 20 ppm threshold.b) The Matlab code for the Bisection method is as follows. The output from all four programswill be shown in Question #8.% assn2 Q4b (Bisection algorithm)tol= 0.0001; % relative accuracymaxit = 1000;% maximum number of iterations% find roots between 0 and 1, 3 and 4, 4 and 5, 7 and 8% first between 0 and 1:[r1, nit1, calls1] = Bisection(0, 1, tol, maxit);if nit1 == -1fprintf("No root found between 0 and 1\n");elsefprintf("Root 1: %f, took %d iterations and %d calls to f(x)\n", r1, nit1, calls1);end% second between 3 and 4[r2, nit2, calls2] = Bisection(3, 4, tol, maxit);if nit2 == -1fprintf("No root found between 3 and 4\n");elsefprintf("Root 2: %f, took %d iterations and %d calls to f(x)\n", r2, nit2, calls2);end% third between 4 and 5[r3, nit3, calls3] = Bisection(4, 5, tol, maxit);if nit3 == -1fprintf("No root found between 4 and 5\n");elsefprintf("Root 3: %f, took %d iterations and %d calls to f(x)\n", r3, nit3, calls3);end% and last between 7 and 8[r4, nit4, calls4] = Bisection(7, 8, tol, maxit);if nit4 == -1fprintf("No root found between 7 and 8\n");elsefprintf("Root 4: %f, took %d iterations and %d calls to f(x)\n", r4, nit4, calls4);endfprintf(’The radon level is too high in the interval [%f, %f]\n’, r1, r2);fprintf(’and in the interval [%f, %f]\n’,r3, r4);function [xr, nit, calls] = Bisection(xl, xu, tol, maxit)eps= 1.e-10;fl= f(xl);fu= f(xu);x0= xl;calls = 2;nit= 1;while nit < maxit%update root estimatexr =(xl + xu)/2;%estimate relative errorif abs(xr) < epse = abs((xr-x0)/eps);

ENGM3356 Numerical Methods and Partial Differential EquationsWinter 2018: pg 3elsee = abs((xr-x0)/xr);end%return if rel error is small enoughif e < tolreturn;end%compute f(xr) and save itfr = f(xr);calls = calls +1;if fr == 0% means we’ve already found the rootreturn;end%revise upper and lower bracketsif fl*fr < 0xu = xr;fu = fr;elseif fu*fr < 0xl = xr;fl = fr;elsenit = -1;return;endx0= xr;nit = nit + 1;endnit = -1;return;endfunction y = f(x)y = -2.18 + 5.0768*x - 2.6870*x.ˆ2 + 0.5129*x.ˆ3 - 0.0318*x.ˆ4;end4) The Matlab code for the False Position method is as follows:% assn2 Q5 (False Position algorithm)tol= 0.0001; % relative accuracymaxit = 1000;% maximum number of iterations% find roots between 0 and 1, 3 and 4, 4 and 5, 7 and 8% first between 0 and 1:[r1, nit1, calls1] = FalsePos(0, 1, tol, maxit);if nit1 == -1