Notes.

When uploading your work, be sure to choose correct pages, without missing any, for each problem. Besides, make sure all pages are in order, no page is rotated, and none of your code or hand-written work is cut out.

According to our homework policies, you are allowed to collaborate for homework assignments in a group of up to 5 students. If you did, please list all your collaborators at the top of your homework, below your name.

Homework must be prepared using a live script (.mlx file). Submission generated using any other means (e.g., Words with MATLAB screenshots) will not be accepted.

Starting from this assignment, you will be asked to write MATLAB functions. Unlike scripts, functions can be included inside your mlx file at the very end of the document. So you no longer need to write an external m-file and print its content using the type function. For more detailed instruction, please read

[Guide] How to Embed Functions in Live Scripts:

Method 1: Using Code Example Environment

Suppose that you have written a function as a separate function m-file, myfoo.m and now you want to include its contents within your homework live script, HW2_Kim.mlx. At the top of MATLAB window, click on Insert > Code Example. It creates a blank box. Any text typed within this box will be formatted in monospaced fonts with syntax highlighting.

Drawbacks:

The box is not an executable code block.

If you make changes to myfoo.m, you need update the code example block accordingly (by copying and pasting).

Method 2: Using TYPE Function

Still suppose that you have a function written as a separate m-file titled myfoo.m. Then you can simply type

type myfoo

within a code block. The contents of myfoo.m will be printed out in your live script.

Note:

We have been using this to print out the contents of external script m-files. This works in the exact same way for function m-files.

Drawbacks:

If you forget to run that block after your final edit, your live script will show an outdated version of your function.

As you complete more and more assignments, your folders will be filled with more m-files and may be difficult to manage unless you are well-organized.

Method 3 (Recommended): Writing Functions at the End of Live Script

Instead of writing a separate m-file, you may just write the function at the very bottom of the document. The functions are local to the live script, and so you are able to call the function anywhere inside the live script. This way, your mlx file is fully self-contained and you no longer need to go back and forth between Editor and Live Editor. See the live script accompanying Lecture 9; all functions defined for the live script are gathered in the last section titled “Functions Used”.

This is my recommended method!

The only caveat I can think of:

Suppose Problem 1(a) asks you to write a function, and so you write it at the end of your live script. You then exports it to a pdf file and find out that Problem 1 solutions are on pp. 1–2, with the function written for 1(a) appearing on p. 7. In this case, when you upload the pdf file to Gradescope, you have to select pages 7, 1, and 2 for Problem 1, in that order.

Guide video link:

Updated on September 26, 2021

Math 3607: Homework 4

Due: 10:00PM, Tuesday, September 28, 2021

TOTAL: 30 points

• Problems marked with v are to be done by hand; those marked with � are to be solved using a

computer.

• Important note. Do not use Symbolic Math Toolbox. Any work done using sym or syms will receive

NO credit.

• Another important note. Starting from this assignment, you will be asked to write MATLAB

functions. Instead of writing an external function m-file, include all your functions at the end of your

live script.

1. (Sliders moving along grooves; adapted from LM 2.1–12 and Sample HW01) The mechanical

device shown in Figure 1 consists of two grooves in which sliders slide. These sliders are

connected to a straight rod.

`

L

A

B

Cpx,yq

W

O

Figure 1: The bronze square is a piece

of metal with two grooves cut out of

it as shown. There are sliders at the

points A and B which slide in these

grooves. The slider at A can only slide

vertically, and the one at B can only

slide horizontally. There is a straight

rod attached at A and B, which ex-

tends to C. As the point C moves

around the block, it traces out a closed

curve.

(a) v Analytically, determine the curve which is traced out by C in one rotation.

Suggestion. Let px,yq be the coordinates of the point C. Express the variables x,y in

terms of L, `, and θ, where θ P r0, 2πq is the angle from the part of the horizontal

groove which is to the right of B to the rod BC.

