如何用matlab做高精度计算？【第三辑】(完)
如何用matlab做高精度计算？【第三辑】(完)
高精度计算是一种程序设计的算法。由于中央处理器的字长限制，如32位CPU中一个整数最大只能取值4,294,967,295（=2^32-1），因此在超范围数值计算中，往往要采用模拟手段。通常使用分离字符的方法来处理数字数组。 维基百科【高精度计算】
在一、二辑中，给大家介绍了如何使用matlab自带工具箱以及大神John D'Errico开发的工具箱实现高精度计算。本辑作为用matlab做高精度计算的压轴辑，将给大家介绍一款效率远超前面两辑中所介绍的工具箱的高精度计算神器 —— Multiprecision Computing Toolbox for MATLAB (AdvanpixMCT)。
正式开始之前先通过官方的测试对比表来看看AdvanpixMCT工具相较于自带VPA、Maple以及Mathematica效率：
(截自www.advanpix.com)
(截自www.advanpix.com)
从上面两张表不难看出，无论是矩阵计算还是常规计算，AdvanpixMCT都力压另外三者，计算效率可以说是远远超越matlab的VPA以及Maple，大幅领先Mathematica，我想这也是为什么AdvanpixMCT会成付费工具箱的原因。
AdvanpixMCT提供的计算支持涵盖如下领域：
- 实数和复数、全矩阵和稀疏矩阵、多维数组
- 初等和特殊数学函数
- 线性方程组的求解器（包括直接和迭代稀疏求解器）
- 矩阵分析函数和因式分解
- 特征值和特征向量，包括广义和大规模问题。
- 奇异值分解
- 非线性方程组的求解器（使用Levenberg-Marquardt和其他信任区域方法进行fsolve）
- 数值积分（包括自适应quadgk和全套高斯正交）
- 优化和多项式
- 常微分方程求解器
- 数据分析和傅里叶变换
- 数论函数
前两辑中关于如何定义和使用高精度计算工具箱已经讲得非常多了，AdvanpixMCT的使用跟它们并无太多差异。在安装好免费7天试用版后，可以通过以下代码块中示例代码进行测试。
>> mp.Digits(34); % Setup default precision to 34 decimal digits (quadruple).
>> format longG % Toolbox shows all digits in case of 'long' formats.
% Simple expressions evaluation in quadruple precision:
>> x = mp('pi/4')
x =
0.7853981633974483096156608458198757
>> y = mp('sqrt(2)/2')
y =
0.707106781186547524400844362104849
>> A = repmat(x,10,10); % Create mp-matrix by scalar replication.
>> norm(y-cos(A)) % Calculations are done with 34 digits precision.
ans =
9.629649721936179265279889712924637e-35
% Assemble matrix row by row:
% (Butcher tableau for Gauss-Legendre method: https://goo.gl/izuFWg)
>> a1 = mp('[ 5/36 2/9-sqrt(15)/15 5/36-sqrt(15)/30 ]');
>> a2 = mp('[ 5/36+sqrt(15)/24 2/9 5/36-sqrt(15)/24 ]');
>> a3 = mp('[ 5/36+sqrt(15)/30 2/9+sqrt(15)/15 5/36 ]');
>> a = [a1;a2;a3]; % Concatenate rows into final matrix
% Create mp-matrices by conversion from double
>> A = mp(rand(5));
>> B = mp(eye(5));
>> [V,D] = eig(A,B);
>> norm(A*V - B*V*D,1)/(norm(A,1) * norm(B,1))
ans =
8.1433805952291741154581605524001229e-34
% Create sparse matrix with accumulation using triplets:
>> i = [6 6 6 5 10 10 9 9]';
>> j = [1 1 1 2 3 3 10 10]';
>> v = mp('[100 202 173 305 410 550 323 121]')'; % Prepare vector of nonzeros
>> S = sparse(i,j,v) % Form quadruple precision sparse matrix
S =
(6,1) 475
(5,2) 305
(10,3) 960
(9,10) 444高精度计算真的有那么重要吗？在某些情况下，还非得使用高精度计算才好使，比如处理病态特征值问题，目前唯一可靠的办法就是通过扩展计算精度来的达到较准确的计算。下面通过AdvanpixMCT提供的案例一起来看看精度对处理病态特征值问题到底又多重要，这里选用特征值敏感的Grcar矩阵来作为演示。Grcar矩阵是只含有-1，0，1三种元素的特征矩阵，在matlab中可以通过调用galleray函数实现Grcar矩阵的生成，如8*8的Grcar矩阵：
gallery('grcar',8)
ans =
1 1 1 1 0 0 0 0
-1 1 1 1 1 0 0 0
0 -1 1 1 1 1 0 0
0 0 -1 1 1 1 1 0
0 0 0 -1 1 1 1 1
0 0 0 0 -1 1 1 1
0 0 0 0 0 -1 1 1
0 0 0 0 0 0 -1 1理论表明，Grcar矩阵与其转置矩阵应该具有相同的特征值。下面就通过简单的代码来看看什么叫差之毫厘、谬以千里。
% MATLAB浮点数精度
A = -gallery('grcar',150);
figure('Color','w');
plot(eig(A),'k.'), hold on, axis equal
plot(eig(A'),'ro')
% MATLAB VPA高精度计算(34位精度)
digits(34);
A = vpa(-gallery('grcar',150));
eA = eig(A);
cA = eig(A');
figure('Color','w');
plot(eA,'k.'), hold on, axis equal
plot(cA,'ro')
（a) 浮点数计算结果
(b) VPA高精度计算结果
图(a)、(b)中，黑色点是150*150大小Grcar矩阵特征值图，红色是其转置矩阵特征值图。尽管Grcar矩阵的条件数cond(A) = cond(A') = 3.6106不高，但是使用双精度浮点数计算依然导致了极大的误差产生。
基本的介绍到此就接近尾声了，感兴趣的伙伴可以自己根据帮助文档研究更多应用场景。
AdvanpixMCT作为强大MATLAB高精度计算工具箱，其售价也不算便宜，学术单用户版售价249美元，而企业单用户版售价996美元。对于学生用户，可通过发送邮件获取专属的5折优惠。即便五折，换成人民币也是830多元左右。本来也是想入手的，但是想想还是有点舍不得，800元都可以买Endnote 20版的个人版授权了(PS：咱已入坑)。AdvanpixMCT采用的是VMProtect强加密方法以及一些特殊的文件关联方式，安装之后如超过7天试用期，即便卸载之后重新安装依然无法再次使用，除非重装系统。
不过咱也想到一种笨笨的办法，那就通过虚拟机快照的方式来实现。首先，在虚拟机中安装好MATLAB并关机，再者创建虚拟机快照，之后再开机，然后安装AdvanpixMCT，使用期结束后，使用快照将虚拟机系统还原，然后再次安装AdvanpixMCT，又可以试用7天，如此循环往复，直至千秋万代。温馨提示：AdvanpixMCT试用版目前是没有任何功能限制的，即全部功能都是能够试用。
参考资料：www.advanpix.com
腾讯云开发者
