Symabolic variables and the vpa function.
Afficher commentaires plus anciens
I am trying to solve for principal stresses by solving the nonlinear equation. I get three solutions, sigma 1, sigma 2, and sigma 3 as a vector like I want. However they are giant long expressions in terms of J1, J2 and J3. How do I get a single number from these. vpa() gets me a shorter equation, but still in terms of the j values. Double does not work because they are symbolic. here is my code:
StressState= [120 61 70 -45 49 53];
sigmax=StressState(1);
sigmay=StressState(2);
sigmaz=StressState(3);
tauxy=StressState(4);
tauyz=StressState(5);
tauzx=StressState(6);
___________________________________________________________
function [stresses] = PrincipalStresses(sigmax,sigmay,sigmaz,tauxy,tauyz,tauzx)
J1=(sigmax*sigmay*sigmaz)+(2*tauxy*tauyz*tauzx)-(sigmax*(tauyz^2))-(sigmay*(tauzx^2))-(sigmaz*(tauxy^2));
J2=(tauxy^2)+(tauyz^2)+(tauzx^2)-(sigmax*sigmay)-(sigmay*sigmaz)-(sigmaz*sigmax);
J3=sigmax+sigmay+sigmaz;
eqn='(x^3)-(J3*(x^2))-(J2*x)-J1=0';
stresses=vpa(solve(eqn,'x'));
__________________________________________________________
[stresses] = PrincipalStresses(sigmax,sigmay,sigmaz,tauxy,tauyz,tauzx);
If I do:
s=stresses(1)
I get
s= 0.33333333333333333333333333333333*J3 + (0.5*J1 + ((0.037037037037037037037037037037037*J3^3 + 0.16666666666666666666666666666667*J2*J3 + 0.5*J1)^2 - 1.0*(0.11111111111111111111111111111111*J3^2 + 0.33333333333333333333333333333333*J2)^3)^(1/2) + 0.037037037037037037037037037037037*J3^3 + 0.16666666666666666666666666666667*J2*J3)^(1/3) + (0.11111111111111111111111111111111*J3^2 + 0.33333333333333333333333333333333*J2)/(0.5*J1 + ((0.037037037037037037037037037037037*J3^3 + 0.16666666666666666666666666666667*J2*J3 + 0.5*J1)^2 - 1.0*(0.11111111111111111111111111111111*J3^2 + 0.33333333333333333333333333333333*J2)^3)^(1/2) + 0.037037037037037037037037037037037*J3^3 + 0.16666666666666666666666666666667*J2*J3)^(1/3)
If I don't use vpa() I get:
s= J3/3 + (J1/2 + J3^3/27 + (J2*J3)/6 + ((J3^3/27 + (J2*J3)/6 + J1/2)^2 - (J3^2/9 + J2/3)^3)^(1/2))^(1/3) + (J3^2/9 + J2/3)/(J1/2 + J3^3/27 + (J2*J3)/6 + ((J3^3/27 + (J2*J3)/6 + J1/2)^2 - (J3^2/9 + J2/3)^3)^(1/2))^(1/3)
What am I missing? I know I should be getting actual numbers but I can only get these as answers. Thanks in advance.
Réponses (1)
Walter Roberson
le 3 Mar 2014
eqn= (x^3)-(J3*(x^2))-(J2*x)-J1;
No quotations, and the "=0" removed.
Catégories
En savoir plus sur Assembly dans Centre d'aide et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!