0

I am trying to just do some stress calculations that match an example problem I found in a textbook.

The eigenvalue method for a sympy Matrix is giving me an answer that I am not expecting. It is including imaginary numbers. Not sure what I am doing wrong. ANy help would be greatly appreciated

I have looked through the documentation, tried rounding it, or seeing if there were any additional keywords I was missing

from sympy import Matrix
def stressMatrix(x=0,y=0,z=0,xy=0,xz=0,yz=0):

    sigma_x, sigma_y, sigma_z, tau_xy, tau_xz, tau_yz =sympy.var('sigma_x, sigma_y, sigma_z, tau_xy, tau_xz, tau_yz')
    stressVars = {sigma_x: x, sigma_y: y, sigma_z: z, tau_xy:xy, tau_xz: xz, tau_yz: yz}
    stressTensor = sympy.Matrix([[sigma_x,tau_xy,tau_xz],[tau_xy,sigma_y,tau_yz],[tau_xz,tau_yz,sigma_z]])

    return stressTensor.subs(stressVars)

t = stressMatrix(x=-19,y=4.6,z=-8.3,xy=-4.7,xz=6.45,yz=11.8)

eigenValues = t.eigenvals()

I expect an output of 11.618, -9.001, -25.316

But I am getting:

{-227/30 + 411091/(3600*(52767683/216000 + sqrt(7409803141670898)*I/72000)**(1/3)) + (52767683/216000 + sqrt(7409803141670898)*I/72000)**(1/3): 1, -227/30 + 411091/(3600*(-1/2 + sqrt(3)*I/2)*(52767683/216000 + sqrt(7409803141670898)*I/72000)**(1/3)) + (-1/2 + sqrt(3)*I/2)*(52767683/216000 + sqrt(7409803141670898)*I/72000)**(1/3): 1, -227/30 + (-1/2 - sqrt(3)*I/2)*(52767683/216000 + sqrt(7409803141670898)*I/72000)**(1/3) + 411091/(3600*(-1/2 - sqrt(3)*I/2)*(52767683/216000 + sqrt(7409803141670898)*I/72000)**(1/3)): 1}
smichr
  • 16,948
  • 2
  • 27
  • 34

1 Answers1

2

These eigenvalues are real. You can see that if you evaluate them numerically:

In [25]: e1, e2, e3 = t.eigenvals()                                                                                               

In [26]: e1.evalf()                                                                                                               
Out[26]: -25.3162883162696 - 0.e-22⋅ⅈ

In [27]: e2.evalf()                                                                                                               
Out[27]: 11.6178013689878 - 0.e-19⋅ⅈ

In [28]: e3.evalf()                                                                                                               
Out[28]: -9.00151305271824 + 0.e-22⋅ⅈ

SymPy otherwise shows the roots using expressions that involve complex numbers:

In [29]: e1                                                                                                                       
Out[29]: 
                                                                                  ________________________________
  227                           411091                           ⎛  1   √3⋅ⅈ⎞    ╱ 52767683   √7409803141670898⋅ⅈ 
- ─── + ────────────────────────────────────────────────────── + ⎜- ─ + ────⎟⋅3 ╱  ──────── + ─────────────────── 
   30                         ________________________________   ⎝  2    2  ⎠ ╲╱    216000           72000        
             ⎛  1   √3⋅ⅈ⎞    ╱ 52767683   √7409803141670898⋅ⅈ                                                     
        3600⋅⎜- ─ + ────⎟⋅3 ╱  ──────── + ───────────────────                                                     
             ⎝  2    2  ⎠ ╲╱    216000           72000  

This is because there is no other way to express these roots using radicals due to Casus irreducibilis: https://en.wikipedia.org/wiki/Casus_irreducibilis

Oscar Benjamin
  • 12,649
  • 1
  • 12
  • 14
  • `e1.evalf(chop=True)` could be used to [discard small imaginary parts](https://stackoverflow.com/a/53432111/190597). – unutbu Sep 11 '19 at 00:34