{"id":4813,"date":"2013-02-08T17:26:55","date_gmt":"2013-02-08T16:26:55","guid":{"rendered":"http:\/\/www.walkingrandomly.com\/?p=4813"},"modified":"2013-02-11T10:24:49","modified_gmt":"2013-02-11T09:24:49","slug":"fun-with-mathematicas-reduce-and-floating-point-arithmetic","status":"publish","type":"post","link":"https:\/\/walkingrandomly.com\/?p=4813","title":{"rendered":"Fun with Mathematica’s Reduce and floating point arithmetic"},"content":{"rendered":"

I was at a seminar yesterday where we were playing with Mathematica and wanted to solve the following equation<\/p>\n

1.0609^t == 1.5<\/pre>\n
You punch this into Mathematica as follows:<\/p>\n
Solve[1.0609^t == 1.5]<\/pre>\n
Mathematica returns the following<\/p>\n
During evaluation of In[1]:= Solve::ifun: Inverse functions are being used by Solve, \r\nso some solutions may not be found; use Reduce for complete solution information. >>\r\n\r\nOut[1]= {{t -> 6.85862}}<\/pre>\n
I have got the solution I expect but Mathematica suggests that I’ll get more information if I use Reduce. So, lets do that.<\/p>\n
In[2]:= Reduce[1.0609^t == 1.5, t]\r\nDuring evaluation of In[2]:= Reduce::ratnz: Reduce was unable to solve the system with inexact \r\ncoefficients. The answer was obtained by solving a corresponding exact system and numericizing the result. >>\r\nOut[20]= C[1] \\[Element] Integers && \r\n  t == -16.9154 (-0.405465 + (0. + 6.28319 I) C[1])<\/pre>\n
Looks complex and a little complicated! To understand why complex numbers have appeared in the mix you need to know that Mathematica always considers variables to be complex unless you tell it otherwise. So, it has given you the infinite number of complex values of t that would satisfy the original equation. No problem, let’s just tell Mathematica that we are only interested in real values of p.<\/p>\n
In[3]:= Reduce[1.0609^t == 1.5, t, Reals]\r\n\r\nDuring evaluation of In[3]:= Reduce::ratnz: Reduce was unable to solve the system with inexact \r\ncoefficients. The answer was obtained by solving a corresponding exact system and numericizing the result. >>\r\n\r\nOut[3]= t == 6.85862<\/pre>\n
Again, I get the solution I expect. However, Mathematica still feels that it is necessary to give me a warning message. Time was ticking on so I posed the question on Mathematica Stack Exchange and we moved on.<\/p>\n
http:\/\/mathematica.stackexchange.com\/questions\/19235\/why-does-mathematica-struggle-with-solving-this-equation<\/a><\/p>\n
At the time of writing, the lead answer says that \u2018Mathematica is not “struggling” with your equation. The message is simply FYI — to tell you that, for this equation, it prefers to work with exact quantities rather than inexact quantities (reals)\u2019<\/p>\n
I’d accept this as an explanation except for one thing; the message say\u2019s that it is UNABLE to solve the equation as I originally posed it and so needs to proceed by solving a corresponding exact system. That implies that it has struggled a great deal, given up and tried something else.<\/p>\n
This would come as a surprise to anyone with a calculator who would simply do the following manipulations<\/p>\n
1.0609^t == 1.5<\/p>\n
Log[1.0609^t] == Log[1.5]<\/p>\n
T*Log[1.0609] == Log[1.5]<\/p>\n
T= Log[1.5]\/Log[1.0609]<\/p>\n
Mathematica evalutes this as<\/p>\n
T=6.8586187084788275<\/p>\n
This appears to solve the equation exactly.\u00a0 Plug it back into Mathematica (or my calculator) to get<\/p>\n
In[4]:= 1.0609^(6.8586187084788275)\r\n\r\nOut[4]= 1.5<\/pre>\n
I had no trouble dealing with inexact quantities, and I didn\u2019t need to \u2018solve a corresponding exact system and numericize the result\u2019.\u00a0 This appears to be a simple problem. So, why does Mathematica bang on about it so much?<\/p>\n
Over to MATLAB for a while<\/strong><\/p>\n
Mathematica is complaining that we have asked it to work with inexact quantities.\u00a0 How could this be? 1.0609, 6.8586187084788275 and 1.5 look pretty exact to me! However, it turns out that as far as the computer is concerned some of these numbers are far from exact.<\/p>\n
When you input a number such as 1.0609 into Mathematica, it considers it to be a double precision<\/a> number and 1.0609 cannot be exactly represented as such.\u00a0 The closest Mathematica, or indeed any numerical system that uses 64bit IEEE arithmetic, can get is 4777868844677359\/4503599627370496 which evaluates to 1.0608999999999999541699935434735380113124847412109375.  I wondered if this is why Mathematica was complaining about my problem.<\/p>\n
At this point, I switch tools and use MATLAB for a while in order to investigate further.\u00a0 I do this for no reason other than I know the related functions in MATLAB a little better.\u00a0 MATLAB’s sym command can give us the rational number that exactly represents what is actually stored in memory (Thanks to Nick Higham for pointing this out to me).<\/p>\n
>> sym(1.0609,'f')\r\n\r\nans =\r\n4777868844677359\/4503599627370496<\/pre>\n
We can then evaluate this fraction to whatever precision we like using the vpa command:<\/p>\n
>> vpa(sym(1.0609,'f'),100)\r\nans =\r\n1.0608999999999999541699935434735380113124847412109375\r\n\r\n>> vpa(sym(1.5,'f'),100)\r\nans =\r\n1.5\r\n\r\n>> vpa(sym( 6.8586187084788275,'f'),100)\r\nans =\r\n6.8586187084788274859192824806086719036102294921875<\/pre>\n
So, 1.5 can be represented exactly in 64bit double precision arithmetic but 1.0609 and 6.8586187075 cannot.\u00a0 Mathematica is unhappy with this state of affairs and so chooses to solve an exact problem instead.\u00a0 I guess if I am working in a system that cannot even represent the numbers in my problem (e.g. 1.0609) how can I expect to solve it?<\/p>\n
Which Exact Problem?<\/strong>
\nSo, which exact equation might Reduce be choosing to solve?\u00a0 It could solve the equation that I mean:<\/p>\n
(10609\/10000)^t == 1500\/1000<\/pre>\n
which does have an exact solution and so Reduce can find it.<\/p>\n
(Log[2] - Log[3])\/(2*(2*Log[2] + 2*Log[5] - Log[103]))<\/pre>\n
Evaluating this gives 6.858618708478698:<\/p>\n
(Log[2] - Log[3])\/(2*(2*Log[2] + 2*Log[5] - Log[103])) \/\/ N \/\/ FullForm\r\n\r\n6.858618708478698`<\/pre>\n
Alternatively, Mathematica could convert the double precision number 1.0609 to the fraction that exactly represents what’s actually in memory and solve that.<\/p>\n
(4777868844677359\/4503599627370496)^t == 1500\/1000<\/pre>\n
This also has an exact solution:<\/p>\n
(Log[2] - Log[3])\/(52*Log[2] - Log[643] - Log[2833] - Log[18251] - Log[143711])<\/pre>\n
which evaluates to\u00a06.858618708478904:<\/p>\n
(Log[2] - Log[3])\/(52*Log[2] - Log[643] - Log[2833] - Log[18251] - Log[143711]) \/\/ N \/\/ FullForm\r\n\r\n6.858618708478904`<\/pre>\n
Let’s take a look at the exact number Reduce is giving:<\/p>\n
Quiet@Reduce[1.0609^t == 1.5, t, Reals] \/\/ N \/\/ FullForm\r\n\r\nEqual[t, 6.858618708478698`]<\/pre>\n
So, it looks like Mathematica is solving the equation I meant to solve and evaluating this solution it at the end.<\/p>\n
Summary of solutions<\/strong>
\nHere I summarize the solutions I’ve found for this problem:<\/p>\n
    \n
  • 6.8586187084788275 – Pen,Paper+Mathematica for final evaluation.<\/li>\n
  • 6.858618708478698 – Mathematica solving the exact problem I mean and evaluating to double precision at the end.<\/li>\n
  • 6.858618708478904 – Mathematica solving the exact problem derived from what I really asked it and evaluating at the end.<\/li>\n<\/ul>\n
    What I find fun is that my simple minded pen and paper solution seems to satisfy the original equation better than the solutions arrived at by more complicated means.\u00a0 Using MATLAB’s vpa again I plug in the three solutions above, evaluate to 50 digits and see which one is closest to 1.5:<\/p>\n
    >> vpa('1.5 - 1.0609^6.8586187084788275',50)\r\nans =\r\n-0.00000000000000045535780896732093552784442911868195148156712149736\r\n>> vpa('1.5 - 1.0609^6.858618708478698',50)\r\nans =\r\n0.000000000000011028236861872639054542278886208515813177349441555\r\n>> vpa('1.5 - 1.0609^6.858618708478904',50)\r\nans =\r\n-0.0000000000000072391029234017787617507985059241243968693347431197<\/pre>\n
    So, ‘my’ solution is better. Of course, this advantage goes away if I evaluate (Log[2] – Log[3])\/(2*(2*Log[2] + 2*Log[5] – Log[103])) to 50 decimal places and plug that in.<\/p>\n
    \r\n>> sol=vpa('(log(2) - log(3))\/(2*(2*log(2) + 2*log(5) - log(103)))',50)\r\nsol =\r\n6.858618708478822364949699852597847078441119153527\r\n>> vpa(1.5 - 1.0609^sol',50)\r\nans =\r\n0.0000000000000000000000000000000000000014693679385278593849609206715278070972733319459651\r\n<\/pre>\n","protected":false},"excerpt":{"rendered":"
    I was at a seminar yesterday where we were playing with Mathematica and wanted to solve the following equation 1.0609^t == 1.5 You punch this into Mathematica as follows: Solve[1.0609^t == 1.5] Mathematica returns the following During evaluation of In[1]:= Solve::ifun: Inverse functions are being used by Solve, so some solutions may not be found; […]<\/p>\n","protected":false},"author":3,"featured_media":0,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"jetpack_post_was_ever_published":false,"footnotes":"","jetpack_publicize_message":"","jetpack_publicize_feature_enabled":true,"jetpack_social_post_already_shared":false,"jetpack_social_options":{"image_generator_settings":{"template":"highway","default_image_id":0,"font":"","enabled":false},"version":2}},"categories":[8,11],"tags":[],"class_list":["post-4813","post","type-post","status-publish","format-standard","hentry","category-mathematica","category-matlab"],"jetpack_publicize_connections":[],"jetpack_featured_media_url":"","jetpack_shortlink":"https:\/\/wp.me\/p3swhs-1fD","jetpack_sharing_enabled":true,"_links":{"self":[{"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=\/wp\/v2\/posts\/4813","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=\/wp\/v2\/users\/3"}],"replies":[{"embeddable":true,"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=%2Fwp%2Fv2%2Fcomments&post=4813"}],"version-history":[{"count":21,"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=\/wp\/v2\/posts\/4813\/revisions"}],"predecessor-version":[{"id":4829,"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=\/wp\/v2\/posts\/4813\/revisions\/4829"}],"wp:attachment":[{"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=%2Fwp%2Fv2%2Fmedia&parent=4813"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=%2Fwp%2Fv2%2Fcategories&post=4813"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=%2Fwp%2Fv2%2Ftags&post=4813"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}