如何用matlab做高精度计算?【第一辑】

a9873242c49701962a07d7436df1c0d5.png

高精度计算是一种程序设计的算法。由于中央处理器的字长限制,如32位CPU中一个整数最大只能取值4,294,967,295(=2^32-1),因此在超范围数值计算中,往往要采用模拟手段。通常使用分离字符的方法来处理数字数组。

维基百科【高精度计算】

对于跟咱一样的普通使用者而言,往往并不关心如何去实现高精度计算,更不会去研究相应的算法。咱这里讲的高精度计算也指的是计算过程中保持数据的精度不丢失。因为内容较多,计划分成三辑进行分享。

第一辑主要介绍matlab自带的高精度计算工具;第二辑主要介绍来自于File Exchange中的两款高精度计算工具箱;第三辑主要介绍一款收费的高精度计算工具箱Multiprecision Computing Toolbox。

matlab自带的高精度计算工具主要依赖于Symbolic Math Toolbox工具箱。因此,想要使用matlab自带的高精度计算工具,务必需要安装Symbolic Math Toolbox工具箱。

两个关键函数 —— vpadigits

digits函数使用来控制vpa的计算精度,使用前,按如下方式设置想要的精度即可:digits(num),num为设置的精度位数。示例如下:

digits(45)  % 
>> vpa(pi)
ans = 3.1415926535897932384626433832795028841971694
>> length('3.1415926535897932384626433832795028841971694')
ans = 45

在计算过程中,精度越高所花费的时间会相应增加,精度越低所花费的时间会相应减少,因此需要在精度与计算时间间做一个权衡。默认情况下,MATLAB使用16位精度。而使用vpa可获得更高的精度,vpa的默认精度为32位,即在未使用digits进行精度位数定义时,digits的返回值为32。如重启matlab后,在命令窗口用π来测试:

>> digits
 Digits = 32
>> dpi = vpa(pi)
dpi = 3.1415926535897932384626433832795
>> length(char(dpi))
ans = 33            % 此处为33是因为有小数点

正如上面所讲matlab的默认浮点数是16位精度,若超过16位则使用四舍五入法进行截断,如下图所示。

b1917932cb0d30ea547bb3f6cf66ec48.png

那怎么来实现让MATLAB存储高精度数值呢?就得用到咱们主角vpa函数了,使用vpa处理的数据自动为sym型数据。而sym本身也是创建符号变量、表达式、函数、矩阵等函数,且其所创变量、表达式、函数、矩阵同样为sym型数据。

一、创建高精度数据

小伙伴们可能会想直接使用vpa或sym是不是就可以搞定了,那咱们就来看看下面的示例,定义一个高精度小数:3.141514546465512487984541。

>> a = sym(3.141514546465512487984541)
a = 7074061870420537/2251799813685248
>> vpa(a)
 ans = 3.1415145464655123141994863544824
>> b = vpa(3.141514546465512487984541)
b = 3.1415145464655123141994863544824
>> vpa(b)
ans = 3.1415145464655123141994863544824

从结果不难看出,无论直接使用sym还是vpa都未能正确创建咱们想要的数字,这是为什么呢?究其原因,就是MATLAB默认浮点数精度在作怪,超过16位精度的数字就会自动被截断,而其后被无效数值填充,导致无法得出正确结果。那要怎么解决此问题呢?问题根源既然出在数值上,那咱不用数值不就可以了吗?要不试试看?咱直接将数值定义成字符串,然后再使用sym或vpa。

>> a = sym('3.141514546465512487984541')
a = 3.141514546465512487984541
>> b = vpa('3.141514546465512487984541')
b = 3.141514546465512487984541

显然咱的方案是正确的。那有的小伙伴就问了,用num2str来转换上面的数值可以吗?不用怀疑,肯定是不会对的,不信的小伙伴可以私下试试。

二、sym型数据数值化

使用vpa函数可以将sym型数据数值化,通过char函数可将数值化的sym型数据转成可用的字符串型数据。如下:

>> p = sym(pi);
>> a = sym('3.141514546465512487984541');
>> b = sym(1/2);
>> f = a*sin(b*p);
>> val = vpa(f)
val = 3.141514546465512487984541
% vpa有效的保留计算精度,f和a完全相等

三、更改数据的精度

除了使用digits定义精度位数外,使用vpa也可以直接定义计算精度,如取pi的前128位:

>> pi128 = vpa(pi,128)
pi128 = 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446

四、符号计算结果数值化

虽然符号计算的结果是精确的,但可能是一种不便于阅读的形式,如使用solve求解高次多项式的根,它给出的不是具体解的值,而是用root来表示根。

>> syms x
y = solve(x^3 - x + 1, x)
y = root(z^3 - z + 1, z, 1)
    root(z^3 - z + 1, z, 2)
    root(z^3 - z + 1, z, 3)
>> vvpa = vpa(y)
vvpa =
0.66235897862237301298045442723905 - 0.56227951206230124389918214490937i
0.66235897862237301298045442723905 + 0.56227951206230124389918214490937i
                                      -1.3247179572447460259609088544781

由于篇幅不宜过长,更多有关vpa的应用小伙伴可以自行研究使用。接下来会在另外两辑中介绍第三方高精度计算工具。欲知后事如何,且看下回分解!

参考资料:

[1] www.mathworks.com/help/symbolic/vpa.html

[2] https://zh.m.wikipedia.org/zh-hans/高精度计算

dc10770bc211a801befbaa147de33565.png

如需转载,请在公众号中回复“转载”获取授权!

文章出处登录后可见!

已经登录?立即刷新

共计人评分,平均

到目前为止还没有投票!成为第一位评论此文章。

(0)
扎眼的阳光的头像扎眼的阳光普通用户
上一篇 2023年4月5日
下一篇 2023年4月5日

相关推荐