How to skip the division of zero in three summation term.

F=n-M-N;
syms n M N
T= symsum(symsum(symsum((1/F),n,1,3),M,1,3),N,1,3);

Hi everybody,

I am trying to solve a problem, which involving three summation terms in the Temperature, denoted as T. Each of the sum range from 1 to 3.

However, there are certain value gives division by zero condition. How can I skip this condition while it add up others non-zero terms and give me a value ?

Desperately need the helps. Thanks.

Answer 1

[Walter Roberson]

syms n M N
F = n-M-N;
Fz = piecewise(F == 0, 0, 1/F);
T = symsum( symsum( symsum( Fz, n, 1, 3), M, 1, 3), N, 1, 3)

[Walter Roberson]

Better:

[n, M, N] = ndgrid(1:3, 1:3, 1:3);
F = n-M-N;
F(F == 0) = [];
T = sum(1./F);

[Chee Hao Hor]

Hi Mr Walter,

1. Do you mean that instead of using summation, we use matrix ?

[n, M, N] = ndgrid(1:3, 1:3, 1:3);

1. What does the third line mean ?

F(F == 0) = [];

Also, I dont really understand this line meaning

Fz = piecewise(F == 0, 0, 1/F);

Sorry, I am still new in Matlab. Really thanks, Looking forward for the explaine.

[Walter Roberson]

If you had vectors a b c, then

[A, B, C] = ndgrid(a, b, c);

would have the property that

A(J,K,L) = a(J)
B(J,K,L) = b(K)
C(J,K,L) = c(L)

Together these express all possible combinations of the inputs, in matrix form. A expresses values from the first vector, B expresses values from the second vector, C expresses values from the third vector. If you were then to do D = A + B.*C then somewhere in the output you would express every possible combinations of inputs: D(J,K,L) = a(J) + b(K) .* c(L) and this is true even if a, b, c are different lengths. D = A + B.*C would calculate all of the possible answers, simultaneously.

With n, M, n being 3D matrices, F = n-M-N is a 3D matrix result that inside it has all of the possible answers for selecting one entry from n, one entry from M, and one entry from N.

Now, some of those answers are 0, such as F(3,2,1) which would be 3-2-1 . You do not want to include those locations where F == 0 in the sum. So you delete them out of F.

F == 0 tests all of the entries in F and returns a logical array the same sizes as F, that is true where F is 0, and is false where F is not 0. Then F(F==0) is "logical indexing": it selects the locations of F precisely where F==0 is true -- so we have gone from the test F==0 to selecting the locations in F where the test is true. And what we do with those locations in F(F==0)=[] is delete them. In MATLAB, assigning an explicit [] to something means to delete it. (This does not hold true for a calculated empty array, only for literal [] or literal '' (two apostrophes in a row -- the empty character vector))

Thus, after F(F==0)=[]; then F has no more 0 elements, with only the combinations that led to 0 having been removed and all of the other results left in place. We can then proceed to take 1 divided by what is left and sum those to get the solution.

Your original framing of the problem with symsum() is a symbolic approach that only works with the Symbolic Toolbox. But since we know we are using the symbolic toolbox, we can take advantage of other tools in the symbolic toolbox: instead of summing 1/F, we can sum the results of a different expression. The different expression we will choose to sum is:

if F is 0, then use 0 as the expression value
otherwise if F was not 0, use 1/F as the expression value

Here we take advantage of the fact that when you are just taking a sum, then instead of removing a problem term from the sum, we can transform the problem term to 0, because 0 does not affect the sum. (The situation would be different if you were calculating mean(), as then the extra 0 would affect the count not just the total.)

[Chee Hao Hor]

Hi Mr Walter,

I am now understand the concept of the code line.

Actually, I am trying to find and plot an extremely cumbersome analytical solution I derived, this is the part of the condition I face.

I will try to modify and see whether it work or not. Would like to ask, are you interests to have a look and check after I had completed code ? Perhaps might take a few days to complete. If yes, could you leave your email to me ?

Anyway, I want to thank you sincerely, and wish you have a very good day.

[Walter Roberson]

If you post it here when it is done, then I might take a look at it.

[Chee Hao Hor]

Hi Walter,

I tried the piecewise, it works, but not for ndgird.

I stuck somewhere, suppose, I want to see both value of b = 1,2 how to code for this b variable ?

