追記

Jupyter Notebookを作っておいたので参考にしてください。
https://nbviewer.jupyter.org/gist/tehen-ni-ishi/625e749cb3bbd9cc50c2cb4594d23666

概要

sympy1.1で`sympy.physics.continuum_mechanics.beam`なる梁問題のソルバーが追加されていたらしいです。
Release Notes for 1.1 · sympy/sympy Wiki · GitHub

• Implementation of Singularity Functions to solve Beam Bending problems (Sampad Kumar Saha's GSoC project).
• ...
• Support for Singularity Functions and 2D Beam Problem Solver.

Continuum Mechanics — SymPy 1.1.1 documentation

そこでたわみ角とたわみに境界条件をいい感じにキメて不静定な問題も解けるソルバーが欲しかったので、作りました。

実装

statically_indeterminate_beam.py
```from sympy.core import Symbol, Eq
from sympy.integrals import integrate
from sympy.series import limit
from sympy.solvers import linsolve

from sympy.physics.continuum_mechanics.beam import Beam

def solve_for_reaction_loads(self: Beam, *reactions: Symbol) -> None:
x = self.variable
l = self.length
shear_curve = limit(self.shear_force(), x, l)
moment_curve = limit(self.bending_moment(), x, l)

C1 = Symbol('C1')
slope_curve = integrate(self.bending_moment(), x) + C1

C2 = Symbol('C2')
deflection_curve = integrate(slope_curve, x) + C2

bc_eqs = []
for position, value in self.boundary_conditions['slope']:
eqs = Eq(slope_curve.subs(x, position), value)
bc_eqs.append(eqs)
for position, value in self.boundary_conditions['deflection']:
eqs = Eq(deflection_curve.subs(x, position), value)
bc_eqs.append(eqs)

equations = [Eq(shear_curve), Eq(moment_curve), *bc_eqs]
solutions = next(iter(linsolve(equations, [*reactions, C1, C2])))
self._reaction_loads = dict(zip(reactions, solutions))

# solution = solve([Eq(shear_curve), Eq(moment_curve), *bc_eqs], [*reactions, C1, C2])
# self._reaction_loads = {symbol: solution[symbol] for symbol in reactions}

def _main():
from sympy.core import symbols
from sympy.plotting import plot

l = symbols('l', positive=True, real=True)
RA, RB, RC, RE, RF = symbols('R_A, R_B, R_C, R_E, R_F')
MA, MF = symbols('M_A, M_F')
P = symbols('P')
E, I = symbols('E, I')
(-RA, 0, -1),
(-RB, l, -1),
(-RC, (1 + 2) * l, -1),
(-RE, (1 + 2 + 2 + 2) * l, -1),
(-RF, (1 + 2 + 2 + 2 + 3) * l, -1),
(P, (1 + 2 + 2) * l, -1),
(MA, 0, -2),
(MF, 10 * l, -2),
]

b = Beam(10 * l, E, I)
b.bc_deflection = [
(0, 0),
(l, 0),
((1 + 2) * l, 0),
((1 + 2 + 2 + 2) * l, 0),
((1 + 2 + 2 + 2 + 3) * l, 0),
]
b.bc_slope = [
(0, 0),
((1 + 2 + 2 + 2 + 3) * l, 0),
]

b.solve_for_reaction_loads(RA, RB, RC, RE, RF, MA, MF)

# こことんでもなく雑な数字
constants = {l: 1, P: 1, E: 1, I: 1}
plot(
# b.shear_force().subs(constants),
b.bending_moment().subs(constants),
# b.slope().subs(constants),
# b.deflection().subs(constants),
(b.variable, 0, b.length.subs(constants))
)

if __name__ == '__main__':
_main()

```

あとがき

いろいろ問題を入れて元の`solve_for_reaction_loads`の実装との互換性や、不静定梁の問題の整合性はある程度確認はしていますが、万が一間違った結果が出て単位を落としたり設計したものが破損しても免責とさせて頂きます。