(b) � Using the previous result, plot the trajectory of C in one rotation for ` “ 2 and L “ 7.

2. (Spiral triangles to spiral polygons; adapted from LM 5.9–7, 6.8–34) � The following script1

generates spirals using equilateral triangles as shown in the figure below.

1It is slightly modified from the code included in Lecture 9 slides. Note the introduction of a new variable d_rot,

which is accountable for the rotation of the innermost triangle.

1

m = 21; d_angle = 4.5; d_rot = 90;

th = linspace(0, 360, 4) + d_rot;

V = [cosd(th);

sind(th)];

C = colormap(hsv(m));

s = sind(150 – abs(d_angle))/sind(30);

R = [cosd(d_angle) -sind(d_angle);

sind(d_angle) cosd(d_angle)];

hold off

for i = 1:m

if i > 1

V = s*R*V;

end

plot(V(1,:), V(2,:), ’Color’, C(i,:))

hold on

end

set(gcf, ’Color’, ’w’)

axis equal, axis off

Figure 2: Spiral triangles with

m “ 21 and θ “ 4.5˝.

(a) Write a function named spiralgon by modifying the script so that it generates spirals

using m regular n-gons for any n ě 3. Your function must be written at the end of

your homework live script (.mlx) file. Begin the function with the following header and

comments.

function V = spiralgon(n, m, d_angle, d_rot)

% SPIRALGON plots spiraling regular n-gons

% input: n = the number of vertices

% m = the number of regular n-gons

% d_angle = the degree angle between successive n-gons

% (can be positive or negative)

% d_rot = the degree angle by which the innermost n-gon

% is rotated

% output: V = the vertices of the outermost n-gon

….

(b) Run the statements below to generate some aesthetic shapes.

clf

subplot(2, 2, 1), spiralgon(3, 41, 4.5, -90);

subplot(2, 2, 2), spiralgon(4, 37, -2.5, 45);

subplot(2, 2, 3), spiralgon(5, 61, 3, -90);

subplot(2, 2, 4), spiralgon(8, 91, -4, 22.5);

Note. Copy the five lines, paste them inside a single code block, and run it. This code

block must precede your function(s).

3. (Machine epsilon; adapted from LM 9.3–3(a)) � Recall that the number in the computer

which follows 1 is 1 ` eps , which can be verified in MATLAB by

>> format long

2

>> (1 + 0.5*eps) – 1

ans =

0

>> (1 + 0.51*eps) – 1

ans =

2.220446049250313e-16

In the same manner:

(a) Verify that the number in the computer which follows 8 is 8 ` 8 eps by numerically

calculating 8 ` 4 eps and 8 ` 4.01 eps .

(b) Verify that the number in the computer which precedes 16 is 16 ´ 8 eps by numerically

calculating 16 ´ 4.01 eps and 16 ´ 4 eps .

(c) What are the numbers in the computer that precedes and follows 210 “ 1024, respec-

tively? Verify your claims in MATLAB by carrying out appropriate calculations.

Note. Begin with format long as shown in the example above. This is needed only once

before the beginning of part (a).

Note. Answer each part of the problem in a single code block. No external script needs to

be written.

4. (Catastrophic cancellation; LM 9.3–10) We revisit the function from Problem 3 of Homework

3. Consider the function

fpxq “

$

’

&

’

%

ex ´ 1

x

if x ‰ 0

1 if x “ 0,

where we are interested in exploring the catastrophic cancellation which occurs as x Ñ 0

since ex Ñ 1 as x Ñ 0.

(a) v Use the Taylor series expansion of ex to prove that f is continuous at 0.

(b) � Now calculate fpxq numerically for x “ 10´k where k P Nr1, 20s in three slightly

different ways:

i. Calculate fpxq as written.

ii. Calculate it as

f1pxq “

ex ´ 1

log ex

, for x ‰ 0.

(You and I know that analytically f1pxq ” fpxq for all nonzero x – but MATLAB

doesn’t.)

iii. MATLAB has a function which analytically subtracts 1 from the exponential to avoid

catastrophic cancellation before the result is calculated numerically. So define the

function f2pxq to be the same as fpxq except that ex ´ 1 is replaced by expm1(x).

Tabulate the results using disp or fprintf. The table should have four columns with

the first being x, the second using fpxq, the third using f1pxq, and the fourth using f2pxq,

with all shown to full accuracy. Do it as efficiently as you can, without using a loop.

Note. Write your code for this part in a single code block. No external script needs to

be written.

3

(c) v Comment on the results obtained in the previous part. Explain why certain methods

work well while others do not.

5. (Inverting hyperbolic cosine; FNC 1.3.6) The function

x “ coshptq “

et ` e´t

