I will appreciate any suggestion on how I could have a solution to this.
Ältere Kommentare anzeigen
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 Kommentare
Isaac
am 12 Aug. 2015
Isaac
am 12 Aug. 2015
Isaac
am 12 Aug. 2015
Walter Roberson
am 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
am 12 Aug. 2015
Should
SXY = 1258[1 2 3];
read
SXY = 1258*[1 2 3];
Walter Roberson
am 13 Aug. 2015
C0(j) has not been assigned a meaning.
Isaac
am 13 Aug. 2015
Isaac
am 13 Aug. 2015
Antworten (3)
Walter Roberson
am 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 Kommentare
Walter Roberson
am 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
am 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
am 13 Aug. 2015
0 Stimmen
Kategorien
Mehr zu Ordinary Differential Equations finden Sie in Hilfe-Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!