matlab 6.5和2010a里的vpa以及pi
最初是matlab中文論壇的網(wǎng)友他尋歡發(fā)現(xiàn)的,，在matlab 7.0, 7.1, 7.4和2010a下運行vpa(pi,100)分別得到如下結(jié)果
3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068
3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068
3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068
3.141592653589793115997963468544185161590576171875
2010a的結(jié)果跟前幾個版本在小數(shù)點后15位就不一樣了
因為手頭沒有7.1-7.4的版本，倒是為了用老版的netcdf工具箱還留著個6.5,，就在6.5上試了一下,，確實跟2010a的結(jié)果不一樣。很有意思的事情,，翻了一下vpa的help,，兩個版本還是有點區(qū)別的。
這是2010a的,，6.5的沒有紅色的那部分
VPA Variable precision arithmetic.
R = VPA(S) numerically evaluates each element of the double matrix
S using variable precision floating point arithmetic with D decimal
digit accuracy, where D is the current setting of DIGITS.
The resulting R is a SYM.
VPA(S,D) uses D digits, instead of the current setting of DIGITS.
D is an integer or the SYM representation of a number.
It is important to avoid the evaluation of an expression using double
precision floating point arithmetic before it is passed to VPA.
For example,
phi = vpa((1+sqrt(5))/2)
first computes a 16-digit approximation to the golden ratio, then
converts that approximation to one with d digits, where d is the current
setting of DIGITS. To get full precision, use unevaluated string or
symbolic arguments,
phi = vpa('(1+sqrt(5))/2')
or
s = sym('sqrt(5)')
phi = vpa((1+s)/2);
Additional examples:
vpa(pi,780) shows six consecutive 9's near digit 770 in the
decimal expansion of pi.
vpa(hilb(2),5) returns
[ 1., .50000]
[.50000, .33333]
See also double, digits.
Overloaded methods:
sym/vpa
從這里來看實際上是vpa(pi,100)這個用法不是很標(biāo)準(zhǔn),，應(yīng)該是vpa('pi',100)，這樣算出來在6.5和2010a下都是
3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068
到這一步上面的問題基本算結(jié)束了,，但是為什么2010a下vpa(pi,100)的結(jié)果是不正確的,，問題出在pi上還是vpa上？
首先想到的是用其他的常數(shù)對比一下,，比如說自然對數(shù)的底e
6.5的
>> vpa(exp(1),100)
ans =
2.718281828459045534884808148490265011787414550781250000000000000000000000000000000000000000000000000
>> vpa('exp(1)',100)
ans =
2.718281828459045235360287471352662497757247093699959574966967627724076630353547594571382178525166427
2010a的
>> vpa(exp(1),100)
ans =
2.71828182845904553488480814849026501178741455078125
>> vpa('exp(1)',100)
ans =
2.718281828459045235360287471352662497757247093699959574966967627724076630353547594571382178525166427
vpa('exp(1)',100)兩個是一樣的,，vpa(exp(1),100)則只差了一串0。后來覺得這個里面還有個常數(shù)1,，不放心,，又試了一下vpa(eps)，兩個結(jié)果仍然一樣,，只是顯示格式略有區(qū)別,，而vpa('eps')的結(jié)果則是
ans =
eps
被系統(tǒng)按照普通的字符串處理了...
看起來是pi的問題，但pi是內(nèi)置函數(shù),，兩個版本pi的help連標(biāo)點符號都一樣,。后來不死心,，去看vpa的源代碼，頗有柳暗花明的感覺
6.5的,，真簡練
if nargin < 2
d = digits;
end;
r = vpa(sym(s),d);
2010a的
eng = symengine;
if strcmp(eng.kind,'maple')
if nargin == 1
r = mapleengine('vpa',s);
else
r = mapleengine('vpa',s,d);
end
if isa(r,'maplesym')
r = sym(r);
end
else
if nargin < 2
d = digits;
end;
if ischar(s)
ss = evalin(symengine,s);
else
ss = sym(s,'f');
end
r = vpa(ss,d);
end
外層if結(jié)構(gòu)的第一塊不用看,，紅色的部分是關(guān)鍵，當(dāng)輸入是char的時候兩個基本的等價(只不過似乎2010a支持多種符號運算引擎,，所以多了evalin那行,，不過默認安裝好像仍然只有一個可以選)，但對于不是字符的情況,，6.5下用的是
sym(s)
而2010a下則是
sym(s,'f')
這是sym的help里相關(guān)的說明
S = sym(A,flag) converts a numeric scalar or matrix to symbolic form.
The technique for converting floating point numbers is specified by
the optional second argument, which may be 'f', 'r', 'e' or 'd'.
The default is 'r'.
'f' stands for 'floating point'. All values are transformed from
double precision to exact numeric values N*2^e for integers N and e.
'r' stands for 'rational'. Floating point numbers obtained by
evaluating expressions of the form p/q, p*pi/q, sqrt(p), 2^q and 10^q
for modest sized integers p and q are converted to the corresponding
symbolic form. This effectively compensates for the roundoff error
involved in the original evaluation, but may not represent the floating
point value precisely. If no simple rational approximation can be
found, the 'f' form is used.
sym(s)相當(dāng)于sym(s,'r')這是一種有理數(shù)的近似,，從help來看，對于可以用兩個不太大的p和q以p/q, p*pi/q, sqrt(p), 2^q, 10^q的形式表示的數(shù),，r要比f精確,，如果找不到這樣的形式,，則r等價于f,，至此，問題全部解決,，只是不明白2010a為什么舍棄了精度較高的r而選擇f,。想要改回來的朋友可以自行修改vpa.m，命令行下open就行了,。
個人對符號運算不算熟,，有疏漏之處歡迎指正。
最后留個小尾巴,，其實6.5和2010a的sym函數(shù)也是不一樣的
當(dāng)flag是r的時候沒什么區(qū)別,，f的時候一個用的是2的指數(shù)的形式，一個是大整數(shù)的分?jǐn)?shù)
6.5
>> sym(1/3,'f')
ans =
'1.5555555555555'*2^(-2)
2010a
>> sym(1/3,'f')
ans =
6004799503160661/18014398509481984
有興趣可以去研究下~