2

can be inverted to yield a formula for acoshpxq:

t “ log

´

x ´

a

x2 ´ 1

¯

. (‹)

In MATLAB, let t=-4:-4:-16 and x=cosh(t).

(a) v� Find the condition number of the problem fpxq “ acoshpxq by hand. (You may

use Equation (‹), or look up a formula for f 1 in a calculus book.) Then evaluate κf at

the elements of x in MATLAB.

(b) � Evaluate the right-hand side of Equation (‹) using x to approximate t. Record the

accuracy of the answers (by displaying absolute and/or relative errors), and explain.

(Warning: Use format long to get enough digits or use fprintf with a suitable

format.)

(c) � An alternate formula for acoshpxq is

t “ ´2 log

˜

c

x ` 1

2

`

c

x ´ 1

2

¸

. (:)

Apply Equation (:) to x and record the accuracy as in part (b). Comment on your

observation.

(d) v Based on your experiments, which of the formulas (‹) and (:) is unstable? What is

the problem with that formula?

Note. Write your code for each of parts (a), (b), and (c) in a single code block. No external

script needs to be written.

4

function x = ieee2fp(ieee)

% calculates the floating-point number from its ieee representation

% input: ieee = 8-character hexidecimal string or

% 64-bit binary string

% = (structure) when representation is in scientific notation

% +- 1.xxxxx… 2^exponent or

% +- 0.xxxxx… 2^(-1023)

% .sign = +1 or -1 as a fp number, not text

% .exponent = the exponent shown above, must be in [-1023, 1024]

% .mantissa = 13-character hexidecimal string or

% 52-bit binary string

bin2hex = [(0:9)+’0′, (0:5)+’a’];

if ischar(ieee)

if length(ieee) == 16

if any( (ieee < '0' | ieee > ‘f’) | (ieee > ‘9’ & ieee < 'a') )
error(['***** hexidecimal string not all 0-9 or a-f: ', ieee])
end
x = hex2num(ieee);
elseif length(ieee) == 64
if any( ieee < '0' | ieee > ‘1’ )

error([‘***** binary string not all 0-1: ‘, ieee])

end

Ieee = reshape(ieee, 4, 16)’;

Ieee_dec = bin2dec(Ieee);

ieee_hex = bin2hex(Ieee_dec+1);

x = hex2num(char(ieee_hex));

else

error([‘***** string not correct length: ‘, ieee])

end

else

if ~isstruct(ieee)

error([‘***** input is not a structure variable: ‘, num2str(ieee)])

end

if ieee.sign == +1

exponent_bin(1) = ‘0’;

elseif ieee.sign == -1

exponent_bin(1) = ‘1’;

else

error([‘***** incorrrect sign: ‘, num2str(ieee.sign)])

end

if mod(ieee.exponent, 1) ~= 0 || ieee.exponent < -(2^10 - 1) || ieee.exponent > 2^10

error([‘***** incorrect exponent: ‘, …

num2str(ieee.exponent)])

end

exponent_bin(2:12) = dec2bin(ieee.exponent+(2^10-1), 11);

if length(ieee.mantissa) == 52

if any( ieee.mantissa<'0' | ieee.mantissa>‘1’ )

error([‘***** binary string not all 0-1: ‘, ieee])

end

x = ieee2fp([exponent_bin, ieee.mantissa]);

elseif length(ieee.mantissa) == 13

