I will appreciate any suggestion on how I could have a solution to this.
Mostrar comentarios más antiguos
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 comentarios
Isaac
el 12 de Ag. de 2015
Isaac
el 12 de Ag. de 2015
Isaac
el 12 de Ag. de 2015
Walter Roberson
el 12 de Ag. de 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
el 12 de Ag. de 2015
Should
SXY = 1258[1 2 3];
read
SXY = 1258*[1 2 3];
Walter Roberson
el 13 de Ag. de 2015
C0(j) has not been assigned a meaning.
Isaac
el 13 de Ag. de 2015
Isaac
el 13 de Ag. de 2015
Respuestas (3)
Walter Roberson
el 13 de Ag. de 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 comentarios
Walter Roberson
el 13 de Ag. de 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
el 14 de Ag. de 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
el 13 de Ag. de 2015
0 votos
Categorías
Más información sobre Ordinary Differential Equations en Centro de ayuda y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!