Ciao a tutti!
Sto cercando di calcolare la freccia di un albero di trasmissione mediante il metodo della linea elastica. In parole povere, questo metodo si basa sulla risoluzione di un sistema di equazioni differenziali del primo e secondo ordine e sfruttando determinate condizioni iniziali.
Per fare ciò sto usando Sympy poichè, per ricavare reazioni vincolari e equazioni delle azioni interne necessarie ai vari calcoli, mi è comodo utilizzare variabili simboliche.
Per risolvere questo problema ho provato con questo codice:


import sympy as sym

v1, v2, v3, v4, v5, r1, r2, r3, r4, r5 = sym.symbols('v1, v2, v3, v4, v5, r1, r2, r3, r4, r5', cls = sym.Function) # deflection function

# parametri vari
M_ab_resultant = sym.sqrt((Mx1_act_ab.args[1].subs(equilibrium_system_slow[0]).subs(shaft_dimensions).subs(transmission_data))**2 + (My1_act_ab.args[1].subs(equilibrium_system_slow[0]).subs(shaft_dimensions).subs(transmission_data))**2)
M_bd_resultant = sym.sqrt((Mx1_act_bd.args[1].subs(equilibrium_system_slow[0]).subs(shaft_dimensions).subs(transmission_data))**2 + (My1_act_bd.args[1].subs(equilibrium_system_slow[0]).subs(shaft_dimensions).subs(transmission_data))**2)

# equazioni diff primo ordine
rotation1_01 = sym.Eq(young_modulus*moment_inertia01*sym.Derivative(r1(xi), xi, 1), M_ab_resultant)
rotation1_12 = sym.Eq(young_modulus*moment_inertia12*sym.Derivative(r2(xi), xi, 1), M_bd_resultant)
rotation1_23 = sym.Eq(young_modulus*moment_inertia23*sym.Derivative(r3(xi), xi, 1), M_bd_resultant)
rotation1_34 = sym.Eq(young_modulus*moment_inertia34*sym.Derivative(r4(xi), xi, 1), M_bd_resultant)
rotation1_45 = sym.Eq(young_modulus*moment_inertia45*sym.Derivative(r5(xi), xi, 1), M_bd_resultant)

# equazioni diff secondo ordine
deflection1_01 = sym.Eq(young_modulus*moment_inertia01*sym.Derivative(v1(xi), xi, 2), M_ab_resultant)
deflection1_12 = sym.Eq(young_modulus*moment_inertia12*sym.Derivative(v2(xi), xi, 2), M_bd_resultant)
deflection1_23 = sym.Eq(young_modulus*moment_inertia23*sym.Derivative(v3(xi), xi, 2), M_bd_resultant)
deflection1_34 = sym.Eq(young_modulus*moment_inertia34*sym.Derivative(v4(xi), xi, 2), M_bd_resultant)
deflection1_45 = sym.Eq(young_modulus*moment_inertia45*sym.Derivative(v5(xi), xi, 2), M_bd_resultant)

# lista equazioni
ODE_system = [rotation1_01, rotation1_12, rotation1_23, rotation1_34, rotation1_45, deflection1_01, deflection1_12, deflection1_23, deflection1_34, deflection1_45]

# altri parametri..
dist_1 = bearings_data[b1_width]*0.5
dist_2 = dist_1 + 6 + 0.5*gear_general_data[g1_width]
dist_3 = dist_2 + (shaft_dimensions[d2] - gear_general_data[g2_width])
dist_4 = dist_3 + 6 + gear_general_data[g2_width]
dist_5 = dist_4 + bearings_data[b2_width]*0.5

# condizioni al contorno
ics = {v1(0):0, v1(dist_1) : v2(dist_1), v2(dist_2):v3(dist_2), v3(dist_3):v4(dist_3), v4(dist_4):v5(dist_4), v5(dist_5):0, r1(dist_1) : r2(dist_1), r2(dist_2):r3(dist_2), r3(dist_3):r4(dist_3), r4(dist_4):r5(dist_4)}

# risoluzione sistema
elastic_line_conf_1 = sym.dsolve(ODE_system, [r1(xi), r2(xi), r3(xi), r4(xi), r5(xi), v1(xi), v2(xi), v3(xi), v4(xi), v5(xi)], ics=ics)


Ovviamente non va, e mi mostra il seguente errore:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-86-26e73208c1d3> in <module>
     27 ics = {v1(0):0, v1(dist_1) : v2(dist_1), v2(dist_2):v3(dist_2), v3(dist_3):v4(dist_3), v4(dist_4):v5(dist_4), v5(dist_5):0, r1(dist_1) : r2(dist_1), r2(dist_2):r3(dist_2), r3(dist_3):r4(dist_3), r4(dist_4):r5(dist_4)}
     28
---> 29 elastic_line_conf_1 = sym.dsolve(ODE_system, [r1(xi), r2(xi), r3(xi), r4(xi), r5(xi), v1(xi), v2(xi), v3(xi), v4(xi), v5(xi)], ics=ics)

c:\python37\lib\site-packages\sympy\solvers\ode.py in dsolve(eq, func, hint, simplify, ics, xi, eta, x0, n, **kwargs)
    599         match['eq'] = eq
    600         if len(set(order.values()))!=1:
--> 601             raise ValueError("It solves only those systems of equations whose orders are equal")
    602         match['order'] = list(order.values())[0]
    603         def recur_len(l):

ValueError: It solves only those systems of equations whose orders are equal


Il messaggio è chiaro, purtroppo.
Cosa ho sbagliato? E' possibile utilizzare le comodità del risolutore di sympy, dsolve? (nella documentazione non ho trovato una risposta)
C'è un'altra via? Vorrei evitare Scipy e tenermi le variabili in simbolico.

P.S. : se servono chiarimenti o ulteriori parti di codice chiedete pure!