if any( (ieee.mantissa<'0'|ieee.mantissa>‘f’) | …

(ieee.mantissa>’9’&ieee.mantissa<'a') )
error(['***** hexidecimal mantissa not all 0-9 or a-f: ', ...
ieee.mantissa])
end
Ieee = reshape(exponenet_bin, 4, 3)';
Ieee_dec = bin2dec(Ieee);
ieee_hex = bin2hex(Ieee_dec+1);
x = hex2num([char(ieee_hex), ieee.mantissa]);
else
error(['***** incorrect mantissa length: ', num2str(ieee.mantissa)])
end
end
end
function ieee = fp2ieee(x)
% calculates the IEEE representation of a floating-point number in binary
% input: x = floating-point number
% output: ieee.bits = (structure) IEEE representation, i.e., 64 bits of 0's and 1's
% .split = IEEE representation separated for easier reading
% .sci = scientific notation:
% +- 1.xxxxx... 2^exponent or
% +- 0.xxxxx... 2^(-1023)
% .sign = +1 if x >= 0

% = -1 if x < 0
if ~isfinite(x)
error('***** x is not a finite number')
end
x_hex = num2hex(x);
x_exp_dec = hex2dec(x_hex(1:3));
ieee.sign = +1;
x_sign_char = '+';
x_sign_bin = '0';
if x_exp_dec >= 2^11

x_exp_dec = x_exp_dec – 2^11;

ieee.sign = -1;

x_sign_char = ‘-‘;

x_sign_bin = ‘1’;

end

x_real_exp = x_exp_dec – (2^10 – 1);

x_mantissa = dec2bin(hex2dec(x_hex(4:16)), 52);

ieee.bits = [dec2bin(hex2dec(x_hex(1:3)), 12), x_mantissa];

if x_exp_dec > 0

ieee.sci = [x_sign_char, ‘1.’, x_mantissa];

else

ieee.sci = [x_sign_char, ‘0.’, x_mantissa];

end

ieee.sci = [ieee.sci, ‘ x 2^’, num2str(x_real_exp)];

ieee.split = [x_sign_bin, ‘|’, split(dec2bin(x_exp_dec, 11)), ‘|’, …

split(x_mantissa)];

end

function split_string = split(string)

if length(string) == 0

return

end

