Computational Methods in physics,, Matlab

Queens College, CUNY Physics Department

Physics 661 – Computational Methods in Physics Spring 2019 – Dr. David Goldberg

Page 1/3

HOMEWORK ASSIGNMENT 2 Due: Wednesday, 2/20/19, before 4PM

Bisection Method Part 1 – Specific Tolerance

• Create a function, bisect01, that outputs the zero of a function, 𝑓𝑓(𝑥𝑥), using the bisection method that uses a while loop. The function should throw an error if, 𝑓𝑓(𝑎𝑎)𝑓𝑓(𝑏𝑏) < 0, is not satisfied.

• Bisect01 should be called by: bisect01(f, a, b, tol), where: o f: the function handle to 𝑓𝑓(𝑥𝑥). o a: the lower bound of the interval for the root. o b: the lower bound of the interval for the root. o tol: the tolerance to how close to zero is sufficient

• For example, to find the root of 𝑓𝑓(𝑥𝑥) = cos 𝑥𝑥 − 𝑥𝑥, on the interval [0, 1], to a tolerance of 10−10, you would define the function by either a function file:

function y = func01(x) y = cos(x) - x; end

o Or defining it in the command window by: func01 = @(x) cos(x) - x;

o The bisection function would then be called by: bisect01(func01, 0, 1, 1e-10)

Part 2 – Specific Number of Decimal Places • Create a function, bisect02, that outputs the zero of a function, 𝑓𝑓(𝑥𝑥), using the bisection method.

The function should throw an error if, 𝑓𝑓(𝑎𝑎)𝑓𝑓(𝑏𝑏) < 0, is not satisfied. • Bisect02 should be called by: bisect02(f, a, b, p), where:

o f: the function handle to 𝑓𝑓(𝑥𝑥). o a: the lower bound of the interval for the root. o b: the lower bound of the interval for the root. o p: the number of decimal places for the accuracy of the root

• For example, to find the root of 𝑓𝑓(𝑥𝑥) = cos 𝑥𝑥 − 𝑥𝑥, on the interval [0, 1], to 6 decimal places, you would define the function as in Part 1, then called it by: bisect02(func01, 0, 1, 6)

Physics 661 Spring 2019

Dr. David Goldberg

[email protected]

Page 2/3

Part 3 – default Options Explore the possibility of defining default options for the tolerance of number of decimal places. Look up the function nargin in the MATLAB help for examples on this. The help document uses the switch statement as an alternative to the if statement. Use this information to modify your previous functions as follows:

• Set a default tolerance of bisect01 to a specific value, say 10−10, but allow the user to change it. Then function could then be called in two different ways:

o As in Part 1 where the tolerance is set by the last argument: bisect01(f, a, b, tol)

o By the same function just leaving out the last argument: bisect01(f, a, b)

where the function will use a default value for tol. • Set a default number of decimal places of bisect02 to a specific value, say 6, but allow the user to

change it. Then function could then be called in two different ways: o As in Part 2 where the number of decimal places is set by the last argument:

bisect01(f, a, b, p) o By the same function just leaving out the last argument:

bisect01(f, a, b) where the function will use a default value for p.

For your reference, the script below is what was developed in class:

clear f = @(x) x.^3 + x - 1; a = 0; b = 1; fa = f(a); fb = f(b); tol = 1e-10; n = 200; for i = 1:n c = (a+b)/2; fc = f(c); if abs(fc) < tol break end if fa*fc < 0 b = c; fb = fc; else a = c; fa = fc; end end

Physics 661 Spring 2019

Dr. David Goldberg

[email protected]

Page 3/3

Other Root Finding Methods Part 1 – Creating Root Finding Functions

• Create functions that outputs the result of a following rooting finding algorithms: o Fixed Point Iteration, FPI – Input should include:

 The function handle of the function of which to find a root.  Initial guess for the root.  The number of iterations.

o Newton-Raphson Method – Input should include:  The function handle of the function of which to find a root.  The function handle of the derivative of function of which to find a root.  Initial guess for the root.  The number of iterations.

o Secant Method – Input should include:  The function handle of the function of which to find a root.  Two distinct initial guesses for the root. Note, you function should check that the

two points are distinct and throw an error if not.  The number of iterations.

• All functions should include an appropriate stopping condition.

Part 2 – Data Output Modify the functions in part 1 to include an additional input option to display the data to the command window. For example, you can have a function fpi(f, x0) and you can also allow for an additional input, fpi(f, x0, 'disp'), where if the last input argument is supplied, it will display the relevant data. An example of this functionality, run the following code:

x = 0:0.1:1; y = exp(x); fprintf('%12s %6s %12s\n','Iteration #','x', 'exp(x)') for i = 1:length(x) fprintf('%12u %6.2f %12.8f\n',i,x(i), y(i)) end

More detailed information is available in the MATLAB help documentation for the fprintf function: https://www.mathworks.com/help/matlab/ref/fprintf.html#input_argument_d119e322652

Your submission should include the m-files of your functions and a documentation file (docx, pdf, etc.) that includes:

• Examples of different functions used and their output. • Compare the different methods for the same functions and initial guess. • If possible, compare the output to the expected value of the root. • Give examples of when your functions behave well and as well as badly.

Note: the modified functions from part 2 should be helpful for this documentation.