Question

matlab问题,我编写的程序有人知道错哪了吗?详见问题补充

程序:
syms z;
r=250;
l=500;
for N=1:30
c=N/pi/r^2;
b=acos((l-z)./sqrt(l^2+r^2-2*l.*z));
a=r^2.*(acos(z/r)-z.*sqrt(r^2-z.^2)/r)+(l^2+r^2-2*l*z).*(b-sqrt((r^2-z.^2)./(l^2+r^2-2*l.*z)));
F=@(z)(2.*z.*l*c*b*exp(-c*a));
E(N)=quad(F,0,r);
end
plot(E,'-ok')
问题:
??? Undefined function or method 'isfinite' for input arguments of type 'sym'.
Error in ==> quad at 83
if ~isfinite(y(1))
有人能告诉我,我写的程序哪里错误吗?

Answer 1

原因

在前面的代码中，a是一个关于z的符号表达式，这样导致匿名函数F也是一个符号表达式，而quad*系列函数是用来计算数值积分的，不能处理sym类型的表达式。

 

解决方法

1、把循环体斗缓余中最后两句

    F=@(z)(2.*z.*l*c*b*exp(-c*a));
    E(N)=quad(F,0,r);

改成

    F=2.*z.*l*c*b*exp(-c*a);
    E(N)=quadgk(@(z)double(subs(F)),0,r);

这种方法使用符号运算，效率比较低（在我的电脑上计算一步哪李平均需要大约10秒）。另外空滚，注意积分函数使用quadgk，用quad可能算不出。

 

2、把最后两句改成（和上一种方法仅最后一句不同）

    F=2.*z.*l*c*b*exp(-c*a);
    E(N)=quadgk(eval(['@(z)' vectorize(F)]),0,r);

这种方法效率比上一种高很多，我在同一台电脑上测试，全部计算完只需要1秒多一点的时间。但在2008b上会遇到一个莫名其妙的问题——pi被表示成大写的PI，导致计算出错，已知在2007b上没问题，其它版本的情况不清楚。

 

3、计算出的结果是随N按指数规律增大的，大约在N=16就超出了双精度浮点数的表示范围而得到无穷大。如果把纵轴按照对数分度，画出的曲线大致如下：

你确定这是你想要的结果？