if length(string) <= 4
split_string = string;
else
split_string = [split(string(1:end-4)), ' ', string(end-3: ...
end)];
end
end
function floatgui(callbackarg)
%FLOATGUI Show structure of floating point numbers.
% The set of positive model floating point numbers is determined
% by three parameters: t, emin, and emax. It is the set of rational
% numbers of the form x = (1+f)*2^e where f = (integer)/2^t,
% 0 <= f < 1, e = integer, and emin <= e <= emax.
%
% IEEE 754 double precision has t = 52, emin = -1022, emax = 1023.
% Copyright 2014 Cleve Moler
% Copyright 2014 The MathWorks, Inc.
% Initialize parameters
if nargin == 0
t = 3;
emin = -4;
emax = 2;
logscale = 0;
Fpos=[50 300 900 250];
else
t = round(get(findobj('tag','t'),'value'));
emin = round(get(findobj('tag','emin'),'value'));
emax = round(get(findobj('tag','emax'),'value'));
logscale = get(findobj('style','check'),'value');
Fpos=get(gcf,'pos');
end
% Position figure window
shg
clf reset
set(gcf,'pos',Fpos,'name','floatgui', ...
'resize','on','defaultuicontrolunits','normalized',...
'numbertitle','off','menubar','none')
% Generate and plot floating point numbers
f = (0:2^t-1)/2^t;
F = [];
for e = emin:emax
F = [F (1+f)*2^e];
end
for x = F
text(x,0,'|','fontunits','normalized','fontsize',0.3)
end
% Set axes
set(gca,'pos',[.05 .6 .9 .2],'fontunits','normalized','fontsize',0.22)
if logscale
set(gca,'xscale','log')
xmin = 1/2^(-emin+.5);
xmax = 2^(emax+1.5);
else
set(gca,'xscale','linear')
xmin = 0;
xmax = 2^(emax+1);
end
axis([xmin xmax -1 1])
% Set tick marks
fmin = min(F);
fmax = max(F);
xtick = 1;
xticklab = {'1'};
if fmin < 1
xtick = [1/2 xtick];
xticklab = ['1/2' xticklab];
end
if logscale & (fmin < 1/4)
xtick = [1/4 xtick];
xticklab = ['1/4' xticklab];
end
if fmin < 1/2
xtick = [fmin xtick];
xticklab = [['1/' int2str(1/fmin)] xticklab];
end
if 2 < fmax
xtick = [xtick 2];
xticklab = [xticklab '2'];
end
if 4 < fmax
xtick = [xtick 4];
xticklab = [xticklab '4'];
end
if max(xtick) < fmax
xtick = [xtick fmax];
if fmax == round(fmax)
fmaxlab = int2str(fmax);
else
over = 2^(emax+1);
fmaxlab = [int2str(over) '-1/' int2str(1/(over-fmax))];
end
xticklab = [xticklab fmaxlab];
end
set(gca,'xtick',xtick,'xticklabel',xticklab,'xminortick','off','ytick',[])
% Create uicontrols
uicontrol('style','slider','tag','emin','value',emin, ...
'min',-8,'max',0,...
'pos',[0.15 0.26 0.13 0.07],'sliderstep',[1/8 1/8], ...
'callback','floatgui(1)');
uicontrol('style','slider','tag','t','value',t, ...
'min',0,'max',8,...
'pos',[0.435 0.26 0.13 0.07],'sliderstep',[1/8 1/8], ...
'callback','floatgui(1)');
uicontrol('style','slider','tag','emax','value',emax, ...
'min',0,'max',8,...
'pos',[0.72 0.26 0.13 0.07],'sliderstep',[1/8 1/8], ...
'callback','floatgui(1)');
uicontrol('style','text','string',['emin = ' int2str(emin)], ...
'pos',[0.15 0.35 0.13 0.07],'fontunits','normalized','fontsize',0.7)
uicontrol('style','text','string',['t = ' int2str(t)], ...
'pos',[0.435 0.35 0.13 0.07],'fontunits','normalized','fontsize',0.7)
uicontrol('style','text','string',['emax = ' int2str(emax)], ...
'pos',[0.72 0.35 0.13 0.07],'fontunits','normalized','fontsize',0.7)
uicontrol('style','check','string','log scale','value',logscale, ...
'pos',[0.435 0.15 0.13 0.07],'fontunits','normalized','fontsize',0.7, ...
'callback','floatgui(1)');
uicontrol('style','push','pos',[0.88 0.1 0.07 0.07], ...
'fontunits','normalized','fontsize',0.7, ...
'string','close','callback','close(gcf)')
% eps
if fmax > 1

eps = 2^(-t);

text(1,0,’|’,’color’,’r’,’fontunits’,’normalized’,’fontsize’,0.3)

text(1+eps,0,’|’,’color’,’r’,’fontunits’,’normalized’,’fontsize’,0.3)

