Sympy solve函数求解梁力学方程组得错误解，不满足约束

我尝试用Sympy求解特殊荷载作用下梁的力学行为方程组，运行solve求解未知量时没有错误警告，但将求得的变量代回后，结果不满足设定的方程组。绘制y1、y2、y3、y4的挠度曲线时也出现异常，检查发现y2在x=10处、y4在x=20处不满足约束条件。我不清楚原因，因为求解后几乎没做额外操作。

相关代码如下：

import numpy as np
from sympy import *
#lots of setup
E = 30*10**9 
I = 5.5*10**9/(1000**4)
w1 = 30000
p1 = 200000
ma = 90000
p2 = 15000
x,r1,m,r2,r3,c1,c2,c3,c4,c5,c6,c7,c8 = symbols('x,r1,m,r2,r3,c1,c2,c3,c4,c5,c6,c7,c8')
m1 = -m+r1*x-(-w1*x**3/36+w1*x**2/2)
t1p = integrate(m1,x)+c1
y1p = integrate(t1p,x)+c2
t1 = t1p/E/I/1.5
y1 = y1p/E/I/1.5
m2 = -m+r1*x-(w1*6/2)*(x-2)-p1*(x-6)
t2p = integrate(m2,x)+c3
y2p = integrate(t2p,x)+c4
t2 = t2p/E/I
y2 = y2p/E/I
m3 = m2+r2*(x-10)+ma
t3p = integrate(m3,x)+c5
y3p = integrate(t3p,x)+c6
t3 = t3p/E/I
y3 = y3p/E/I
m4 = m3-p2*(x-15)
t4p = integrate(m4,x)+c7
y4p = integrate(t4p,x)+c8
t4 = t4p/E/I
y4 = y4p/E/I

#Equations to be solved

eq1 = Eq(w1*6/2-p1-p2+r1+r2+r3,0)
eq2 = Eq(m-w1*6/2*2-p1*6-ma-p2*15+r2*10+r3*20,0)
eq3 = Eq(y1.subs(x,0),0)
eq4 = Eq(t1.subs(x,0),0)
eq5 = Eq(y2.subs(x,10),0)
eq6 = Eq(y4.subs(x,20),0)
eq7 = Eq(y1.subs(x,6),y2.subs(x,6))
eq8 = Eq(t1.subs(x,6),t2.subs(x,6))
eq9 = Eq(y2.subs(x,10),y3.subs(x,10))
eq10 = Eq(t2.subs(x,10),t3.subs(x,10))
eq11 = Eq(y3.subs(x,15),y4.subs(x,15))
eq12 = Eq(t3.subs(x,15),t4.subs(x,15))
sol = solve((eq1,eq2,eq3,eq4,eq5,eq6,eq7,eq8,eq9,eq10,eq11,eq12),(r1,m,r2,r3,c1,c2,c3,c4,c5,c6,c7,c8))
# print(sol)
# print(y1)
# print(y2)
# print(y3)
# print(y4)
# m1 = m1.subs({m:sol[m], r1:sol[r1]})
# m2 = m2.subs({m:sol[m], r1:sol[r1]})
# m3 = m3.subs({m:sol[m], r1:sol[r1], r2:sol[r2]})
# m4 = m4.subs({m:sol[m], r1:sol[r1], r2:sol[r2]})
# y1 = y1.subs({m:sol[m], r1:sol[r1], c1:sol[c1], c2:sol[c2]})
y2 = y2.subs({m:sol[m], r1:sol[r1], c3:sol[c3], c4:sol[c4]})
# y3 = y3.subs({m:sol[m], r1:sol[r1], r2:sol[r2], c5:sol[c5], c6:sol[c6]})
y4 = y4.subs({m:sol[m], r1:sol[r1], r2:sol[r2], c7:sol[c7], c8:sol[c8]})

#check parameters are satisfied
print(y2.subs(x,10),y4.subs(x,20))

内容的提问来源于stack exchange，提问作者Nicholas Marchak