For Loop for Symbolic

2 views (last 30 days)
Syed Abdul Rafay
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
Unable to perform assignment because value of type 'sym' is not convertible to 'double'.

Caused by:
Error using symengine
Unable to convert expression containing symbolic variables into double array. Apply 'subs' function first to substitute values for variables.
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
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.
Syed Abdul Rafay
Syed Abdul Rafay on 26 Oct 2023
I will try it. Actually I tried a lot of different approaches but it was not giving me int values. I will try your approach to see if it works.

Sign in to comment.

Accepted Answer

Star Strider
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
A = 
.
  2 Comments
Syed Abdul Rafay
Syed Abdul Rafay on 26 Oct 2023
Thanks for sharing. I see where I had the problem. Thanks for your help.
Star Strider
Star Strider on 26 Oct 2023
My pleasure!
If my Answer helped you solve your problem, please Accept it!
.

Sign in to comment.

More Answers (0)

Categories

Find more on MATLAB in Help Center and File Exchange

Tags

Products


Release

R2022b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!