I will appreciate any suggestion on how I could have a solution to this.
Show older comments
SX=1000*[1 2 3];
SY=2000*[1.5 2 3];
SXY = 1258[1 2 3];
a = [0.3 0.6 0.9];
syms rb
for j=1:1:3
if pwmid(j)<=pwc(j)
SRR(j)=0.5*(SX(j)+SY(j)).*(1-(a(j).^2)/rb^2)+0.5*(SX(j)-SY(j)).*(1+(3*a(j).^4/rb^4)-(4*a(j).^2/rb^2))*cos(2*thbkso(j))...
+SXY(j).*(1+(3*a(j).^4/rb^4)-(4*a(j).^2/rb^2))*sin(2*thbkso(j))+(a(j).^2/rb^2).*pwmid(j);
STT(j)=0.5*(SX(j)+SY(j)).*(1+(a(j).^2)/rb^2)-0.5*(SX(j)-SY(j)).*(1+(3*a(j).^4/rb^4))*cos(2*thbkso(j))...
-SXY(j).*(1+(3*a(j).^4/rb^4))*sin(2*thbkso(j))-(a(j).^2/rb^2).*pwmid(j);
SRT(j)=(0.5*(SX(j)-SY(j)).*sin(2*thbkso(j))+SXY(j).*cos(2*thbkso(j))).*(1-(3*a(j).^4/rb^4)+(2*a(j).^2/rb^2));
SIGMA1A(j)=0.5*(STT(j)+SRR(j))+0.5*((STT(j)-SRR(j)).^2+4*SRT(j).^2).^0.5;
SIGMA3A(j)=0.5*(STT(j)+SRR(j))-0.5*((STT(j)-SRR(j)).^2+4*SRT(j).^2).^0.5;
C0FUN(j)=SIGMA1A(j)-SIGMA3A(j);
rbsoln{j}=double(vpasolve(C0FUN(j)==C0(j),rb));
cell(rbsoln);
rw(j)=rbsoln{j}(1);
rbkt_art(j) = rbsoln{j}(1)-a(j);
else
rw(j)=a(j);
rbkt_art(j)=rbkt_int(j);
end
end
8 Comments
Isaac
on 12 Aug 2015
Isaac
on 12 Aug 2015
Isaac
on 12 Aug 2015
Walter Roberson
on 12 Aug 2015
I do not understand what your question is? vpasolve() is going to give you a solution if it can find one. If it is not the "right" solution then you need some way of characterizing the "right" solution. For example is there a range of values that the "right" solution will always be in?
Walter Roberson
on 12 Aug 2015
Should
SXY = 1258[1 2 3];
read
SXY = 1258*[1 2 3];
Walter Roberson
on 13 Aug 2015
C0(j) has not been assigned a meaning.
Isaac
on 13 Aug 2015
Isaac
on 13 Aug 2015
Answers (3)
Walter Roberson
on 13 Aug 2015
All solutions to those equations are strictly imaginary for the parameters you give.
For example, for j = 1, the solutions are
(3/10)*sqrt(5)*sqrt(roots([+8105,-9500,+4790,-1164,+117]))
and the negatives of those.
2 Comments
Walter Roberson
on 13 Aug 2015
Please explain what you mean when you said you were concerned about solve or vpasolve "not giving favorable results" ?
If you want all of the results, then you may have to use solve() instead of vpasolve(), and you might have to double() the result of solve() to get numeric values. I do not have the Symbolic Toolbox so I cannot check exactly what would be returned.
Walter Roberson
on 14 Aug 2015
I could have made a mistake along the way, but if I got it right then:
for j = 1 : 3 A = 32 * SXY(j) * sin(thbkso(j)) * cos(thbkso(j))^3 * SX(j) - 32 * SXY(j) * sin(thbkso(j)) * cos(thbkso(j))^3 * SY(j) - 16 * SXY(j) * sin(thbkso(j)) * cos(thbkso(j)) * SX(j) + 16 * SXY(j) * sin(thbkso(j)) * cos(thbkso(j)) * SY(j) - C0(j)^2 + SX(j)^2 - 2 * SX(j) * SY(j) + 4 * SXY(j)^2 + SY(j)^2;
B = - 32 * cos(thbkso(j))^4 * SX(j)^2 * a(j)^2 + 64 * cos(thbkso(j))^4 * SX(j) * SY(j) * a(j)^2 + 128 * cos(thbkso(j))^4 * SXY(j)^2 * a(j)^2 - 32 * cos(thbkso(j))^4 * SY(j)^2 * a(j)^2 - 8 * sin(thbkso(j)) * cos(thbkso(j)) * SX(j) * SXY(j) * a(j)^2 - 8 * sin(thbkso(j)) * cos(thbkso(j)) * SXY(j) * SY(j) * a(j)^2 + 16 * sin(thbkso(j)) * cos(thbkso(j)) * SXY(j) * a(j)^2 * pwmid(j) + 28 * cos(thbkso(j))^2 * SX(j)^2 * a(j)^2 - 64 * cos(thbkso(j))^2 * SX(j) * SY(j) * a(j)^2 + 8 * cos(thbkso(j))^2 * SX(j) * a(j)^2 * pwmid(j) - 128 * cos(thbkso(j))^2 * SXY(j)^2 * a(j)^2 + 36 * cos(thbkso(j))^2 * SY(j)^2 * a(j)^2 - 8 * cos(thbkso(j))^2 * SY(j) * a(j)^2 * pwmid(j) - 2 * SX(j)^2 * a(j)^2 + 8 * SX(j) * SY(j) * a(j)^2 - 4 * SX(j) * a(j)^2 * pwmid(j) + 16 * SXY(j)^2 * a(j)^2 - 6 * SY(j)^2 * a(j)^2 + 4 * SY(j) * a(j)^2 * pwmid(j);
C = 128 * sin(thbkso(j)) * cos(thbkso(j))^3 * SX(j) * SXY(j) * a(j)^4 - 128 * sin(thbkso(j)) * cos(thbkso(j))^3 * SXY(j) * SY(j) * a(j)^4 + 48 * cos(thbkso(j))^4 * SX(j)^2 * a(j)^4 - 96 * cos(thbkso(j))^4 * SX(j) * SY(j) * a(j)^4 - 192 * cos(thbkso(j))^4 * SXY(j)^2 * a(j)^4 + 48 * cos(thbkso(j))^4 * SY(j)^2 * a(j)^4 - 48 * sin(thbkso(j)) * cos(thbkso(j)) * SX(j) * SXY(j) * a(j)^4 + 80 * sin(thbkso(j)) * cos(thbkso(j)) * SXY(j) * SY(j) * a(j)^4 - 32 * sin(thbkso(j)) * cos(thbkso(j)) * SXY(j) * a(j)^4 * pwmid(j) - 40 * cos(thbkso(j))^2 * SX(j)^2 * a(j)^4 + 96 * cos(thbkso(j))^2 * SX(j) * SY(j) * a(j)^4 - 16 * cos(thbkso(j))^2 * SX(j) * a(j)^4 * pwmid(j) + 192 * cos(thbkso(j))^2 * SXY(j)^2 * a(j)^4 - 56 * cos(thbkso(j))^2 * SY(j)^2 * a(j)^4 + 16 * cos(thbkso(j))^2 * SY(j) * a(j)^4 * pwmid(j) + 7 * SX(j)^2 * a(j)^4 - 18 * SX(j) * SY(j) * a(j)^4 + 4 * SX(j) * a(j)^4 * pwmid(j) - 8 * SXY(j)^2 * a(j)^4 + 15 * SY(j)^2 * a(j)^4 - 12 * SY(j) * a(j)^4 * pwmid(j) + 4 * a(j)^4 * pwmid(j)^2;
E = 288 * SXY(j) * a(j)^8 * sin(thbkso(j)) * cos(thbkso(j))^3 * SX(j) - 288 * SXY(j) * a(j)^8 * sin(thbkso(j)) * cos(thbkso(j))^3 * SY(j) - 144 * SXY(j) * a(j)^8 * sin(thbkso(j)) * cos(thbkso(j)) * SX(j) + 144 * SXY(j) * a(j)^8 * sin(thbkso(j)) * cos(thbkso(j)) * SY(j) + 9 * SX(j)^2 * a(j)^8 - 18 * SY(j) * a(j)^8 * SX(j) + 36 * SXY(j)^2 * a(j)^8 + 9 * SY(j)^2 * a(j)^8;
sols_plus = sqrt( roots([A, B, C, 0, E]) );
sols{j} = [sols_plus; -sols_plus];
endI am not certain of these coefficients; I am concerned that the previous solution did not have a 0 in the x^1 position but this does.
Isaac
on 13 Aug 2015
0 votes
Categories
Find more on Ordinary Differential Equations 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!