MatLab - variable precision arithmetic

Question

I have a brief question regarding the vpa command one may use to evaluate symbolic expressions in MatLab.

My textbook says the following:

"You need to be careful when you use functions such as sqrt on numbers, which by default result in a double-precision floating-point number. You need to pass such input to vpa as a symbolic string for correct evaluation: vpa('sqrt(5)/pi')."

I don't quite understand the jargon here. Why is it that for most inputs I get the exact same answer whether I type vpa(input) or vpa('input'), but not for square roots? For instance, if I type vpa(sin(pi/4)) or vpa('sin(pi/4)'), I get the exact same answers, but if I type the give problem above as vpa(sqrt(5)/pi), I do not get the same answer as when I type vpa('sqrt(5)/pi').

If someone could explain this in a little more detail than what my book does above, I would be very grateful!

Futurity

you might be interested to read this question: In MATLAB, are variables REALLY double-precision by default? – Shannan

Answer 1

I'm no MatLab expert, but without quotes, you're passing the result of sqrt(5)/pi into vpa():

  vpa(sqrt(5)/pi)
= vpa(0.7117625434171772)

With quotes, you're passing the expression sqrt(5)/pi (unevaluated and in exact form) into vpa() and then telling MatLab to compute sqrt(5)/pi with variable precision.

Portemonnaie

Thanks! That makes sense. However, why is it that when I enter simply sin(pi/4) into vpa, that it doesn't matter whether or not I use quotes? – Futurity

It should matter. As @Ben Volgt pointed out, sin(pi / 4) = sin(45 deg) = sqrt(2) / 2, which is irrational and won't be exactly approximated by double precision floats. – Portemonnaie

After reading the docs it sounds like vpa returns a double-precision result, in which case only the precision of intermediate values is improved. – Thorlay

Great. Much appreciated! I think I get this now :) – Futurity

Answer 2

Never assume that a number like vpa(sin(pi/4)) is exact to the full precision, because MATLAB will generally compute the number inside the call to vpa using floating point arithmetic, so only accurate to about 16 digits.

However, it appears that it is correct here. For example, we know that

sin(pi/4) == sqrt(2)/2

Lets test that result. I'll use 100 digits of precision, comparing both vpa and my own HPF tools.

>> vpa(sin(pi/4),100)
ans =
0.7071067811865475244008443621048490392848359376884740365883398689953662392310535194251937671638207864

>> vpa(sqrt(sym(2))/2,100)
ans =
0.7071067811865475244008443621048490392848359376884740365883398689953662392310535194251937671638207864

>> sqrt(hpf(2,100))/2
ans =
0.7071067811865475244008443621048490392848359376884740365883398689953662392310535194251937671638207864

>> sin(hpf('pi',100)/4)
ans =
0.7071067811865475244008443621048490392848359376884740365883398689953662392310535194251937671638207864

So, my guess is the parser has recognized the input as something the symbolic toolbox can compute more accurately. As I said before, be careful though. What is sin(pi/12)?

>> vpa(sin(pi/12),100)
ans =
0.25881904510252073947640383266843855381011962890625

>> vpa('sin(pi/12)',100)
ans =
0.2588190451025207623488988376240483283490689013199305138140032073150569747488019969223679746942496655

>> vpa(sin(sym('pi')/12),100)
ans =
0.2588190451025207623488988376240483283490689013199305138140032073150569747488019969223679746942496655

>> sin(hpf('pi',100)/12)
ans =
0.2588190451025207623488988376240483283490689013199305138140032073150569747488019969223679746942496655

See that in the first case, the parser did not save us. In the others, I forced MATLAB to compute the correct value. In fact, a bit of effort will give us the value for sin(pi/12), as sqrt(2)*(sqrt(3) - 1)/4.

>> DefaultNumberOfDigits 100
>> (sqrt(hpf(3)) - 1)*sqrt(hpf(2))/4
ans =
0.2588190451025207623488988376240483283490689013199305138140032073150569747488019969223679746942496655

The point is, do not trust the parser to save you here.

Edit: As a test of Amro's comment, I respectfully state that MATLAB IS doing something of interest here. See that vpa is able to return the correct first 100 digits of pi, even when passed pi as a double precision number. Since pi (as a double) is not correct past about the 16th decimal digit, there is something fishy going on.

>> vpa(pi,100)
ans =
3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068

>> vpa('pi',100)
ans =
3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068

