python - Sympy, I can't solve an equation symbolically. It keeps trying to give real solutions - Stack Overflow

admin管理员组

文章数量:1391969

import sympy as sym

# Python 3.13
# Used variables in question
x = sym.Symbol("x")
I1 = sym.Symbol("I₁")
L = sym.Symbol("L")
pi = sym.Symbol("Π")
v2 = sym.Symbol("v2")
P = sym.Symbol("P")
E = sym.Symbol("E")

I = I1 * (1 - x/L)

# For the part a of the question one in Assignment 3.
question_part = "a"
if question_part == "a":
    func_v = (1 - sym.cos((pi*x)/(2*L))) * v2
    moment = P * L
    second_derivative_v = sym.diff(func_v, x, x)
    internal_strain_e_func = 0.5 * moment *second_derivative_v / (E * I)
    internal_strain_e = sym.integrate(internal_strain_e_func, (x, 0, L))
    external_e = P * v2
    tip_deflection_eq = sym.Eq(internal_strain_e, external_e)
    tip_deflection = sym.solve(tip_deflection_eq, v2, exclude=[I1, L, pi, P, E], implicit=True)
    sym.pprint(tip_deflection)

Code above is for virtual displacement method for a beam deflection question. I can technically do it by hand however, I get a lot of questions that requires for me to solve integrals and I do make mistakes when I do it by hand. So in this one I thought I would do it by sympy and its okey up until the solve() method. For some reason I can't get a symbolically written solution for eq. The eq. I am trying to solve is:

Edit: I am sorry everyone v2 is not supposed be a part in function v. The question gave it to me and I took it without thinking about it although I did learn stuff that I didn't know and would probably ask later so not much of a waste of your time. Thanks!

import sympy as sym

# Python 3.13
# Used variables in question
x = sym.Symbol("x")
I1 = sym.Symbol("I₁")
L = sym.Symbol("L")
pi = sym.Symbol("Π")
v2 = sym.Symbol("v2")
P = sym.Symbol("P")
E = sym.Symbol("E")

I = I1 * (1 - x/L)

# For the part a of the question one in Assignment 3.
question_part = "a"
if question_part == "a":
    func_v = (1 - sym.cos((pi*x)/(2*L))) * v2
    moment = P * L
    second_derivative_v = sym.diff(func_v, x, x)
    internal_strain_e_func = 0.5 * moment *second_derivative_v / (E * I)
    internal_strain_e = sym.integrate(internal_strain_e_func, (x, 0, L))
    external_e = P * v2
    tip_deflection_eq = sym.Eq(internal_strain_e, external_e)
    tip_deflection = sym.solve(tip_deflection_eq, v2, exclude=[I1, L, pi, P, E], implicit=True)
    sym.pprint(tip_deflection)

Code above is for virtual displacement method for a beam deflection question. I can technically do it by hand however, I get a lot of questions that requires for me to solve integrals and I do make mistakes when I do it by hand. So in this one I thought I would do it by sympy and its okey up until the solve() method. For some reason I can't get a symbolically written solution for eq. The eq. I am trying to solve is:

Edit: I am sorry everyone v2 is not supposed be a part in function v. The question gave it to me and I took it without thinking about it although I did learn stuff that I didn't know and would probably ask later so not much of a waste of your time. Thanks!

• python
• sympy
• equation

Share Improve this question Follow edited Mar 12 at 18:18 President James K. Polk 42.1k2929 gold badges109109 silver badges145145 bronze badges asked Mar 12 at 11:48 Hür Doğan ÜNLÜHür Doğan ÜNLÜ 5544 bronze badges 4

• Well, L-x is L(1-x/L), so remove a factor of 1-x/L from top and bottom. After which, the remaining cosine integral is easily done by hand. – lastchance Commented Mar 12 at 12:05
• Yes I know it can be done by hand but i will need to a lot more of these types of integrals and I would like to learn how can I do them by code. – Hür Doğan ÜNLÜ Commented Mar 12 at 12:12
• 2 You might want to check again how you are building your equation, because the code is actually building a different equation than the one shown on your image. – Davide_sd Commented Mar 12 at 12:15
• Can you clarify your load distribution, please: it's difficult to infer from your code. Also, if your second moment of area is really given by I1(1-x/L) then it will be zero at the end of the beam and (a) it will be very unusual; (b) it will potentially cause problems in the denominator of the integrand. – lastchance Commented Mar 12 at 21:07

