이 질문을 팔로우합니다.
- 팔로우하는 게시물 피드에서 업데이트를 확인할 수 있습니다.
- 정보 수신 기본 설정에 따라 이메일을 받을 수 있습니다.
How to skip the division of zero in three summation term.
조회 수: 6 (최근 30일)
이전 댓글 표시
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.
채택된 답변
Walter Roberson
2019년 4월 22일
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)
댓글 수: 23
Walter Roberson
2019년 4월 22일
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
2019년 4월 22일
편집: Chee Hao Hor
2019년 4월 22일
Hi Mr Walter,
- Do you mean that instead of using summation, we use matrix ?
[n, M, N] = ndgrid(1:3, 1:3, 1:3);
- 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
2019년 4월 22일
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
2019년 4월 22일
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.
Chee Hao Hor
2019년 4월 23일
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
2019년 4월 23일
Do you really mean to use / after b.^2 at the begining of A? Or do you mean to use ./ instead ?
Chee Hao Hor
2019년 4월 23일
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
2019년 4월 25일
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
2019년 4월 25일
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
2019년 4월 25일
Code attached that completes all of the calcuations in less than 1 second.
Chee Hao Hor
2019년 4월 25일
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
2019년 4월 25일
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
2019년 4월 25일
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
2019년 4월 25일
SumT = table(Tn(:), T1n(:), T2n(:), 'VariableNames', {'T', 'T1', 'T2'});
Chee Hao Hor
2019년 4월 25일
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
2019년 4월 25일
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
2019년 4월 25일
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
2019년 4월 25일
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
2019년 4월 25일
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
2019년 4월 25일
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
2019년 4월 26일
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.
추가 답변 (0개)
참고 항목
카테고리
Help Center 및 File Exchange에서 Calculus에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!오류 발생
페이지가 변경되었기 때문에 동작을 완료할 수 없습니다. 업데이트된 상태를 보려면 페이지를 다시 불러오십시오.
웹사이트 선택
번역된 콘텐츠를 보고 지역별 이벤트와 혜택을 살펴보려면 웹사이트를 선택하십시오. 현재 계신 지역에 따라 다음 웹사이트를 권장합니다:
또한 다음 목록에서 웹사이트를 선택하실 수도 있습니다.
사이트 성능 최적화 방법
최고의 사이트 성능을 위해 중국 사이트(중국어 또는 영어)를 선택하십시오. 현재 계신 지역에서는 다른 국가의 MathWorks 사이트 방문이 최적화되지 않았습니다.
미주
- América Latina (Español)
- Canada (English)
- United States (English)
유럽
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
아시아 태평양
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)