syms n M N
b=[1,2];
A = (b.^2/(((b.^2)+(M.*pi).^2).*((b.^2)+(N.*pi).^2)))*((1-cos((n+M+N).*pi))./((n+M+N).*pi)+(1-cos((n-M-N).*pi))./((n-M-N).*pi)+(1-cos((n+M-N).*pi))./((n+M-N).*pi)+(1-cos((n-M+N).*pi))./((n-M+N).*pi));
Az = piecewise( ((n-M-N).*pi) == 0 , 0,((n+M-N).*pi) == 0 , 0,((n-M+N).*pi) == 0 , 0, A);
T2 = symsum( symsum( symsum( Az,M , 1, 2), N, 1, 2), n, 1, 2);
double (T2)

It only works when b is equal to a single value.

[Walter Roberson]

Do you really mean to use / after b.^2 at the begining of A? Or do you mean to use ./ instead ?

[Chee Hao Hor]

Since it is elementary, perhaps it should be ./. My mistake.

[Chee Hao Hor]

Hi Walter,

I am actually trying to plot this solution as in the attachment. I need to plot a graph, lets say, y vs T, and with varying b, should I just use normal operator or elementary wise ? Please advise.

[Chee Hao Hor]

Good Morning Mr Walter,

I had tried to run this code, it works fine for the T and T1, which gaves me answer in less than 3 minutes time. However, when it solve for T2, the software keep running for 10 hours but still no solution. Could you give me a hand ? Looking forward your reply. Thanks.

syms n
y = [0:0.01:1];
b =1; t=1; br=1; Pr=1;
C1=((1-(-1)^n)/(2*n*pi))*((1-exp(-(n*pi)^2 *t ))/(n*pi)^2 +((n*pi)^2*exp(-(n*pi)^2*t)-((n*pi)^2*cos(2*b*Pr*t))-(2*b*Pr*sin(2*b*Pr*t)))/((4*b^2*Pr^2)+(n*pi)^4 ));
T= symsum ((2*br*C1*sin(n.*pi.*y)),n,1,2);
double(T)


syms n M 
y = [0:0.01:1];
b =1; t=1; Pr= 1; br=1;
C2=(((M.*pi).^2./((8.*b.^2.*Pr.^2)+2.*(n.*pi).^4)).*( (2.*b.*Pr.*exp(-(n.*pi).^2.*t))+((n.*pi).^2.*sin(2.*b.*Pr.*t))- (2.*b.*Pr.*cos(2.*b.*Pr.*t))));

C3=-((b./((8.*(n.*pi).^2.*b.^2.*Pr.^2)+ (2.*(n.*pi)^6))).*( (4.*b.^2.*Pr.^2.*exp(-(n.*pi).^2.*t))+ (2.*b.*Pr.*(n.*pi).^2.*sin(2.*b.*Pr.*t))+((n.*pi).^4.*cos(2.*b.*Pr.*t))-(4.*b.^2.*Pr.^2)-(n.*pi).^4 ));

C4=(((M.*pi).^2./(Pr.^2.*((M.*pi).^4+ b.^2 ))).* ( (b.*Pr.*exp(-(n.*pi).^2.*t))- exp(-(M.*pi).^2.*Pr.*t).*((((M.*pi).^2.*Pr-(n.*pi).^2).*sin(b.*Pr.*t))+ (b.*Pr.*cos(b.*Pr.*t)) )));

A = (2.*b./(b.^2+(M.*pi).^2)).* (((1-cos((n+M).*pi))./((n+M).*pi)) + ((1-cos((n-M).*pi))./((n-M).*pi))).*(C2+C3+C4).*sin(n.*pi.*y);
Az = piecewise( ((n-M).*pi) == 0 , 0, A);
T1 = symsum( symsum( 2.*br.*Az,M , 1, 2), n, 1, 2);
double(T1)


syms n M N
y = [0:0.01:1];
b= 1; t=1;Pr=1; br=1;
C5=(((M.*N.*pi.^2 ).^2./((8.*b^2.*Pr.^2.*(n.*pi).^2 )+ (2.*(n.*pi).^6))).*((2.*b.*Pr.*(n.*pi).^2.*sin(2.*b.*Pr.*t))+((n.*pi).^4.*cos(2.*b.*Pr.*t))+(4.*b.^2.*Pr.^2)+(n.*pi).^4-(exp(-(n.*pi).^2.*t).*(4.*b.^2.*Pr.^2+2.*(n.*pi).^4 ))));

C6=(((b.*(M.*pi).^2)./((8.*b.^2.*Pr.^2)+ (2.*(n.*pi).^4))).*(((n.*pi).^2.*sin(2.*b.*Pr.*t))- (2.*b.*Pr.*cos(2.*b.*Pr.*t))+ (2.*b.*Pr.*exp(-(n.*pi).^2.*t))));


