Symabolic variables and the vpa function.
Show older comments
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.
Answers (1)
Walter Roberson
on 3 Mar 2014
eqn= (x^3)-(J3*(x^2))-(J2*x)-J1;
No quotations, and the "=0" removed.
Categories
Find more on Assembly 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!