3.5.B.sagews

from sympy import *
half=S(1)/2

# SI units:
m  = Symbol("m",  positive=True)
s  = Symbol("s",  positive=True)
kg = Symbol("kg", positive=True)


# Given:
g = var("g")
Mass, mass  = var("M, m")
r, theta = var("r, theta")
c, l = var("c, l")
theta=half*M*r*r

# 8 Unknowns:
# x = x
# xpp = x''
# phipp = phi''
# psipp = psi''
x, xpp, phipp, psipp, xc = var("x, ddotx, ddotphi, ddotpsi, x_c")
S1, S2, Fc = var("S1, S2, F_c")

eq1 = Eq(Mass*xpp, S2 - S1 - Mass*g + Fc)
eq2 = Eq(theta*psipp, -S2*r + Fc*r)
eq3 = Eq(4*l*l*mass*phipp, -2*mass*g*l + S1*l)
eq4 = Eq(xc, 2*x)
eq5 = Eq(phipp, xpp/l)
eq6 = Eq(psipp, xpp/r)
eq7 = Eq(Fc, -c*xc)

eqs=[eq1,eq2,eq3,eq4,eq5,eq6,eq7]
unknowns=[x, xpp, phipp, psipp, xc, S1, S2, Fc]

# don't solve for x:
del unknowns[0]

sol = solve(eqs, unknowns, dict=True)

# print solution for x''
# x'' = ...
pprint("x''=")
pprint(sol[0][xpp])

x''= -(2*M*g + 8*c*x + 4*g*m) ------------------------- 3*M + 8*m