B1 = (b.^2./(((b.^2)+(M.*pi).^2).*((b.^2)+(N.*pi).^2))).*((1-cos((n+M+N).*pi))./((n+M+N).*pi)+(1-cos((n-M-N).*pi))./((n-M-N).*pi)+(1-cos((n+M-N).*pi))./((n+M-N).*pi)+(1-cos((n-M+N).*pi))./((n-M+N).*pi));
B1z = piecewise( ((n-M-N).*pi) == 0 , 0,((n+M-N).*pi) == 0 , 0,((n-M+N).*pi) == 0 , 0, B1);
T2 = symsum( symsum( symsum( 2.*br.*B1z.*(C5+C6).*sin(n.*pi.*y),N , 1, 2), M, 1, 2), n, 1, 2);
double (T2)

[Walter Roberson]

Why do you do

piecewise( ((n-M).*pi) == 0, 0, A)

when you can just do

piecewise( n == M, 0, A)

and likewise

piecewise( ((n-M-N).*pi) == 0, etc)

is clearer

piecewise( n-M-N == 0, etc)

[Walter Roberson]

Code attached that completes all of the calcuations in less than 1 second.

[Chee Hao Hor]

Mr Walter,

YEAH ! It works for the rest of the parts !! very happy.

Suppose I want to export the result into excel file, Tn, T1n and T2n take 1 coloumn , I wrote this, but it say Dimensions of matrices being concatenated are not consistent. Please advise. :)

SumT = table(['Tn';'T1n';'T2n'],[Tn,T1n,T2n]);
writetable(SumT,'testdata.xls','Sheet',1)

[Walter Roberson]

SumT = table({'Tn';'T1n';'T2n'}, {Tn;T1n;T2n});

Note that this would be for a table with three rows. Although the table would only have two columns, writetable() would expand this out to 102 columns, most with names similar to Var2_38

Writing the data out to rows like this is probably not nearly as useful as if you had written to columns, which you would do with

SumT = table(Tn, T1n, T2n);

This would lead to three columns, each with meaningful names.

[Chee Hao Hor]

The first one works, all i need to do is to manually transport it.

However, when i tried the second code, it gaves me empty excel data. I dont know why.

[Walter Roberson]

SumT = table(Tn(:), T1n(:), T2n(:), 'VariableNames', {'T', 'T1', 'T2'});

[Chee Hao Hor]

Thanks Walter, now i am somewhere on it.

Actually, could I know what makes the long code too so long time to run as I wrote previously ? Which line is actually causing the problem ? All I could observed is the symsum( ...) disappears on yours code and it run so fast. I wish to learn from you.

[Walter Roberson]

The long stage is the way you use y. Notice in my versions I leave y symbolic until the end, so I am only dealing with two equations at a time until I sum() them to get down to one. So it all comes down to a single equation in y that has already had the piecewise resolved. Then I subs in the vector of specific values as the last step.

[Chee Hao Hor]

Orhs, I get it now.

Also, I try to put value of 5 5 5, instead of 2 2 2 for n M N, the results is NaN, what does this mean ? It is okay when it is below the value of 4, for n M N, respectively.

T2 = subs( sum( subs( sum( subs( sum( subs( 2.*br.*Bz, n, 1:5) ), M, 1:5 ) ), N, 1:5 ) ), y, Y );

[Walter Roberson]

My suspicion is that exp() is overflowing. I will check sometime tomorrow as I do not have access to matlab at the moment.

[Chee Hao Hor]

Alright Walter, I looking forward from you tomorrow, i put the code here, easy for you to access tomorrow. I skipped the part which has no problem.

syms n M N y real
assume([n M N y] >= 0);
Y = 0:0.01:1;
b =1; t=2; br=5; Pr=1;

