the matrix addition should come 0 but it is not happening

1 visualización (últimos 30 días)
Pranshu
Pranshu el 28 de Jun. de 2024
Comentada: Aquatris el 28 de Jun. de 2024
i am adding a matrix in for loop
Dtb=zeros(3,3);
Dtrb=Dtb;
Drrb=Dtb;
Dts=zeros(2,2);
for N=1:nlayer
h1=-h/2+(N-1)*th;
h2=h1+th;
ah1=h2-h1;
bh1=(h2^2-h1^2)/2;
ch1=(h2^3-h1^3)/3;
% -----------
if(N <= nlayer/2 || N == 1)
theta=(1+(-1)^N)*pi/4;
else
theta = -1*(-1+(-1)^N)*pi/4;
end
% --------------
theta = rad2deg( theta );
t1=(cosd(theta))^4;
t2=(sind(theta))^4;
t3=(sind(theta)*cosd(theta))^2;
t4=sind(theta)*(cosd(theta))^3;
t5=cosd(theta)*(sind(theta))^3;
%**********************************************************
c111=c11*t1 + 2*(c12+2*c66)*t3 + c22*t2;
c121=(c11 + c22 - 4*c66)*t3 + c12*(t1+t2);
c221=c11*t2 + 2*(c12+2*c66)*t3 + c22*t1;
c161=(c11-c12-2*c66)*t4 + (c12-c22+2*c66)*t5;
c261=(c11-c12-2*c66)*t5 + (c12-c22+2*c66)*t4;
c661=(c11+c22-2*c12-2*c66)*t3 + c66*(t1+t2);
c441=c44*(cos(theta))^2 + c55*(sin(theta))^2;
c451=(c55-c44)*cos(theta)*sin(theta);
c551=c55*(cos(theta))^2 + c44*(sin(theta))^2;
Cb=[c111 c121 c161; c121 c221 c261; c161 c261 c661];
Cs=[c551 c451;c451 c441];
Dtrb = Dtrb + bh1*Cb;
end
but the resulting Dtrb atrix is not coming 0,
as you can see from fig: matrix for iteration 10 is bh1*Cb matrix which is added to Dtrb matrix and in figure both are same but opp in sign so if we add it should come 0 but as you can see in fig the result is not 0 but it is something * e-11 coming
  1 comentario
Avni Agrawal
Avni Agrawal el 28 de Jun. de 2024
Hi, can you share the constant values used for the above mentioned code?

Iniciar sesión para comentar.

Respuestas (2)

Aquatris
Aquatris el 28 de Jun. de 2024
Editada: Aquatris el 28 de Jun. de 2024
Given the difference is 1e-12, it is probably floating point issue. Floating points list significant bit is generally unreliable.
x = 1/3; % should be 0.3333333...
sprintf('%0.55f',x) % not really
ans = '0.3333333333333333148296162562473909929394721984863281250'
y = 1e10;
z = x + y - y; % lets add a large number to x and then substract the same number
sprintf('%0.55f',z) % should be the same as x but not really due to floating point precision
ans = '0.3333339691162109375000000000000000000000000000000000000'
  2 comentarios
Pranshu
Pranshu el 28 de Jun. de 2024
then how to get rid of this thing?
Aquatris
Aquatris el 28 de Jun. de 2024
You cannot get rid of it. It is a limitation of representing real numbers in a fixed space. Welcome to the world of computer arithmetic :D

Iniciar sesión para comentar.


SACHIN KHANDELWAL
SACHIN KHANDELWAL el 28 de Jun. de 2024
It would be difficult to debug the issue without sufficient code files. I suspect that the problem is due to precision limitations of floating point representation. The MATLAB values appear to be the same, but when you add them, discrepancies might occur.
Following code might be helpful to understand the above issue :
a1 = -3.45671;
a2 = 1.34576;
b1 = -4.23678;
b2 = 0;
matA = [a1,a2;b1,b2]
matA = 2x2
-3.4567 1.3458 -4.2368 0
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
a1 = 3.4567145;
a2 = -1.3457689;
b1 = 4.23679345;
b2 = 0;
matB = [a1,a2;b1,b2]
matB = 2x2
3.4567 -1.3458 4.2368 0
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
matA+matB
ans = 2x2
1.0e-04 * 0.0450 -0.0890 0.1345 0
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
The issue arises due to the precision limitations of floating-point representation in MATLAB. Even though "matA" and "matB" appear to be the same up to four decimal places, they are not exactly equal due to the additional precision in "matB".
Additionally, you may want to use 'format' function in MATLAB. Kindly, refer the following MathWorks documentation for more information :
Hope this is helpful!
Thanks.
  1 comentario
Pranshu
Pranshu el 28 de Jun. de 2024
your matrix values are different ,but my values are exactly same as you can see in fig, this B matrix is bh1*Cb

Iniciar sesión para comentar.

Categorías

Más información sobre Logical en Help Center y File Exchange.

Community Treasure Hunt

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

Start Hunting!

Translated by