2

I am unable to evaluate and plot a simple (known) function in sympy within Ipython notebook:

y(x) = ( F / 6EI )( x^3 - 3Lx^2 ) where:

  • F = 10^6
  • E = 200E9
  • I = (1/12)(0.5*1^3)
  • L = 3

I define the expression using symbols() and substitute known values (using subs()). 

import sympy
from sympy import symbols

sympy.init_printing()
from IPython.display import Math, Image
from IPython.display import display


F, E, L, I, x = symbols('F, E, L, I, x')    #Definition of sympy symbols

y = F/(6*E*I) * (x**3 - 3*L*x**2)    #Define beam deflection equation
display(Math("y=" + latex(y)))

y = y.subs({F:10**6, E:200E9, L:3, I:(1/12)*(0.5)*(1**3)})    #Substitute known values to define specific deflection curve
display(Math("y=" + latex(y)))

However, it results in:

y(x) = (infinity)x^3 -(infinity)x^2

This is clearly incorrect, as the coefficients should be rational numbers - this is the well understood cantilevered beam deflection equation. It should evaluate to:

y(x) = 2E-5*x^3 - 1.8E-4*x^2 where:

  • y(x=0) = 0
  • y(x=L) = FL^3 / 3EI.

Why does sympy produce this result and how can I modify my solution to attain the correct solution?

As noted in the comments, the above code works correctly for Python 3 but fails in Python 2.7 (the version I am running).

Thomas K
  • 39,200
  • 7
  • 84
  • 86
OnStrike
  • 748
  • 1
  • 6
  • 22
  • Quantity `y` (before and after `subs`) is the correct one when I run your code. Maybe the issue is with the display functions (which you do not provide and I cannot reproduce) – Stelios Sep 15 '16 at 20:26
  • @Stelios Thanks for your interest. The `display()` function is imported from Ipython - Ive edited to reflect this. However,I still get the same result without using `display()`: You used `print y`? – OnStrike Sep 15 '16 at 20:36
  • Yes, `print(y)` (Python 3) gives the expected result. *However*, in Python 2, `print y` gives the "infinity" result! – Stelios Sep 15 '16 at 20:42
  • @Stelios Thank you! I am using Python 2.7 - this gives me a workaround at least, and ill add that to the question.. – OnStrike Sep 15 '16 at 20:47
  • 1
    You substitute `I:(1/12)*(0.5)*(1**3)`. In Python 2, integer division always gives you an integer (rounding down), so `1/12` is 0, I is 0, **6EI** is 0, and you've got a division by zero. Put `from __future__ import division` at the top to get the Python 3 behaviour (`1/12` giving a float) in Python 2. Or use Sympy's Rational class - http://docs.sympy.org/latest/gotchas.html#python-numbers-vs-sympy-numbers – Thomas K Sep 16 '16 at 08:36
  • @ThomasK Beautiful, thank you. I previously tried `float((1/12)*(0.5*1**3)` with no effect, but I now understand it why that didnt work. Ive opted to use the `Rational()` class. Post that as the answer if youd like. – OnStrike Sep 16 '16 at 15:18

1 Answers1

4

If you are using Python 2, 1/12 results in 0, due to integer division. You can use from __future__ import division, or use Python 3 (I recommend using Python 3), and it will return a float.

Better, though, with SymPy, is to use a rational. If you use Rational(1, 12) you will get an exact result. In general, your answers with SymPy will be more accurate if you use exact numbers (rational numbers), and then use evalf() to convert the expression to a float at the end.

asmeurer
  • 86,894
  • 26
  • 169
  • 240