This calculate error is why?

Question

Xianjie 2012-11-10

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/53399-this-calculate-error-is-why

When I calculate some integration, normally I use the for loop the code is as following:

M=10;N=10;P=4;
    for m1=-P:M
      for n1=-P:N
          s=(m1+P)*(N+1+P)+n1+P+1;
          for m=-P:M
              for n=-P:N
                  t=(m+P)*(N+1+P)+n+P+1;
                  A2(s,t)=intefai1(a,m1,m);
              end
          end
      end
  end
for m1=-P:M
    for n1=-P:N
        s=(m1+P)*(N+1+P)+n1+P+1;
        for m=-P:M
            for n=-P:N
                t=(m+P)*(N+1+P)+n+P+1;
                B2(s,t)=intefai1(b,n1,n);
            end
        end
    end
end
function y = intefai1(a,m1,m)
% fai(m1,x)*fai(m,x)
if m1 >= 0
    if m >= 0
        y = alpha1(a,m1,m);
    else
        y = -gam(a,m1,abs(m));
    end
else
    if m >= 0
        y = -gam(a,m,abs(m1));
    else
        y = betass(a,abs(m),abs(m1));
    end
end
function y=gam(a,m,n)
% cos(m*pi*x/a)*sin(n*pi*x/a)
if m==0
    if n==m
        y=0;
    else
        y=-a*(-1+(-1)^n)/(n*pi);
    end
else
    if n==m
        y=0;
    else
        y=a*n*(-1+(-1)^(m+n))/(pi*(m^2-n^2));
    end
end
function y=betass(a,m,n)
% sin(m*pi*x/a)*sin(n*pi*x/a)
y=a/2*dell(m,n)*(1-dell(0,m))*(1-dell(0,n));
function y=alpha1(a,m,n)
% cos(m*pi*x/a)*cos(n*pi*x/a)
y=a*(dell(0,m)*dell(0,n)+1/2*dell(m,n)*(1-dell(0,m)*dell(0,n)));

When I use logical operation like the following

    function y=betass(a,m,n)
    % sin(m*pi*x/a)*sin(n*pi*x/a)
    y = a/2*(m==n).*(~(m==0)).*(~(n==0));
   function y=alpha1(a,m,n)
   y=a*((~m&~n)+1/2*(m==n).*(m|n));
function y=gam(a,m,n)
% cos(m*pi*x/a)*sin(n*pi*x/a)
 y = a*n.*((-1).^(m+n)-1)./(m.^2-n.^2)/pi;
 y(isnan(y)) = 0;
M=10;N=10;P=4;
[xx,yy]=ndgrid(-P:M,-P:N);
cc = -P:M;
dd = kron(cc',ones(N+P+1,1));
mm1 = repmat(dd,1,(N+P+1)*(M+P+1));
nn1 = repmat(xx,N+P+1,M+P+1);
ee = kron(cc,ones(1,N+P+1));
mm = repmat(ee,(M+P+1)*(N+P+1),1);
nn = repmat(yy,(N+P+1),(M+P+1));
A1=intefai1(a,mm1,mm);
B1=intefai1(b,nn1,nn);

Like this, the A1, B1, A2 and B2 are not the same fully. Some element is different. The difference is 6E-19. Why only some element is different, when mm=0, mm1=5,mm1=11, the element is different. Thank you for your sugestion.

2 个评论
显示无隐藏无

Matt Fig 2012-11-10

Please properly format your code so it is readable.

Xianjie 2012-11-11

Ok. Thank you. Now I format my code. Can you give me some sugestions? Thank you for your concern.

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Walter Roberson 2012-11-10

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/53399-this-calculate-error-is-why#answer_65010

http://matlab.wikia.com/wiki/FAQ#Why_is_0.3_-_0.2_-_0.1_.28or_similar.29_not_equal_to_zero.3F

2 个评论
显示无隐藏无

Xianjie 2012-11-11

Thank you for your advice. But only some function can encounter this phenomenon.

Walter Roberson 2012-11-11

Any time you compute the same quantity in two different ways, the problem can arise.

请先登录，再进行评论。

Answer 2

Dr. Seis 2012-11-12

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/53399-this-calculate-error-is-why#answer_65064

在 MATLAB Online 中打开

Instead of determining whether things are equal... determine whether the difference between two numbers are within a tolerance level. For example, instead of

if a==b
   % do something
end

try

if abs(a-b) < 1e-10 
   % do something
end

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Dr. Seis 2012-11-12

Incidentally, if my suggestion for trying to vectorize you code did not help in your previous question, feel free to ask follow-up questions. That is how we learn more about what you are looking for and how you learn more about the hints we suggest.

For example, I would have benefited from knowing the information you provided here in this question for answering your previous question.

请先登录，再进行评论。