This calculate error is why?

1 次查看（过去 30 天）
Xianjie
Xianjie 2012-11-10
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
Matt Fig 2012-11-10
Please properly format your code so it is readable.
Xianjie
Xianjie 2012-11-11
Ok. Thank you. Now I format my code. Can you give me some sugestions? Thank you for your concern.

请先登录，再进行评论。

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

回答（2 个）

Walter Roberson
Walter Roberson 2012-11-10
  2 个评论
Xianjie
Xianjie 2012-11-11
Thank you for your advice. But only some function can encounter this phenomenon.
Walter Roberson
Walter Roberson 2012-11-11
Any time you compute the same quantity in two different ways, the problem can arise.

请先登录，再进行评论。

Dr. Seis
Dr. Seis 2012-11-12
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 个评论
Dr. Seis
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.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Startup and Shutdown 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by