For Loop for Symbolic
2 views (last 30 days)
Show older comments
Syed Abdul Rafay
on 25 Oct 2023
Commented: Star Strider
on 26 Oct 2023
syms lambda
syms x
syms i
R=4;
rho_L= 2702;
rho_r= 5700;
E_L= 70e+9;
E_r= 200e+9;
nu=3;
A=zeros(R+1,R+1);
for r=0:R
for m=0:R
M1= (rho_L/(exp(nu)-1))*(symsum(((-1)^(r-i))*((factorial(r))/(factorial(i)*nu^(r-i+1))),i,0,r));
M2= (rho_r/(exp(nu)-1))*(symsum(((-1)^(r-i))*((factorial(r))/(factorial(i)*nu^(r-i+1))),i,0,r));
M3= (rho_L/(exp(nu)-1))*(symsum(((-1)^(r+1-i))*((factorial(r+1))/(factorial(i)*nu^(r-i+2))),i,0,r+1));
M4= (rho_r/(exp(nu)-1))*(symsum(((-1)^(r+1-i))*((factorial(r+1))/(factorial(i)*nu^(r-i+2))),i,0,r+1));
M5= (rho_L/(exp(nu)-1))*(symsum(((-1)^(r+2-i))*((factorial(r+2))/(factorial(i)*nu^(r-i+3))),i,0,r+2));
M6= (rho_r/(exp(nu)-1))*(symsum(((-1)^(r+2-i))*((factorial(r+2))/(factorial(i)*nu^(r-i+3))),i,0,r+2));
M7= (rho_L/(exp(nu)-1))*(symsum(((-1)^(r+3-i))*((factorial(r+3))/(factorial(i)*nu^(r-i+4))),i,0,r+3));
M8= (rho_r/(exp(nu)-1))*(symsum(((-1)^(r+3-i))*((factorial(r+3))/(factorial(i)*nu^(r-i+4))),i,0,r+3));
M9= ((E_L*nu)/(exp(nu)-1))*(symsum(((-1)^(r-i))*((factorial(r))/(factorial(i)*nu^(r-i+1))),i,0,r));
M10= ((E_r*nu)/(exp(nu)-1))*(symsum(((-1)^(r-i))*((factorial(r))/(factorial(i)*nu^(r-i+1))),i,0,r));
M11= ((E_L*nu^2)/(exp(nu)-1))*(symsum(((-1)^(r-i))*((factorial(r))/(factorial(i)*nu^(r-i+1))),i,0,r));
M12= ((E_r*nu^2)/(exp(nu)-1))*(symsum(((-1)^(r-i))*((factorial(r))/(factorial(i)*nu^(r-i+1))),i,0,r));
M13= ((E_L*exp(nu))/(exp(nu)-1))*(symsum(((-1)^(r+m-i))*((factorial(r+m))/(factorial(i)*nu^(r+m-i+1))),i,0,r));
M14= ((E_r*exp(nu))/(exp(nu)-1))*(symsum(((-1)^(r+m-i))*((factorial(r+m))/(factorial(i)*nu^(r+m-i+1))),i,0,r));
J1= ((rho_L/(r+1))*(1+(1/(exp(nu)-1))))-((M1-M2)*exp(nu))-((rho_r/(exp(nu)-1))*(1/(r+1)));
J2= ((rho_L/(r+2))*(1+(1/(exp(nu)-1))))-((M3-M4)*exp(nu))-((rho_r/(exp(nu)-1))*(1/(r+2)));
I1=(-M9*exp(x*nu)*(x^r))+(M10*exp(x*nu)*(x^r));
I2=(M11*exp(x*nu)*(x^r))+(M12*exp(x*nu)*(x^r));
I3=((rho_L*(x^(r+1)))/(r+1))-(M1*exp(x*nu)*(x^r))+((rho_L*(x^(r+1)))/((r+1)*(exp(nu)-1)))-((rho_r/(exp(nu)-1))*((x^(r+1))/(r+1)))+((M2*exp(x*nu)*(x^r)));
I4=((rho_L*(x^(r+2)))/(r+2))-(M1*exp(x*nu)*(x^(r+1)))+((rho_L*(x^(r+2)))/((r+2)*(exp(nu)-1)))-((rho_r/(exp(nu)-1))*((x^(r+2))/(r+2)))+((M2*exp(x*nu)*(x^(r+1))));
I5=((rho_L*(x^(r+3)))/(r+3))-(M1*exp(x*nu)*(x^(r+2)))+((rho_L*(x^(r+3)))/((r+3)*(exp(nu)-1)))-((rho_r/(exp(nu)-1))*((x^(r+3))/(r+3)))+((M2*exp(x*nu)*(x^(r+2))));
I6=((rho_L*(x^(r+4)))/(r+4))-(M1*exp(x*nu)*(x^(r+3)))+((rho_L*(x^(r+4)))/((r+4)*(exp(nu)-1)))-((rho_r/(exp(nu)-1))*((x^(r+4))/(r+4)))+((M2*exp(x*nu)*(x^(r+3))));
H1=(I1*(r+3)+I2*x-((lambda^2)/6)*I3+((lambda^2)/2)*I4*(x^2)-((lambda^2)/2)*I5*x-I6)*(x^m);
H1_int= int(H1,x,0,1);
H1_int= vpa(H1_int);
H2 = ((J1*(lambda^2))/(6*(m+4)))-((J2*(lambda^2))/(2*(m+3)));
H2=vpa(H2);
G= (E_L/(r+m+1))-M13+((E_L/(exp(nu)-1))*(1/(r+m+1)))+M14;
G=vpa(G);
A(m,r)=H1_int+H2+G;
end
end
A
I have three symbolic terms since I want equations interms of lambda and x. The value of x inserted when I get I's values after r values and m values. I am creating a matrix A with equations having lambda then I take determinant of it to get an equation in terms of lambda and I will use that equation to get value of lambda. I am also sharng what I mean by matrix A. The matrix will have elements like following ;
A11=583418350034.71711450099633761995 - 404.2522486820071715581475695452*lambda^2; %r=0 m=0
A12=342955100491.05339629025357860931 - 236.16385981512705886473948261464*lambda^2; %r=1 m=0
A13=265052784655.08577750712854323127 - 168.18177781842519452660605420978*lambda^2;
A14=203322762437.79122083006315019118 - 137.99080394103775945839083607247*lambda^2;
A15=179379868828.76827703612092258873 - 108.01397559315213169193183004632*lambda^2;
A21=370495028698.81627008913055964828 - 272.77098597891027337922111086527*lambda^2;
A22=228100010500.85073852655079716474 - 172.36591762980006298921802813592*lambda^2;
A23=174909572535.14489241114607164103 - 125.31310187054420562351504054388*lambda^2;
A24=130525707221.65027901348657766572 - 103.54975334547425216012368068556*lambda^2;
A25=110410873567.05769598554052019297 - 81.283292521249163213494102067785*lambda^2;
A31=318836918206.42606727416420657554 - 205.00780517033281356081424897212*lambda^2;
A32=228116511483.27211950393951440909 - 135.3394493480940806583611823317*lambda^2;
A33=206373496726.45624130029358463062 - 99.811322967330165385124898047064*lambda^2;
A34=183606236766.48017821939715890692 - 82.968801004892704989067098085304*lambda^2;
A35=188218588289.44473276912635161394 - 65.334491840422903824342987190786*lambda^2;
A41=249285839982.9324301749480792828 - 163.76728864719879912159069265404*lambda^2;
A42=152346533827.00642443852525232634 - 111.25013657653381685476087220797*lambda^2;
A43=95652277632.223435252212147492575 - 82.914382160944203068253893669008*lambda^2;
A44=33728441196.492275381613993607189 - 69.262645398935919927356866449461*lambda^2;
A45=-54.707899833285924851812473702749*lambda^2 - 30313230834.758141502064257607773;
A51=238011371413.10690304954757739067 - 136.11115268606901224149243477624*lambda^2;
A52=203108325234.9890834260110228119 - 94.364231425385059797811493443391*lambda^2;
A53=248481839194.19255791356538082301 - 70.901034522277689237116974329652*lambda^2;
A54=322167870570.16045465551166213521 - 59.469506840450354534300476840716*lambda^2;
A55=489611884237.36508541953085427943 - 47.10375093128415795710589282516*lambda^2;
and Matrix A will be;
[A11 A12 A13 A14 A15;
A21 A22 A23 A24 A25;
A31 A32 A33 A34 A35;
A41 A42 A43 A44 A45;
A51 A52 A53 A54 A55];
The I take determinant of it once I have matrix A. I want to run for loop since changing r values will give columns values and changing m values will give row values. but I am having the erros
Error in Intgration (line 71)
A(m,r)=H1_int+H2+G;
Caused by:
Error using symengine
Unable to convert expression containing symbolic variables into double array. Apply 'subs' function first to substitute values
for variables.
3 Comments
Dyuman Joshi
on 25 Oct 2023
Part #2
Also, you can parameterize the expression of the variables using symsum() instead of repeatively calculating them each iteration.
Placing independent operations outside for loops is a good technique for improving performance.
syms r k
rho = [rho_L rho_r];
%for M1 to M8
%provide the value of k accordingly
f1(k,r) = (rho(2-rem(k,2))/(exp(nu)-1))*(symsum(((-1)^(r+k-i))*((factorial(r+k))/(factorial(i)*nu^(r-i+k+1))),i,0,r+k));
%for M9 to M12
%same common expression, multiply the extra term accordingly
f2(r) = (1/(exp(nu)-1))*(symsum(((-1)^(r-i))*((factorial(r))/(factorial(i)*nu^(r-i+1))),i,0,r));
You can use the above examples to parameterize for the other variables.
Another thing, you can directly add H1_int, H2 and G and assign the value to A(m,r), then use vpa() on A outside of the loop.
Accepted Answer
Star Strider
on 25 Oct 2023
I am not exactly certain what the problem is, however creating A as a sym matrix and re-defining the loop indices so that they are valid MATLAB indices works and does not throw any errors —
syms lambda
syms x
syms i
R=4;
rho_L= 2702;
rho_r= 5700;
E_L= 70e+9;
E_r= 200e+9;
nu=3;
% A=zeros(R+1,R+1);
A = sym('A',[R+1 R+1]);
for ri=1:R+1
r = ri-1;
for mi=1:R+1
m = mi-1;
M1= (rho_L/(exp(nu)-1))*(symsum(((-1)^(r-i))*((factorial(r))/(factorial(i)*nu^(r-i+1))),i,0,r));
M2= (rho_r/(exp(nu)-1))*(symsum(((-1)^(r-i))*((factorial(r))/(factorial(i)*nu^(r-i+1))),i,0,r));
M3= (rho_L/(exp(nu)-1))*(symsum(((-1)^(r+1-i))*((factorial(r+1))/(factorial(i)*nu^(r-i+2))),i,0,r+1));
M4= (rho_r/(exp(nu)-1))*(symsum(((-1)^(r+1-i))*((factorial(r+1))/(factorial(i)*nu^(r-i+2))),i,0,r+1));
M5= (rho_L/(exp(nu)-1))*(symsum(((-1)^(r+2-i))*((factorial(r+2))/(factorial(i)*nu^(r-i+3))),i,0,r+2));
M6= (rho_r/(exp(nu)-1))*(symsum(((-1)^(r+2-i))*((factorial(r+2))/(factorial(i)*nu^(r-i+3))),i,0,r+2));
M7= (rho_L/(exp(nu)-1))*(symsum(((-1)^(r+3-i))*((factorial(r+3))/(factorial(i)*nu^(r-i+4))),i,0,r+3));
M8= (rho_r/(exp(nu)-1))*(symsum(((-1)^(r+3-i))*((factorial(r+3))/(factorial(i)*nu^(r-i+4))),i,0,r+3));
M9= ((E_L*nu)/(exp(nu)-1))*(symsum(((-1)^(r-i))*((factorial(r))/(factorial(i)*nu^(r-i+1))),i,0,r));
M10= ((E_r*nu)/(exp(nu)-1))*(symsum(((-1)^(r-i))*((factorial(r))/(factorial(i)*nu^(r-i+1))),i,0,r));
M11= ((E_L*nu^2)/(exp(nu)-1))*(symsum(((-1)^(r-i))*((factorial(r))/(factorial(i)*nu^(r-i+1))),i,0,r));
M12= ((E_r*nu^2)/(exp(nu)-1))*(symsum(((-1)^(r-i))*((factorial(r))/(factorial(i)*nu^(r-i+1))),i,0,r));
M13= ((E_L*exp(nu))/(exp(nu)-1))*(symsum(((-1)^(r+m-i))*((factorial(r+m))/(factorial(i)*nu^(r+m-i+1))),i,0,r));
M14= ((E_r*exp(nu))/(exp(nu)-1))*(symsum(((-1)^(r+m-i))*((factorial(r+m))/(factorial(i)*nu^(r+m-i+1))),i,0,r));
J1= ((rho_L/(r+1))*(1+(1/(exp(nu)-1))))-((M1-M2)*exp(nu))-((rho_r/(exp(nu)-1))*(1/(r+1)));
J2= ((rho_L/(r+2))*(1+(1/(exp(nu)-1))))-((M3-M4)*exp(nu))-((rho_r/(exp(nu)-1))*(1/(r+2)));
I1=(-M9*exp(x*nu)*(x^r))+(M10*exp(x*nu)*(x^r));
I2=(M11*exp(x*nu)*(x^r))+(M12*exp(x*nu)*(x^r));
I3=((rho_L*(x^(r+1)))/(r+1))-(M1*exp(x*nu)*(x^r))+((rho_L*(x^(r+1)))/((r+1)*(exp(nu)-1)))-((rho_r/(exp(nu)-1))*((x^(r+1))/(r+1)))+((M2*exp(x*nu)*(x^r)));
I4=((rho_L*(x^(r+2)))/(r+2))-(M1*exp(x*nu)*(x^(r+1)))+((rho_L*(x^(r+2)))/((r+2)*(exp(nu)-1)))-((rho_r/(exp(nu)-1))*((x^(r+2))/(r+2)))+((M2*exp(x*nu)*(x^(r+1))));
I5=((rho_L*(x^(r+3)))/(r+3))-(M1*exp(x*nu)*(x^(r+2)))+((rho_L*(x^(r+3)))/((r+3)*(exp(nu)-1)))-((rho_r/(exp(nu)-1))*((x^(r+3))/(r+3)))+((M2*exp(x*nu)*(x^(r+2))));
I6=((rho_L*(x^(r+4)))/(r+4))-(M1*exp(x*nu)*(x^(r+3)))+((rho_L*(x^(r+4)))/((r+4)*(exp(nu)-1)))-((rho_r/(exp(nu)-1))*((x^(r+4))/(r+4)))+((M2*exp(x*nu)*(x^(r+3))));
H1=(I1*(r+3)+I2*x-((lambda^2)/6)*I3+((lambda^2)/2)*I4*(x^2)-((lambda^2)/2)*I5*x-I6)*(x^m);
H1_int= int(H1,x,0,1);
H1_int= vpa(H1_int);
H2 = ((J1*(lambda^2))/(6*(m+4)))-((J2*(lambda^2))/(2*(m+3)));
H2=vpa(H2);
G= (E_L/(r+m+1))-M13+((E_L/(exp(nu)-1))*(1/(r+m+1)))+M14;
G=vpa(G);
A(mi,ri)=H1_int+H2+G;
end
end
A
.
2 Comments
Star Strider
on 26 Oct 2023
My pleasure!
If my Answer helped you solve your problem, please Accept it!
.
More Answers (0)
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!