注册 登录  
 加关注
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

信息 灵感 创新

I? =Information,Inspiration,Innovation

 
 
 

日志

 
 
关于我

we are 5. Mathematics, Computation, Programming, Engineering, and Making fun of life.

网易考拉推荐

the vpa function of matlab  

2012-06-02 14:00:12|  分类: M&M |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |
official help:
vpa:Variable precision arithmetic
Syntax:
R = vpa(A)
R = vpa(A) uses variable-precision arithmetic (VPA) to compute each element of A to at least d decimal digits of accuracy, where d is the current setting of digits. The toolbox also increases the internal precision of calculations by several digits (guard digits).

R = vpa(A, d)
R = vpa(A, d) uses at least d significant (nonzero) digits, instead of the current setting of digits. The value d must be a positive integer larger than 1 and smaller than 2^29+1
简单地说,这个函数就是表达式计算,将字符串理解为算式进行求解。
vpa(A)使用默认的有效数位计算,精度可以用命令digits显式,默认值为32,这个值必须介于1到2^29+1
例如:
vpa('1/2+3')
ans=3.5
但是vpa并不能对内存空间中既存在的变量做替换,例如:
>> A=[1 2 3;3 2 1;2 3 1];
>> B=[3 1 2;1 3 2;2 1 3];
>> vpa('A*B')
ans =
A*B
而不是我们期待的值:
>> A*B
ans =
    11    10    15
    13    10    13
    11    12    13
这个函数并不是用来做数值计算的,更多的时候是进行数值代换的,例如编写如下M函数:
function f=fract(n)
    syms a;
    if(n>0)
        f=a+1/fract(n-1);
    else
        f=a;
    end
end
这是一个连分数计算的函数,显然n趋向无穷大的时候,这个函数的值趋向于(sqrt(n^2+4)+n)/2
然后使用命令:
>> q=fract(7)
q =
a + 1/(a + 1/(a + 1/(a + 1/(a + 1/(a + 1/(a + 1/a))))))
>> subs(q,'a',2)
ans =
    2.4142
>> vpa(subs(q,'a',2),10)
ans =
2.414215686
>> vpa(subs(q,'a',2),20)
ans =
2.4142156862745098039
注意计算的准确结果是sqrt(2)+1。
其中subs是在符号表达式中进行表达式替换。

另一个有趣的例子是计算圆周率。大家都知道该值为圆周长和直径的比,通过圆内接正多边形或者圆外切正多边形进行计算,可以推导出:
the vpa function of matlab - Castor - 趁年轻,多折腾
M函数: 
function pi=mypi(n)
    if n>0
        pi=sqrt(2+mypi(n-1));
    else
        pi=0;
    end
end
计算方法:
>> n=20;
>> vpa(2^n*sqrt(2-mypi(n-1)))
ans =
3.1415965537048196857540970841798
计算结果与真实值比较接近的。

  评论这张
 
阅读(2068)| 评论(0)
推荐 转载

历史上的今天

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2016