For plane Euler Bernoulli beams under lateral loading, the unknown displacement function is
. In this case, the potential energy of the system has the following form (See Euler Bernoulli Beam and energy expressions):
View Mathematica Code
(*Exact Solution*)
Clear[q,EI,y,x,c1,c2,c3,c4,L]
y=(q*x^4/24+c1*x^3/6+c2*x^2/2+c3*x+c4)/EI;
Eq1=y/.x-> 0;
Eq2=y/.x-> L;
Eq3=EI*D[y,{x,2}]/.x-> 0;
Eq4=EI*D[y,{x,2}]/.x-> L;
a=Solve[{Eq1==0,Eq2==0,Eq3==0,Eq4==0},{c1,c2,c3,c4}];
y=y/.a[[1]];
theta=D[y,x]/.a[[1]];
M=FullSimplify[EI*D[y,{x,2}]/.a[[1]]];
V=FullSimplify[EI*D[y,{x,3}]/.a[[1]]];
(*Second Degree*)
y2=a2*x^2+a1*x+a0;
s=Solve[{(y2/.x-> 0)==0,(y2/.x->L)==0},{a0,a1}];
y2=y2/.s[[1]];
theta2=D[y2,x];
M2=EI*D[y2,{x,2}];
V2=EI*D[y2,{x,3}];
PE=EI/2*Integrate[(D[y2,{x,2}])^2,{x,0,L}]-Integrate[q*y2,{x,0,L}];
Eq1=D[PE,a2];
a=Solve[{Eq1==0},{a2}];
y2=y2/.a[[1]];
theta2=D[y2,x]/.a[[1]];
M2=FullSimplify[EI*D[y2,{x,2}]/.a[[1]]];
V2=FullSimplify[EI*D[y2,{x,3}]/.a[[1]]];
(*Third Degree*)
y3=a3*x^3+a2*x^2+a1*x+a0;
s=Solve[{(y3/.x-> 0)==0,(y3/.x->L)==0},{a0,a1}];
y3=y3/.s[[1]];
theta3=D[y3,x];
M3=EI*D[y3,{x,2}];
V3=EI*D[y3,{x,3}];
PE=EI/2*Integrate[(D[y3,{x,2}])^2,{x,0,L}]-Integrate[q*y3,{x,0,L}];
Eq1=D[PE,a2];
Eq2=D[PE,a3];
a=Solve[{Eq1==0,Eq2==0},{a2,a3}];
y3=y3/.a[[1]];
theta3=D[y3,x]/.a[[1]];
M3=FullSimplify[EI*D[y3,{x,2}]/.a[[1]]];
V3=FullSimplify[EI*D[y3,{x,3}]/.a[[1]]];
(*fourth Degree*)
y4=a4*x^4+a3*x^3+a2*x^2+a1*x+a0;
s=Solve[{(y4/.x-> 0)==0,(y4/.x->L)==0},{a0,a1}];
y4=y4/.s[[1]];
theta4=D[y4,x];
M4=EI*D[y4,{x,2}];
V4=EI*D[y4,{x,3}];
PE=EI/2*Integrate[(D[y4,{x,2}])^2,{x,0,L}]-Integrate[q*y4,{x,0,L}];
Eq1=D[PE,a2];
Eq2=D[PE,a3];
Eq3=D[PE,a4];
a=Solve[{Eq1==0,Eq2==0,Eq3==0},{a2,a3,a4}];
y4=y4/.a[[1]];
theta4=D[y4,x]/.a[[1]];
M4=FullSimplify[EI*D[y4,{x,2}]/.a[[1]]];
V4=FullSimplify[EI*D[y4,{x,3}]/.a[[1]]];
(*fifth Degree*)
y5=a5*x^5+a4*x^4+a3*x^3+a2*x^2+a1*x+a0;
s=Solve[{(y5/.x-> 0)==0,(y5/.x->L)==0},{a0,a1}];
y5=y5/.s[[1]];
theta4=D[y5,x];
M5=EI*D[y5,{x,2}];
V5=EI*D[y5,{x,3}];
PE=EI/2*Integrate[(D[y5,{x,2}])^2,{x,0,L}]-Integrate[q*y5,{x,0,L}];
Eq1=D[PE,a2];
Eq2=D[PE,a3];
Eq3=D[PE,a4];
Eq4=D[PE,a5];
a=Solve[{Eq1==0,Eq2==0,Eq3==0,Eq4==0},{a2,a3,a4,a5}];
y5=y5/.a[[1]];
theta5=D[y5,x]/.a[[1]];
M5=FullSimplify[EI*D[y5,{x,2}]/.a[[1]]];
V5=FullSimplify[EI*D[y5,{x,3}]/.a[[1]]];
EI = 200000*1000000*4*100000000/1000000000000;
q = -25*1000;
L = 10;
Plot[{y,y2,y4},{x,0,L},AxesLabel-> {"x (m)","y (m)"},PlotLegends->{Style["Exact y",Bold,12],Style["y2",Bold,12],Style["y4",Bold,12]}]
Plot[{theta,theta2,theta4},{x,0,L},AxesLabel-> {"x (m)","theta"},PlotLegends-> {Style["Exact theta",Bold,12],Style["theta 2",Bold,12],Style["theta 4",Bold,12]}]
Plot[{M,M2,M4},{x,0,L},AxesLabel-> {"x (m)","M (N.m.)"},PlotLegends-> {Style["Exact M",Bold,12],Style["M2",Bold,12],Style["M4",Bold,12]}]
Plot[{V,V2,V4},{x,0,L},AxesLabel-> {"x (m)","V (N)"},PlotLegends->{Style["Exact V",Bold,12],Style["V 2",Bold,12],Style["V4",Bold,12]}]
View Python Code
from sympy import *
from numpy import *
import sympy as sp
import numpy as np
import matplotlib.pyplot as plt
sp.init_printing(use_latex = "mathjax")
x, EI, P, q, a0, a1, a2, a3, a4, a5, c1, c2, c3, c4 = symbols("x EI P q a_0 a_1 a_2 a_3 a_4 a_5 c_1 c_2 c_3 c_4")
display("Exact")
y = (q*x**4/24+c1*x**3/6+c2*x**2/2+c3*x+c4)/EI
Eq1 = y.subs(x,0)
Eq2 = y.subs(x,L)
Eq3 = EI * y.diff(x,2).subs(x,0)
Eq4 = EI * y.diff(x,2).subs(x,L)
a = solve((Eq1,Eq2,Eq3,Eq4), (c1,c2,c3,c4))
display("constants: ", a)
y = y.subs({c1:a[c1], c2:a[c2], c3:a[c3], c4:a[c4]})
theta = y.diff(x)
M = EI*y.diff(x,2)
V = EI*y.diff(x,3)
display("y,theta,M,V: ", y,theta,M,V)
display("Second Degree")
y2 = a2*x**2+a1*x+a0
s = solve((y2.subs(x,0), y2.subs(x,L)), (a0, a1))
display("a0, a1", s)
y2 = y2.subs({a0:s[a0], a1:s[a1]})
PE2 = (EI/2)*integrate(y2.diff(x,2)**2, (x,0,L)) - integrate(q*y2, (x,0,L))
display("Potential Energy: ", PE2)
Eq1_2 = PE2.diff(a2)
a = solve(Eq1_2, a2)
y2 = y2.subs(a2,a[0])
theta2 = y2.diff(x)
M2 = EI*y2.diff(x,2)
V2 = EI*y2.diff(x,3)
display("y,theta,M,V: ", y2,theta2,M2,V2)
display("Third Degree")
y3 = a3*x**3+a2*x**2+a1*x+a0
s = solve((y3.subs(x,0), y3.subs(x,L)), (a0, a1))
display("a0, a1", s)
y3 = y3.subs({a0:s[a0], a1:s[a1]})
PE3 = (EI/2)*integrate(y3.diff(x,2)**2, (x,0,L)) - integrate(q*y3, (x,0,L))
display("Potential Energy: ", PE3)
Eq1_3 = PE3.diff(a2)
Eq2_3 = PE3.diff(a3)
a = solve((Eq1_3, Eq2_3), (a2,a3))
y3 = y3.subs({a2:a[a2], a3:a[a3]})
theta3 = y3.diff(x)
M3 = EI*y3.diff(x,2)
V3 = EI*y3.diff(x,3)
display("y,theta,M,V: ", y3,theta3,M3,V3)
display("Fourth Degree")
y4 = a4*x**4+a3*x**3+a2*x**2+a1*x+a0
s = solve((y4.subs(x,0), y4.subs(x,L)), (a0, a1))
display("a0, a1", s)
y4 = y4.subs({a0:s[a0], a1:s[a1]})
PE4 = (EI/2)*integrate(y4.diff(x,2)**2, (x,0,L)) - integrate(q*y4, (x,0,L))
display("Potential Energy: ", PE4)
Eq1_4 = PE4.diff(a2)
Eq2_4 = PE4.diff(a3)
Eq3_4 = PE4.diff(a4)
a = solve((Eq1_4, Eq2_4, Eq3_4), (a2,a3, a4))
y4 = y4.subs({a2:a[a2], a3:a[a3], a4:a[a4]})
theta4 = y4.diff(x)
M4 = EI*y4.diff(x,2)
V4 = EI*y4.diff(x,3)
display("y,theta,M,V: ", y4,theta4,M4,V4)
display("Fifth Degree")
y5 = a5*x**5+a4*x**4+a3*x**3+a2*x**2+a1*x+a0
s = solve((y5.subs(x,0), y5.subs(x,L)), (a0, a1))
display("a0, a1", s)
y5 = y5.subs({a0:s[a0], a1:s[a1]})
PE5 = (EI/2)*integrate(y5.diff(x,2)**2, (x,0,L)) - integrate(q*y5, (x,0,L))
display("Potential Energy: ", PE5)
Eq1_5 = PE5.diff(a2)
Eq2_5 = PE5.diff(a3)
Eq3_5 = PE5.diff(a4)
Eq4_5 = PE5.diff(a5)
a = solve((Eq1_5, Eq2_5, Eq3_5, Eq4_5), (a2,a3,a4,a5))
y5 = y5.subs({a2:a[a2], a3:a[a3], a4:a[a4], a5:a[a5]})
theta5 = y5.diff(x)
M5 = EI*y5.diff(x,2)
V5 = EI*y5.diff(x,3)
display("y,theta,M,V: ", y5,theta5,M5,V5)
# Sub graphed things in
y = y.subs({EI:200000*1000000*4*100000000/1000000000000, L:10, q:-25*100})
y2 = y2.subs({EI:200000*1000000*4*100000000/1000000000000, L:10, q:-25*100})
y4 = y4.subs({EI:200000*1000000*4*100000000/1000000000000, L:10, q:-25*100})
theta = theta.subs({EI:200000*1000000*4*100000000/1000000000000, L:10, q:-25*100})
theta2 = theta2.subs({EI:200000*1000000*4*100000000/1000000000000, L:10, q:-25*100})
theta4 = theta4.subs({EI:200000*1000000*4*100000000/1000000000000, L:10, q:-25*100})
M = M.subs({EI:200000*1000000*4*100000000/1000000000000, L:10, q:-25*100})
M2 = M2.subs({EI:200000*1000000*4*100000000/1000000000000, L:10, q:-25*100})
M4 = M4.subs({EI:200000*1000000*4*100000000/1000000000000, L:10, q:-25*100})
V = V.subs({EI:200000*1000000*4*100000000/1000000000000, L:10, q:-25*100})
V2 = V2.subs({EI:200000*1000000*4*100000000/1000000000000, L:10, q:-25*100})
V4 = V4.subs({EI:200000*1000000*4*100000000/1000000000000, L:10, q:-25*100})
#Graph
fig, ax = plt.subplots(2,2, figsize = (10,6))
plt.setp(ax[0,0], xlabel = "x(m) ", ylabel = "y(m)")
plt.setp(ax[0,1], xlabel = "x(m)", ylabel = "theta")
plt.setp(ax[1,0], xlabel = "x(m)", ylabel = "M(N*m)")
plt.setp(ax[1,1], xlabel = "x(m)", ylabel = "V(N)")
x1 = np.arange(0,10.1,0.1)
y_list = [y.subs({x:i}) for i in x1]
y2_list = [y2.subs({x:i}) for i in x1]
y4_list = [y4.subs({x:i}) for i in x1]
theta_list = [theta.subs({x:i}) for i in x1]
theta2_list = [theta2.subs({x:i}) for i in x1]
theta4_list = [theta4.subs({x:i}) for i in x1]
M_list = [M.subs({x:i}) for i in x1]
M2_list = [M2.subs({x:i}) for i in x1]
M4_list = [M4.subs({x:i}) for i in x1]
V_list = [V.subs({x:i}) for i in x1]
V2_list = [V2.subs({x:i}) for i in x1]
V4_list = [V4.subs({x:i}) for i in x1]
ax[0,0].plot(x1, y_list, 'blue', label = "Exact y")
ax[0,0].plot(x1, y2_list, 'orange', label = "y2")
ax[0,0].plot(x1, y4_list, 'green', label = "y4")
ax[0,0].legend()
ax[0,1].plot(x1, theta_list, 'blue', label = "Exact theta")
ax[0,1].plot(x1, theta2_list, 'orange', label = "theta 2")
ax[0,1].plot(x1, theta4_list, 'green', label = "theta 4")
ax[0,1].legend()
ax[1,0].plot(x1, M_list, 'blue', label = "Exact M")
ax[1,0].plot(x1, M2_list, 'orange', label = "M2")
ax[1,0].plot(x1, M4_list, 'green', label = "M4")
ax[1,0].legend()
ax[1,1].plot(x1, V_list, 'blue', label = "Exact V")
ax[1,1].plot(x1, V2_list, 'orange', label = "V2")
ax[1,1].plot(x1, V4_list, 'green', label = "V4")
ax[1,1].legend()
Exact Solution:
First, the exact solution that would satisfy the equilibrium equations can be obtained. The equilibrium equation as shown in the Euler Bernoulli beam section apply only when
and
are constant. Since
is variable, the differential equation has the form:
Approximate solutions using higher order polynomials can be obtained similarly. The following plots show the convergence of the polynomial functions to the exact solution up to a polynomial of the sixth degree (Mathematica code is given below). In the first plot, it is clear that a polynomial of the second degree is capable of capturing the displacement quite accurately. However, the corresponding approximate solutions for the moment and shear are not as accurate. The accuracy of the approximate solution for the bending moment is higher than that for the shearing force and that is because the moment is related to the second derivative of
while the shear is related to the third derivative of
. It takes up to a polynomial of the seventh degree for
for the approximate shearing force to have a flat distribution similar to the exact solution.
View Mathematica Code
Clear[y, b, L, t, Ii, Ee, P]
b = 1/4;
L = 8;
t = (1/2 - 1/4*X1/L);
Ii = b*t^3/12;
Ee = 20*10^9;
P = 10000;
a = DSolve[{D[Ee*Ii*y''[X1], {X1, 2}] == 0, y[0] == 0, y'[0] == 0, y''[L] == 0, (D[Ee*Ii, X1] /. X1 -> L)*y''[L] + (Ee*Ii /. X1 -> L)*y'''[L] == P}, y[X1], X1, Assumptions -> 0 <= X1 <= L]
y = FullSimplify[Chop[N[y[X1] /. a[[1]]]]]
th = D[y, X1];
M = FullSimplify[Ee*Ii*D[y, {X1, 2}]]
V = FullSimplify[D[M, X1]]
Plot[y, {X1, 0, L}, PlotLabel -> "Displacement y", AxesLabel -> {"X1 (m)", "y (m)"}, PlotLegends -> {"y exact"}]
Plot[M, {X1, 0, L}, PlotLabel -> "Moment M", AxesLabel -> {"X1 (m)", "M (N.m)"}, PlotLegends -> {"M exact"}]
Plot[V, {X1, 0, L}, PlotRange -> {{0, L}, {0, 2 P}}, PlotLabel -> "Shear V", AxesLabel -> {"X1 (m)", "V (N)"}, PlotLegends -> {"V exact"}]
Print["Second Order Polynomial"]
y2 = a2*X1^2;
d2y = D[y2, {X1, 2}];
PE = Integrate[Ee*Ii/2*d2y^2, {X1, 0, L}] - (-P)*y2 /. X1 -> L;
sol = Solve[D[PE, a2] == 0, a2]
y2 = y2 /. sol[[1]]
M2 = FullSimplify[Ee*Ii*D[y2, {X1, 2}]]
V2 = FullSimplify[D[M2, X1]]
Plot[{y, y2}, {X1, 0, L}, PlotLabel -> "Displacement y", AxesLabel -> {"X1 (m)", "y (m)"}, PlotLegends -> {"y exact", "y2"}]
Plot[{M, M2}, {X1, 0, L}, PlotLabel -> "Moment M", AxesLabel -> {"X1 (m)", "M (N.m)"}, PlotLegends -> {"M exact", "M2"}]
Plot[{V, V2}, {X1, 0, L}, PlotRange -> {{0, L}, {0, 2 P}}, PlotLabel -> "Shear V", AxesLabel -> {"X1 (m)", "V (N)"}, PlotLegends -> {"V exact", "V2"}]
Print["Third Order Polynomial"]
y3 = a2*X1^2 + a3*X1^3;
d2y = D[y3, {X1, 2}];
PE = Integrate[Ee*Ii/2*d2y^2, {X1, 0, L}] - (-P)*y3 /. X1 -> L;
sol = Solve[{D[PE, a2] == 0, D[PE, a3] == 0}, {a2, a3}]
y3 = y3 /. sol[[1]]
M3 = FullSimplify[Ee*Ii*D[y3, {X1, 2}]]
V3 = FullSimplify[D[M3, X1]]
Plot[{y, y3}, {X1, 0, L}, PlotLabel -> "Displacement y", AxesLabel -> {"X1 (m)", "y (m)"}, PlotLegends -> {"y exact", "y3"}]
Plot[{M, M3}, {X1, 0, L}, PlotLabel -> "Moment M", AxesLabel -> {"X1 (m)", "M (N.m)"}, PlotLegends -> {"M exact", "M3"}]
Plot[{V, V3}, {X1, 0, L}, PlotRange -> {{0, L}, {0, 2 P}}, PlotLabel -> "Shear V", AxesLabel -> {"X1 (m)", "V (N)"}, PlotLegends -> {"V exact", "V3"}]
Print["Fourth Order Polynomial"]
y4 = a2*X1^2 + a3*X1^3 + a4*X1^4;
d2y = D[y4, {X1, 2}];
PE = Integrate[Ee*Ii/2*d2y^2, {X1, 0, L}] - (-P)*y4 /. X1 -> L;
sol = Solve[{D[PE, a2] == 0, D[PE, a3] == 0, D[PE, a4] == 0}, {a2, a3, a4}]
y4 = y4 /. sol[[1]]
M4 = FullSimplify[Ee*Ii*D[y4, {X1, 2}]]
V4 = FullSimplify[D[M4, X1]]
Plot[{y, y4}, {X1, 0, L}, PlotLabel -> "Displacement y", AxesLabel -> {"X1 (m)", "y (m)"}, PlotLegends -> {"y exact", "y4"}]
Plot[{M, M4}, {X1, 0, L}, PlotLabel -> "Moment M", AxesLabel -> {"X1 (m)", "M (N.m)"}, PlotLegends -> {"M exact", "M4"}]
Plot[{V, V4}, {X1, 0, L}, PlotRange -> {{0, L}, {0, 2 P}}, PlotLabel -> "Shear V", AxesLabel -> {"X1 (m)", "V (N)"}, PlotLegends -> {"V exact", "V4"}]
Print["Fifth Order Polynomial"]
y5 = a2*X1^2 + a3*X1^3 + a4*X1^4 + a5*X1^5;
d2y = D[y5, {X1, 2}];
PE = Integrate[Ee*Ii/2*d2y^2, {X1, 0, L}] - (-P)*y5 /. X1 -> L;
sol = Solve[{D[PE, a2] == 0, D[PE, a3] == 0, D[PE, a4] == 0, D[PE, a5] == 0}, {a2, a3, a4, a5}]
y5 = y5 /. sol[[1]]
M5 = FullSimplify[Ee*Ii*D[y5, {X1, 2}]]
V5 = FullSimplify[D[M5, X1]]
Plot[{y, y5}, {X1, 0, L}, PlotLabel -> "Displacement y", AxesLabel -> {"X1 (m)", "y (m)"}, PlotLegends -> {"y exact", "y5"}]
Plot[{M, M5}, {X1, 0, L}, PlotLabel -> "Moment M", AxesLabel -> {"X1 (m)", "M (N.m)"}, PlotLegends -> {"M exact", "M5"}]
Plot[{V, V5}, {X1, 0, L}, PlotRange -> {{0, L}, {0, 2 P}}, PlotLabel -> "Shear V", AxesLabel -> {"X1 (m)", "V (N)"}, PlotLegends -> {"V exact", "V5"}]
Print["Sixth Order Polynomial"]
y6 = a2*X1^2 + a3*X1^3 + a4*X1^4 + a5*X1^5 + a6*X1^6;
d2y = D[y6, {X1, 2}];
PE = Integrate[Ee*Ii/2*d2y^2, {X1, 0, L}] - (-P)*y6 /. X1 -> L;
sol = Solve[{D[PE, a2] == 0, D[PE, a3] == 0, D[PE, a4] == 0, D[PE, a5] == 0, D[PE, a6] == 0}, {a2, a3, a4, a5, a6}]
y6 = y6 /. sol[[1]]
M6 = FullSimplify[Ee*Ii*D[y6, {X1, 2}]]
V6 = FullSimplify[D[M6, X1]]
Plot[{y, y6}, {X1, 0, L}, PlotLabel -> "Displacement y", AxesLabel -> {"X1 (m)", "y (m)"}, PlotLegends -> {"y exact", "y6"}]
Plot[{M, M6}, {X1, 0, L}, PlotLabel -> "Moment M", AxesLabel -> {"X1 (m)", "M (N.m)"}, PlotLegends -> {"M exact", "M6"}]
Plot[{V, V6}, {X1, 0, L}, PlotRange -> {{0, L}, {0, 2 P}}, PlotLabel -> "Shear V", AxesLabel -> {"X1 (m)", "V (N)"}, PlotLegends -> {"V exact", "V6"}]
Print["Seventh Order Polynomial"]
y7 = a2*X1^2 + a3*X1^3 + a4*X1^4 + a5*X1^5 + a6*X1^6 + a7*X1^7;
d2y = D[y7, {X1, 2}];
PE = Integrate[Ee*Ii/2*d2y^2, {X1, 0, L}] - (-P)*y7 /. X1 -> L;
sol = Solve[{D[PE, a2] == 0, D[PE, a3] == 0, D[PE, a4] == 0, D[PE, a5] == 0, D[PE, a6] == 0, D[PE, a7] == 0}, {a2, a3, a4, a5,a6, a7}]
y7 = y7 /. sol[[1]]
M7 = FullSimplify[Ee*Ii*D[y7, {X1, 2}]]
V7 = FullSimplify[D[M7, X1]]
Plot[{y, y7}, {X1, 0, L}, PlotLabel -> "Displacement y", AxesLabel -> {"X1 (m)", "y (m)"}, PlotLegends -> {"y exact", "y7"}]
Plot[{M, M7}, {X1, 0, L}, PlotLabel -> "Moment M", AxesLabel -> {"X1 (m)", "M (N.m)"}, PlotLegends -> {"M exact", "M7"}]
Plot[{V, V7}, {X1, 0, L}, PlotRange -> {{0, L}, {0, 2 P}}, PlotLabel -> "Shear V", AxesLabel -> {"X1 (m)", "V (N)"}, PlotLegends -> {"V exact", "V7"}]
Plot[{y, y2, y3, y4, y5, y6, y7}, {X1, 0, L}, PlotLabel -> "Displacement y", AxesLabel -> {"X1 (m)", "y (m)"}, PlotLegends -> {"y exact", "y2", "y3", "y4", "y5", "y6", "y7"}]
Plot[{M, M2, M3, M4, M5, M6, M7}, {X1, 0, L}, PlotLabel -> "Moment M", AxesLabel -> {"X1 (m)", "M (N.m)"}, PlotLegends -> {"M exact", "M2", "M3", "M4", "M5", "M6", "M7"}]
Plot[{V,V2,V3,V4,V5,V6,V7}, {X1, 0, L}, PlotRange ->{{0, L}, {0, 2 P}},PlotLabel->"Shear V",AxesLabel -> {"X1 (m)", "V (N)"},PlotLegends -> {"V exact", "V2", "V3", "V4", "V5", "V6", "V7"}]
View Python Code
from sympy import *
from numpy import *
import sympy as sp
import numpy as np
import matplotlib.pyplot as plt
sp.init_printing(use_latex = "mathjax")
x, a2, a3, a4, a5, a6, a7, EI, L, E= symbols("x a_2 a_3 a_4 a_5 a_6 a_7 EI L E")
display("Exact")
b = 1/4
L = 8
t = (1/2-1/4*x/L)
I = b*t**3/12
E = 20*10**9
P = 10000
# Exact solution done in mathematica code
y = (34.8872 + (-2.96688 + 0.036864*x)*x + (-12.5829 + 0.786432*x)*sp.ln(16-x))/(x-16)
M = -80000 + 10000*x
V = M.diff(x)
display("y,M,V: ", y, M, V)
display("Second Order Polynomial")
y2 = a2*x**2
d2y = y2.diff(x,2)
PE2 = integrate(E*I/2*d2y**2, (x,0,L)) - (-P)*y2.subs(x,L)
display("Potential Energy: ", PE2)
sol = solve(PE2.diff(a2), a2)
y2 = y2.subs(a2,sol[0])
M2 = E*I*y2.diff(x,2)
V2 = M2.diff(x)
display("y2,M2,V2: ", y2, M2, V2)
display("Third Order Polynomial")
y3 = a2*x**2+a3*x**3
d2y = y3.diff(x,2)
PE3 = integrate(E*I/2*d2y**2, (x,0,L)) - (-P)*y3.subs(x,L)
display("Potential Energy: ", PE3)
sol = solve((PE3.diff(a2), PE3.diff(a3)), (a2,a3))
y3 = y3.subs({a2:sol[a2], a3:sol[a3]})
M3 = E*I*y3.diff(x,2)
V3 = M3.diff(x)
display("y3,M3,V3: ", y3, M3, V3)
display("Fourth Order Polynomial")
y4 = a2*x**2+a3*x**3+a4*x**4
d2y = y4.diff(x,2)
PE4 = integrate(E*I/2*d2y**2, (x,0,L)) - (-P)*y4.subs(x,L)
display("Potential Energy: ", PE4)
sol = solve((PE4.diff(a2), PE4.diff(a3), PE4.diff(a4)), (a2,a3,a4))
y4 = y4.subs({a2:sol[a2], a3:sol[a3], a4:sol[a4]})
M4 = E*I*y4.diff(x,2)
V4 = M4.diff(x)
display("y4,M4,V4: ", y4, M4, V4)
display("Fifth Order Polynomial")
y5 = a2*x**2+a3*x**3+a4*x**4+a5*x**5
d2y = y5.diff(x,2)
PE5 = integrate(E*I/2*d2y**2, (x,0,L)) - (-P)*y5.subs(x,L)
display("Potential Energy: ", PE5)
sol = solve((PE5.diff(a2), PE5.diff(a3), PE5.diff(a4), PE5.diff(a5)), (a2,a3,a4,a5))
y5 = y5.subs({a2:sol[a2], a3:sol[a3], a4:sol[a4], a5:sol[a5]})
M5 = E*I*y5.diff(x,2)
V5 = M5.diff(x)
display("y5,M5,V5: ", y5, M5, V5)
display("Sixth Order Polynomial")
y6 = a2*x**2+a3*x**3+a4*x**4+a5*x**5+a6*x**6
d2y = y6.diff(x,2)
PE6 = integrate(E*I/2*d2y**2, (x,0,L)) - (-P)*y6.subs(x,L)
display("Potential Energy: ", PE6)
sol = solve((PE6.diff(a2), PE6.diff(a3), PE6.diff(a4), PE6.diff(a5), PE6.diff(a6)), (a2,a3,a4,a5,a6))
y6 = y6.subs({a2:sol[a2], a3:sol[a3], a4:sol[a4], a5:sol[a5], a6:sol[a6]})
M6 = E*I*y6.diff(x,2)
V6 = M6.diff(x)
display("y6,M6,V6: ", y6, M6, V6)
display("Seventh Order Polynomial")
y7 = a2*x**2+a3*x**3+a4*x**4+a5*x**5+a6*x**6+a7*x**7
d2y = y7.diff(x,2)
PE7 = integrate(E*I/2*d2y**2, (x,0,L)) - (-P)*y7.subs(x,L)
display("Potential Energy: ", PE7)
sol = solve((PE7.diff(a2), PE7.diff(a3), PE7.diff(a4), PE7.diff(a5), PE7.diff(a6), PE7.diff(a7)), (a2,a3,a4,a5,a6,a7))
y7 = y7.subs({a2:sol[a2], a3:sol[a3], a4:sol[a4], a5:sol[a5], a6:sol[a6], a7:sol[a7]})
M7 = E*I*y7.diff(x,2)
V7 = M7.diff(x)
display("y7,M7,V7: ", y7, M7, V7)
#Graph
fig, ax = plt.subplots(3,1, figsize = (10,15))
plt.setp(ax[0], xlabel = "x(m) ", ylabel = "y(m)")
plt.setp(ax[1], xlabel = "x(m)", ylabel = "M(N*m)")
plt.setp(ax[2], xlabel = "x(m)", ylabel = "V(N)")
x1 = np.arange(0,8.1,0.1)
y_list = [y.subs({x:i}) for i in x1]
y2_list = [y2.subs({x:i}) for i in x1]
y3_list = [y3.subs({x:i}) for i in x1]
y4_list = [y4.subs({x:i}) for i in x1]
y5_list = [y5.subs({x:i}) for i in x1]
y6_list = [y6.subs({x:i}) for i in x1]
y7_list = [y7.subs({x:i}) for i in x1]
M_list = [M.subs({x:i}) for i in x1]
M2_list = [M2.subs({x:i}) for i in x1]
M3_list = [M3.subs({x:i}) for i in x1]
M4_list = [M4.subs({x:i}) for i in x1]
M5_list = [M5.subs({x:i}) for i in x1]
M6_list = [M6.subs({x:i}) for i in x1]
M7_list = [M7.subs({x:i}) for i in x1]
V_list = [V.subs({x:i}) for i in x1]
V2_list = [V2.subs({x:i}) for i in x1]
V3_list = [V3.subs({x:i}) for i in x1]
V4_list = [V4.subs({x:i}) for i in x1]
V5_list = [V5.subs({x:i}) for i in x1]
V6_list = [V6.subs({x:i}) for i in x1]
V7_list = [V7.subs({x:i}) for i in x1]
ax[0].plot(x1, y_list, 'blue', label = "Exact y")
ax[0].plot(x1, y2_list, 'orange', label = "y2")
ax[0].plot(x1, y3_list, 'green', label = "y3")
ax[0].plot(x1, y4_list, 'red', label = "y4")
ax[0].plot(x1, y5_list, 'purple', label = "y5")
ax[0].plot(x1, y6_list, 'brown', label = "y6")
ax[0].plot(x1, y7_list, 'cyan', label = "y7")
ax[0].legend()
ax[1].plot(x1, M_list, 'blue', label = "Exact M")
ax[1].plot(x1, M2_list, 'orange', label = "y2")
ax[1].plot(x1, M3_list, 'green', label = "y3")
ax[1].plot(x1, M4_list, 'red', label = "y4")
ax[1].plot(x1, M5_list, 'purple', label = "y5")
ax[1].plot(x1, M6_list, 'brown', label = "y6")
ax[1].plot(x1, M7_list, 'cyan', label = "y7")
ax[1].legend()
ax[2].plot(x1, V_list, 'blue', label = "Exact V")
ax[2].plot(x1, V2_list, 'orange', label = "y2")
ax[2].plot(x1, V3_list, 'green', label = "y3")
ax[2].plot(x1, V4_list, 'red', label = "y4")
ax[2].plot(x1, V5_list, 'purple', label = "y5")
ax[2].plot(x1, V6_list, 'brown', label = "y6")
ax[2].plot(x1, V7_list, 'cyan', label = "y7")
ax[2].legend()