Get notified about new tutorials
RECEIVE NEW TUTORIALS

First 15 Minutes Free

310 sessions given

since Jul 09, 2014

since Jul 09, 2014

Likelihood of Reply:
95%

Response Time:
within an hour

Ray Phan

Jan 30, 2015

<p>First and foremost, there is no such thing as a <code>do</code> keyword in MATLAB, so eliminate that from your code. Also, don't use <a href="http://www.mathworks.com/help/matlab/ref/eps.html?refresh=true" rel="nofollow"><code>eps</code></a> as an actual variable. This is a pre-defined function in MATLAB that calculates <a href="http://en.wikipedia.org/wiki/Machine_epsilon" rel="nofollow">machine epsilon</a>, which is also what you are trying to calculate. By creating a variable called <code>eps</code>, you would shadow over the actual function, and so any other functions in MATLAB that require its use will behave unexpectedly, and that's not what you want. </p>
<p>Use something else instead, like <code>macheps</code>. Also, your algorithm is slightly incorrect. You need to check for <code>1.0 + (macheps/2)</code> in your <code>while</code> loop, not <code>1.0 + macheps</code>.</p>
<p>In other words, do this:</p>
<pre><code>macheps = 1;
while 1.0 + (macheps/2) > 1.0
macheps = macheps / 2;
end
</code></pre>
<p>This should give you <code>2.22 x 10^{-16}</code>, which agrees with MATLAB if you type in <code>eps</code> in the command prompt. To double-check:</p>
<pre><code>>> format long
>> macheps
macheps =
2.220446049250313e-16
>> eps
ans =
2.220446049250313e-16
</code></pre>
<hr>
<h1>Bonus</h1>
<p>In case you didn't know, machine epsilon is the upper bound on the relative error due to floating point arithmetic. In other words, this would be the maximum difference expected between a true floating point number and one that is calculated on a computer due to the finite number of bits used to store a floating point number. </p>
<p>If you recall, floating numbers inevitably are represented as binary bits on your computer (or pretty much anything digital). In terms of the <a href="http://en.wikipedia.org/wiki/IEEE_floating_point" rel="nofollow">IEEE 754 floating point standard</a>, MATLAB assumes that all numerical values are of type <code>double</code>, which represents floating point numbers as 64-bits. You can obviously override this behaviour by explicitly casting to another type. With the IEEE 754 floating point standard, for <code>double</code> precision type numbers, there are 52 bits that represent the <strong>fractional</strong> part of the number. </p>
<p>Here's a nice diagram of what I am talking about: </p>
<p><img src="http://upload.wikimedia.org/wikipedia/commons/thumb/a/a9/IEEE_754_Double_Floating_Point_Format.svg/618px-IEEE_754_Double_Floating_Point_Format.svg.png" alt=""></p>
<p><sup>Source: <a href="http://en.wikipedia.org/wiki/Double-precision_floating-point_format" rel="nofollow">Wikipedia</a></sup></p>
<p>You see that there is one bit reserved for the sign of the number, 11 bits that are reserved for exponent base and finally, 52 bits are reserved for the fractional part. This in total adds up to 64 bits. The fractional part is a collection or summation of numbers of base 2, with negative exponents starting from -1 down to -52. The MSB of the floating point number starts with<code>2^{-1}</code>, all the way down to <code>2^{-52}</code> as the LSB. Essentially, machine epsilon calculates the maximum resolution difference for an increase of 1 bit in binary between two numbers, given that they have the <strong>same</strong> sign and the <strong>same</strong> exponent base. Technically speaking, machine epsilon actually equals to <code>2^{-52}</code> as this is the maximum resolution of a single bit in floating point, given those conditions that I talked about earlier.</p>
<p>If you actually look at the code above closely, the division by 2 is <strong>bit-shifting</strong> your number to the right by one position at each iteration, starting at the whole value of 1, or <code>2^{0}</code>, and we take this number and add this to 1. We keep bit-shifting, and seeing what the value is equal to by adding this bit-shifted value by 1, and we go up until the point where when we bit shift to the right, a change is no longer registered. If you bit-shift any more to the right, the value will become <strong>0</strong> due to underflow, and so <code>1.0 + 0.0 = 1.0</code>, and this is no longer <code>> 1.0</code>, which is what the <code>while</code> loop is checking. </p>
<p>Once the <code>while</code> loop quits, it is this threshold that defines machine epsilon. If you're curious, if you punch in <code>2^{-52}</code> in the command prompt, you will get what <code>eps</code> is equal to:</p>
<pre><code>>> 2^-52
ans =
2.220446049250313e-16
</code></pre>
<p>This makes sense as you are shifting one bit to the right 52 times, and the point before the loop stops would be at its LSB, which is <code>2^{-52}</code>. For the sake of being complete, if you were to place a counter inside your <code>while</code> loop, and count up how many times the <code>while</code> loop executes, it will execute exactly 52 times, representing 52 bit shifts to the right:</p>
<pre><code>macheps = 1;
count = 0;
while 1.0 + (macheps/2) > 1.0
macheps = macheps / 2;
count = count + 1;
end
>> count
count =
52
</code></pre>
<p>This tip was originally posted on <a href="http://stackoverflow.com/questions/27511349/MATLAB%20-%20Calculating%20machine%20epsilon/27511886">Stack Overflow</a>.</p>