This calculate error is why?

8 views (last 30 days)
Xianjie
Xianjie on 10 Nov 2012
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 Comments
Matt Fig
Matt Fig on 10 Nov 2012
Please properly format your code so it is readable.
Xianjie
Xianjie on 11 Nov 2012
Ok. Thank you. Now I format my code. Can you give me some sugestions? Thank you for your concern.

Sign in to comment.

Answers (2)

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

Sign in to comment.


Dr. Seis
Dr. Seis on 12 Nov 2012
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 Comment
Dr. Seis
Dr. Seis on 12 Nov 2012
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.

Sign in to comment.

Tags

Community Treasure Hunt

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

Start Hunting!