Asked by 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.

Answer by Walter Roberson
on 10 Nov 2012

Xianjie
on 11 Nov 2012

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

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.

Answer by 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

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.

Opportunities for recent engineering grads.

Apply Today
## 2 Comments

## Matt Fig (view profile)

## Direct link to this comment

https://in.mathworks.com/matlabcentral/answers/53399-this-calculate-error-is-why#comment_110567

## Xianjie (view profile)

## Direct link to this comment

https://in.mathworks.com/matlabcentral/answers/53399-this-calculate-error-is-why#comment_110618

Sign in to comment.