This calculate error is why?

조회 수: 4 (최근 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 CenterFile Exchange에서 Logical에 대해 자세히 알아보기

태그

Community Treasure Hunt

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

Start Hunting!

Translated by