How to reduce the time of solving the lengthy piecewise terms

Hi experts,

I am trying to plot the derived analytical solution, however, when the triple summation terms is >10, the compute time seems could not finished up summing.

I have attached the .m files.

Any advise to further simplified the solving method ?

The triple summation terms start from C5 to C13.

Thanks,

Andy

Answer 1

[Walter Roberson]

On my system, the following version takes approximately 300 seconds in total.

tic;
syms n M N y real
assume([n M N y] >= 0);
Y = 0:0.02:1;
b =10; t=1; br=1; Pr=1;
%t = t*/pr
 
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 = subs( sum( subs( 2.*br.*C1.*sin(n.*pi.*y), n, 1:20 )), y, Y );
Tn = double(T);
display(Tn);
 
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).^4)).* (((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 , 0, A );
 
T1 = subs( sum( subs( sum( subs( 2.*br.*Az, n, 1:20 )), M, 1:20)), y, Y);
T1n = double(T1);
display(T1n);
 
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).^4).*((b.^2)+(N.*pi).^4))).*( ((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 | Pr.*(M.^2+N.^2)== n.^2 | Pr.^2.*((M.*pi).^4 + b.^2)== (2.*Pr.*(M.*n.*pi.^2).^2)-(n.*pi).^4 | Pr.^2.*((N.*pi).^4+ b.^2 )== (2.*Pr.*(N.*n.*pi.^2).^2)-(n.*pi).^4 , 0, B );

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

[ng, Mg, Ng] = ndgrid(1:20);
T2a = subs( 2.*br.*Bz, {n,M,N}, {ng, Mg, Ng});
T2b = sum(T2a(:));
T2 = subs(T2b, y, Y);
toc

tic
T2n = double(T2);
toc

display(T2n);
 
T3n=Tn+T1n+T2n;
 
SumT = table(Y(:),Tn(:), T1n(:), T2n(:),T3n(:), 'VariableNames', {'y','T', 'T1', 'T2','SumT'});
plot(T3n,Y);

[Chee Hao Hor]

Dear Walter,

Good day, sorry for late reply, i was trying to compare and understand the scripts.

As i observed from the updated code, it seems the trouble cames from the %T2 =subs(....), which make the compute time goes insanely long. As I have learned from you, the previous coding leaves y as symbolic until the end, then piecewise and lastly sub in the vector of specific values as the last step.

In this updated script, you are creating a 3D grid "ng,Mg,Ng", eg. in here, each length is 20, each grid division is 1.

Then we sub the vector of the specific values after the piecewise elimination, (while y leaves as symbolic).

After that, we sum up all the vectors, and lastly sub the symbolic y.

1. Did I get the understanding correctly ?
2. The previous script have difficulty to solve for lengthy sum, while this script discrete the lengthy sum so that the compute time is reduce ?

Thanks,

Andy

[Walter Roberson]

You are correct. Your code was substituting in 1:20 on one dimension, taking the sum(), substituting in 1:20 on a different dimension, taking the sum, substituing in 1:20 on a third dimension, and taking the sum. It turns out that taking the sum even once is really time consuming. So I substitute in a 3D grid of n, M, N values, which gives enough information to resolve each of the piecewise values, resulting in 0 being put in place for some of the locations, and resulting in something that is not piecewise() being put in place for the remaining locations. Then when you sum() those, you are not adding anything piecewise, which saves a lot of time, because when you add piecewise, MATLAB wants to create one larger piecewise() instead of leaving it as the sum of two piecewise()

[Chee Hao Hor]

Dear Walter,

Thanks for the explaination, now I understand the concept.

Have a good day ! :)

Andy.