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.