problems in accuracy assessment

1 次查看(过去 30 天)
Dear friends I tried continued fractiona approximation to pi, and got the following rational approximation
1244969988745789040106366155256015114976454553337906033890313/396286255419907262612262286579004636343912658949984351074158
but if I try to find the difference of it with pi, I failed, since
>> sym(1244969988745789040106366155256015114976454553337906033890313/396286255419907262612262286579004636343912658949984351074158)
I got pi directly. How can I assess the difference? Thanks
  2 个评论
John D'Errico
John D'Errico 2014-12-24
Please, learn to use the code button to allow your code to be displayed properly. Otherwise, it is unreadable.
Dingyu Xue
Dingyu Xue 2014-12-24
Thanks John, but it seems I cannot use CODE
Dear friends I tried continued fractiona approximation to pi, and got the following rational approximation
1244969988745789040106366155256015114976454553337906033890313/396286255419907262612262286579004636343912658949984351074158
but if I try to find the difference of it with pi, I failed, since
sym(1244969988745789040106366155256015114976454553337906033890313/396286255419907262612262286579004636343912658949984351074158)
I got pi directly. How can I assess the difference?
Thanks

请先登录,再进行评论。

采纳的回答

John D'Errico
John D'Errico 2014-12-24
编辑:John D'Errico 2014-12-24
Using my HPF toolbox, it looks to be correct to about 120 significant digits.
double(log10(abs(hpf('1244969988745789040106366155256015114976454553337906033890313',1000)/hpf('396286255419907262612262286579004636343912658949984351074158',1000) - hpf('pi',1000))))
ans =
-119.53
Looking back at what you did, this is wrong. Why?
sym(1244969988745789040106366155256015114976454553337906033890313/396286255419907262612262286579004636343912658949984351074158)
ans =
pi
It is wrong because MATLAB will first compute the ratio as a DOUBLE precision number, and only then pass it to sym. sym then recognizes the ratio as a double precision approximation to pi, so it returns that as the result.
If you insist on doing this with the symbolic toolbox, we might try...
log10(abs(vpa(sym('1244969988745789040106366155256015114976454553337906033890313/396286255419907262612262286579004636343912658949984351074158 - pi'),200)))
ans =
-119.52655990941722562288816970381
Of course, if you want a better approximation, say accurate to 250 digits, HPF can provide it too.
[N,D] = rat(hpf('pi',500),hpf('1e-250'))
N =
165213461336854892880967397129207248425612593828416835132438559586109735260365863184040108685528782183687397179966550074377225
D =
52589078074164381557662496030246690228630477463437042437123678996605078304650603819026850043803035007564529614801553827009746
If you need convincing that this ratio is indeed accurate to 250 digits for pi, try this:
double(log10(abs(N/D - hpf('pi',1000))))
ans =
-250.05
  2 个评论
Dingyu Xue
Dingyu Xue 2014-12-24
Thanks John
Excellent!
Actually I did got the rational number with continued fraction, contfrac in MUPAD, with 120 levels. I tried your method without HPF Toolbox, and made continued fraction with 500 levels, the accuracy is about 500 digits. Is there any hidden rule in the continued fractions? The number of levels is almost the same as the number of correct digits?
John D'Errico
John D'Errico 2014-12-24
Well, you cannot have more correct digits in a rational fraction approximation to a number than are in the original number. So if you are working with 120 digit approximations to pi, you can get roughly 120 correct digits. In the HPF tool, you can provide a tolerance for the approximation. If none is provided, it will assume a tolerance.
And no, in terms of rules, if you are asking if there is any rule for the approximation itself, I don't think that there can exist any simple rule for a rational fraction approximation to a number like pi, although I do not know if that has ever been proven. I'm sure Google would tell me though. Here is the continued fraction I generate for a 100 digit approximation to pi.
3 + 1/(7 + 1/(16 + 1/(-294 + 1/(3 + 1/(-4 + 1/(5 + 1/(-15 + 1/(-3 + 1/(2 + 1/(2 + 1/(2 + 1/(2 + 1/(3 + 1/(-85 + 1/(-3 + 1/(2 + 1/(15 + 1/(3 + 1/(14 + 1/(-5 + 1/(-2 + 1/(-6 + 1/(-6 + 1/(-100 + 1/(3 + 1/(2 + 1/(6 + 1/(3 + 1/(6 + 1/(-2 + 1/(-6 + 1/(-9 + 1/(9 + 1/(-3 + 1/(-3 + 1/(-8 + 1/(4 + 1/(-2 + 1/(-13 + 1/(3 + 1/(-5 + 1/(2 + 1/(9 + 1/(-2 + 1/(-3 + 1/(8 + 1/(-2 + 1/(-5 + 1/(-2 + 1/(-2 + 1/(-4 + 1/(3 + 1/(4 + 1/(4 + 1/(17 + 1/(-162 + 1/(-46 + 1/(24 + 1/(-3 + 1/(-3 + 1/(6 + 1/(-3 + 1/(-25 + 1/(4 + 1/(-5 + 1/(4 + 1/(-2)))))))))))))))))))))))))))))))))))))))))))))))))))))))))))))))))))

请先登录,再进行评论。

更多回答(0 个)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by