Add a comment  | 

2 Answers 2

Sorted by: Reset to default 3

Don't declare a symbol called pi. You should use SymPy's sympy.pi if you mean the mathematical number pi.

This is your equation as stated:

In [42]: x, L, v2, E, I1, P, v2 = symbols('x, L, v2, E, I1, P, v2')

In [43]: eq = Eq(Integral((P*(L - x)*(pi/(2*L))**2*cos((pi*x)/(2*L)) * v2)/(E*I1*(1 - x/L)), (x, 0, L
    ...: )), P*v2)

In [44]: eq
Out[44]: 
L                                   
⌠                                   
⎮  2                 ⎛π⋅x⎞          
⎮ π ⋅P⋅v₂⋅(L - x)⋅cos⎜───⎟          
⎮                    ⎝2⋅L⎠          
⎮ ──────────────────────── dx = P⋅v₂
⎮            2 ⎛    x⎞              
⎮    4⋅E⋅I₁⋅L ⋅⎜1 - ─⎟              
⎮              ⎝    L⎠              
⌡                                   
0  

We can use doit to evaluate the integral and then solve for v2:

In [45]: eq.doit()
Out[45]: 
π⋅P⋅v₂       
────── = P⋅v₂
2⋅E⋅I₁       

In [46]: solve(eq.doit(), v2)
Out[46]: [0]

Apparently this is not the answer you were looking for but in general only v2 = 0 can solve the equation for arbitrary values of the parameters.

Maybe you wanted to solve for something else rather than v2?

This is the equation rearranged:

In [47]: eq2 = eq.doit()

In [48]: eq3 = Eq(factor(eq2.lhs - eq2.rhs), 0)

In [49]: eq3
Out[49]: 
-P⋅v₂⋅(2⋅E⋅I₁ - π)     
─────────────────── = 0
      2⋅E⋅I₁ 

At least one of the factors in the numerator on the lhs must be zero to satisfy the equation. If you want a value of v2 that can satisfy this then in general it has to be v2 = 0.

Maybe the equation has not been put together correctly?

Here you go: cantilever beam with UDL w and tip load P, so that you can check your answer against standard results. But you can modify it to whatever you like if you supply the appropriate loading distribution w(x) and second moment of area I(x).

For constant w and P tip deflection would be the combination of deflections due to UDL: wL^4/8EI and to tip load: PL^3/3EI

import sympy

x = sympy.symbols( 'x' )                         # distance along beam
E = sympy.symbols( 'E' )                         # Young's modulus
I = sympy.symbols( 'I' )                         # second moment of area
L = sympy.symbols( 'L' )                         # length of beam
w = sympy.symbols( 'w' )                         # load per unit length
P = sympy.symbols( 'P' )                         # tip load
S = sympy.symbols( 'S' )                         # shear stress (+ve upward to right of cut)
M = sympy.symbols( 'M' )                         # bending moment (+ve anticlockwise to right of cut)

#w = 0    # special case
#P = 0    # special case

S0 = -sympy.integrate( w, ( x, 0, L ) ) - P                              # S at x=0
M0 = -sympy.integrate( w * x, ( x, 0, L ) ) - P * L                      # M at x=0

S = sympy.integrate( w, ( x, 0, x ) ) + S0                               # Shear stress, from dS/dx = w
M = sympy.integrate( -S, ( x, 0, x ) ) + M0                              # Bending moment, from dM/dx = -S

m = x - L                                        # moment distribution due to unit force at tip

tip_deflection = sympy.integrate( M * m / ( E * I ), ( x, 0, L ) )       # Change in energy integral
tip_deflection = sympy.simplify( tip_deflection )                        # Looks horrible otherwise
print(tip_deflection)

Output:

L**3*(L*w/8 + P/3)/(E*I)

本文标签：