C5=(((M.*N.*pi.^2 ).^2./((8.*b^2.*Pr.^2.*(n.*pi).^2 )+ (2.*(n.*pi).^6))).*((2.*b.*Pr.*(n.*pi).^2.*sin(2.*b.*Pr.*t))+((n.*pi).^4.*cos(2.*b.*Pr.*t))+(4.*b.^2.*Pr.^2)+(n.*pi).^4-(exp(-(n.*pi).^2.*t).*(4.*b.^2.*Pr.^2+2.*(n.*pi).^4 ))));
C6=(((b.*(M.*pi).^2)./((8.*b.^2.*Pr.^2)+ (2.*(n.*pi).^4))).*(((n.*pi).^2.*sin(2.*b.*Pr.*t))- (2.*b.*Pr.*cos(2.*b.*Pr.*t))+ (2.*b.*Pr.*exp(-(n.*pi).^2.*t))));
C7= -(((M.*N.*pi.^2 ).^2./((Pr.^2.*((N.*pi).^4 + b.^2 ))- (2.*Pr.*(N.*n.*pi.^2).^2) +(n.*pi).^4)).* ((exp(-(n.*pi).^2.*t).*(Pr.*(N.*pi).^2-(n.*pi).^2))+ (exp(-(N.*pi).^2.*Pr.*t).*(b.*Pr.*sin(b.*Pr.*t)+(((n.*pi).^2- (Pr.*(N.*pi).^2)).*cos(b.*Pr.*t))))));
C8=(((b.*(N./pi).^2)./((8.*b.^2.*Pr.^2)+(2.*(n.*pi).^4))).*(((n.*pi).^2.*sin(2.*b.*Pr.*t))-(2.*b.*Pr.*cos(2.*b.*Pr.*t))+ (2.*b.*Pr.*exp(-(n.*pi).^2.*t)) ));
C9=(((b.^2)./((8.*b.^2.*Pr.^2.*(n.*pi).^2)+(2.*(n.*pi).^6))).*((4.*b.^2.*Pr.^2)+((n.*pi).^4)- (2.*b.*Pr.*(n.*pi).^2.*sin(2.*b.*Pr.*t))-((n.*pi).^4.*cos(2.*b.*Pr.*t))- (4.*b.^2.*Pr.^2 .*exp(-(n.*pi).^2.*t)) ));
C10=(((b.*(N.*pi).^2)./((Pr.^2.*((N.*pi).^4+ b.^2 ))- (2.*Pr.*(N.*n.*pi.^2).^2)+(n.*pi).^4 )).*((exp(-(N.*pi).^2.*Pr.*t).*((((Pr.*(N.*pi).^2)-(n.*pi).^2).*sin(b.*Pr.*t))+ (b.*Pr.*cos(b.*Pr.*t))))- (b.*Pr.*exp(-(n.*pi).^2.*t)) ));
C11=-(((M.*N.*pi.^2 ).^2./((Pr.^2.*((M.*pi).^4+ b.^2 ))- (2.*Pr.*(M.*n.*pi.^2).^2)+(n.*pi).^4 )).*((exp(-(n.*pi).^2.*t).*((Pr.*(M.*pi).^2)-(n.*pi).^2))+ (exp(-(M.*pi).^2.*Pr.*t).*((b.*Pr.*sin(b.*Pr.*t))+(((n.*pi).^2- (Pr.*(M.*pi).^2)).*cos(b.*Pr.*t))))) );
C12=(((b.*(M.*pi).^2)./((Pr.^2.*((M.*pi).^4 + b.^2))- (2.*Pr.*(M.*n.*pi.^2).^2)+(n.*pi).^4)).*(exp(-(M.*pi).^2.*Pr.*t).*((((Pr.*(M.*pi).^2)-(n.*pi).^2).*sin(b.*Pr.*t))+ (b.*Pr.*cos(b.*Pr.*t)))- (b.*Pr.*exp(-(n.*pi).^2.*t)) ));
C13=(((M.*N.*pi.^2).^2./(pi.^2.*((Pr.*(M.^2+N.^2))-n.^2 ))) .*(exp(-(n.*pi).^2.*t)- (exp(-(M.*pi).^2.*Pr.*t).*exp(-(N.*pi).^2.*Pr.*t)) ));

B = (b.^2./(((b.^2)+(M.*pi).^2).*((b.^2)+(N.*pi).^2))).*( ((1-cos((n+M+N).*pi))./((n+M+N).*pi)) + ((1-cos((n-M-N).*pi))./((n-M-N).*pi)) + ((1-cos((n+M-N).*pi))./((n+M-N).*pi)) + ((1-cos((n-M+N).*pi))./((n-M+N).*pi))).*(C5+C6+C7+C8+C9+C10+C11+C12+C13).*sin(n.*pi.*y);
Bz = piecewise( n-M-N == 0 | n+M-N == 0 | n-M+N == 0, 0, B );

T2 = subs( sum( subs( sum( subs( sum( subs( 2.*br.*Bz, n, 1:5) ), M, 1:5 ) ), N, 1:5 ) ), y, Y );
T2n = double(T2);
disp(T2n);

[Walter Roberson]

In C13 you have a division by (pi.^2.*((Pr.*(M.^2+N.^2))-n.^2)) . Your pr is 1, so that becomes (pi.^2 * (M^2 + N^2 - n^2)) . Which gives a divison by 0 if the values satisfy a pythagorian triple, such as n = 5, N = 3, M = 4

[Chee Hao Hor]

I keep on looking at exp term yesterday. I had put this under the piecewise line, it works now !

I guess i have 1 last question. Suppose, solving this code, i can get 1 line in a graph of Y against T, at 1 value of b. To put 3 value of b in a single graph, i copy and paste the code 3 times. I try to change the b, from

b=10; 
%to
b=[1,2,3];

It still give me a value only. Please advise.