if eps < 1
text(1.0,1.5,['eps = 1/' int2str(1/eps)], ...
'fontunits','normalized','fontsize',0.3,'fontweight','bold')
else
text(1.0,1.5,'eps = 1', 'fontunits','normalized','fontsize',0.3,'fontweight','bold')
end
end
% Number of numbers
% Exercise:
% How many "floating point" numbers are in the set?
% Complete this statement.
% text(.9*xmax,2,num2str(???)
Catastrophic Cancellation in log(1+x)
Tae Eun Kim
September 23, 2021
This note explains the catastrophic cancellation observed in Problem 4 of Homework 4.
The evaluation of f(x) is severely affected by catastrophic cancellation for small x because of what
is written at the beginning of the problem. Though identical to f(x) mathematically, the function
f1(x) does a better job because of the following reason.
Let x̂ = ̂(1 + x) − 1 = fl((1 +x) − 1), the floating-point representation of the expression (1 +x) − 1.
Note that the subtraction undergoes catastrophic cancellation for small x. Also note that when
log(1 + x) is evaluated in the computer, the input (1 + x) is formed first and then 1 is subtracted
off from it before it is used in the algorithm based on the Taylor series
log ξ = (ξ − 1) −
1
2
(ξ − 1)2 +
1
3
(ξ − 1)3 − · · · .
Therefore the numerical evaluation of f1(x) can be approximated by
f̂1(x) ≈
x̂ − 12x̂
2 + 13x̂
3 − · · ·
x̂
= 1 −
1
2
x̂ +
1
3
x̂2 − · · · ,
which resembles the series expansion used in part (a). It is clear from the right-hand side that
this implementation does not involve subtraction of two nearby numbers for reasonably small x.
However, when x is sufficiently small, (1 + x) is indistinguishable from 1 on the floating-point
number system, in which case
x̂ = ̂(1 + x) − 1 = 0.
We know that it happens when x is smaller than the machine epsilon eps, which is about 2 × 10−16;
in our experiment, it happens for k ≥ 16, in which case both the numerator and the denominator
of f1(x) are evaluated as zeros, resulting in NaN.
One way to avoid the issue with NaN is to use the series expansion obtained in part (a), namely,
f(x) = 1 −
1
2
x +
1
3
x2 − · · · .
Another way is to use the function log1p as suggested in the problem. This function was specifi-
cally designed to avoid catastrophic cancellation occurring in the evaluation of log(1 + x) for small
x by encoding
log(1 + x) = x −
1
2
x2 +
1
3
x3 − · · ·
instead of using the series expansion for log ξ written above.
1
[Content_Types].xml
_rels/.rels
matlab/document.xml
Homework 1 Math 3607, Spring 2021 [Your Name] Table of Contents Problem 1. (LM 2.1--6)
Problem 2. (LM 2.1--8)
Part (a)
(i)
Part (b)
Part (c)
Problem 3. (LM 2.1--14)
Part (a)
Problem 4. (Temperature Conversion)
Problem 5. (Oblate Spheroid) Problem 1. (LM 2.1--6) Your answer here. % Code goes here Problem 2. (LM 2.1--8) Part (a) (i) (ii) (vi) Part (b) Part (c) Problem 3. (LM 2.1--14) Part (a) Problem 4. (Temperature Conversion) Problem 5. (Oblate Spheroid)
matlab/output.xml
manual code ready 0.4
metadata/coreProperties.xml
2021-01-13T21:38:03Z 2021-01-13T21:38:03Z
metadata/mwcoreProperties.xml
application/vnd.mathworks.matlab.code MATLAB Code R2020b
metadata/mwcorePropertiesExtension.xml
9.9.0.1444674 93b378b9-5b5d-4623-b00d-0e154b3a183c
metadata/mwcorePropertiesReleaseInfo.xml
9.9.0.1524771
R2020b
Update 2
Nov 03 2020
2207788044

#### Why Choose Us

- 100% non-plagiarized Papers
- 24/7 /365 Service Available
- Affordable Prices
- Any Paper, Urgency, and Subject
- Will complete your papers in 6 hours
- On-time Delivery
- Money-back and Privacy guarantees
- Unlimited Amendments upon request
- Satisfaction guarantee

#### How it Works

- Click on the “Place Order” tab at the top menu or “Order Now” icon at the bottom and a new page will appear with an order form to be filled.
- Fill in your paper’s requirements in the "
**PAPER DETAILS**" section. - Fill in your paper’s academic level, deadline, and the required number of pages from the drop-down menus.
- Click “
**CREATE ACCOUNT & SIGN IN**” to enter your registration details and get an account with us for record-keeping and then, click on “PROCEED TO CHECKOUT” at the bottom of the page. - From there, the payment sections will show, follow the guided payment process and your order will be available for our writing team to work on it.