vpa('pi',100) - vpa(pi,100)
ans =
0.0

As a test of that fact, lets look at what HPF finds. HPF actually takes the IEEE 754 value, as stored in a double, then converts it to an HPF number.

>> hpf(pi,100)
ans =
3.141592653589793115997963468544185161590576171875

>> hpf('pi',100)
ans =
3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068

>> hpf('pi',100) - hpf(pi,100)
ans =
0.0000000000000001224646799147353177226065932275001058209749445923078164062862089986280348253421170679821480800000000

So clearly, MATLAB is able to recognize pi as something more than just the double precision value it will be passed in as.

Edit2:

In fact, a bit of play tells me what is happening here. VPA is the tricky one, not the parser. Consider the fraction 7/13. If we build it as a double, then print out the floating point value stored in its full glory, we see it is not truly exact. This is as expected.

>> sprintf('%.100f',7/13)
ans =
0.5384615384615384359179302009579259902238845825195312500000000000000000000000000000000000000000000000

7/13 is a repeating decimal value. Here are the correct digits:

>> vpa('7/13',100)
ans =
0.5384615384615384615384615384615384615384615384615384615384615384615384615384615384615384615384615385

Now, suppose we try to create the same number. Here I'll pass 7/13 in asa double, but I'll make a mistake in the bottom decimal digits

>> sprintf('%.100f',0.538461538461538461777777777)
ans =
0.5384615384615384359179302009579259902238845825195312500000000000000000000000000000000000000000000000

Here we see that vpa catches and corrects the 'error' I've made, recognizing that what I passed in is actually identically the same value as when I passed in 7/13.

>> vpa(0.538461538461538461777777777,100)
ans =
0.5384615384615384615384615384615384615384615384615384615384615384615384615384615384615384615384615385

Of course, if I pass in the value as a string, then vpa gets it wrong.

>> vpa('0.538461538461538461777777777',100)
ans =
0.538461538461538461777777777

This explains why vpa is able to catch and correctly compute vpa(sin(pi/4),100), out to the full precision asked. sin(pi/4) is computed as a double, but then vpa sees it as a number that is the same as a double precision version of sqrt(2)/2.

Be careful of course. For example, vpa is not smart enough to catch this simple shift of pi.

>> vpa(pi + 1,100)
ans =
4.141592653589793115997963468544185161590576171875

>> vpa(pi,100)
ans =
3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068

Ghetto

just to clear any confusion, the parser does not treat the input in vpa(sin(pi/4)) any different from calling any other function myFcn(sin(pi/4))... MATLAB first evaluates the expression in double precision before passing the result to the function. With finite precision, the result is rounded to the nearest decimal value it can represent, which VPA treats as an exact value. At that point it can calculate as much decimal points as you wish... I believe @AndrewJanke gave a good explanation in the question I linked to. – Shannan

Thank you very much! I truly appreciate your detailed answer. – Futurity

@woodchips: it seems you are also correct. When you call vpa with a numeric expression, it internally calls vpa(sym(expr)). Now when sym is called with a number, it tries to convert the number to "rational" form in order to compensate for round-off error (flag = 'r' argument), which explains why VPA "corrected" the error as you showed. If you want to get the same results as your HPF tool, I guess it would be: vpa(sym('pi')-sym(pi,'d'), 100) – Shannan

@woodchips: No wonder I've missed this detail, looking at the documentation in previous versions, it wasn't obvious that such conversion was taking place, and its not until R2012a that the doc mentions this fact clearly... +1 for a really enlightening answer – Shannan

Answer 3

If you're getting the exact same answer, you didn't need variable precision arithmetic to begin with.

However, sin(pi/4) should be exactly sqrt(2)/2, which is irrational. You should not get exactly the same answer from different precision. Perhaps you should check how you're displaying (and rounding off) the result.

Thorlay

Thanks. I get this now :). Appreciate it! – Futurity

Answer 4

The latest doc on numeric to symbolic conversion has the answer.

sym tries to correct the round-off error in floating-point inputs to return the exact symbolic form. Specifically, sym corrects round-off error in numeric inputs that match the forms p/q, pπ/q, (p/q)^(1/2), 2^q, and 10^q, where p and q are modest-sized integers.

Thus, sin(pi/4) is 2^(1/2)/2 or (1/2)^(1/2) so the vpa command recognizes it. However, sqrt(5)/pi is not a recognized input form according to the doc.

Mckelvey