For Loop for Symbolic
3 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
Categories
Find more on Creating and Concatenating Matrices in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!