This calculate error is why?
3 views (last 30 days)
Show older comments
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
Answers (2)
Walter Roberson
on 10 Nov 2012
2 Comments
Walter Roberson
on 11 Nov 2012
Any time you compute the same quantity in two different ways, the problem can arise.
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
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.
See Also
Categories
Find more on Startup and Shutdown in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!