# Problem Derivative (D) of function defined with Piecewise or UnitStep functions - SIMPLE BEAM PROBLEM, 3 SUPPORTS

I am solving a simple structural beam problem using slope-deflection approach - 3 support A(left- fixed), B(middle- pin), C(right - pin). Beam AB - uniform load, Beam BC - concentrated load in middle.

I was successful in getting the complete full length beam moment equations using a piecewise and unitstep function. I was able to get the moment at specific points along the beam length and get plots comparing the 2 approaches.

Now when I need to take the 1 st derive to get a shear plot or evaluate the shear at a specific point, the derived function does not work - no matter if I take the derivative of either moment function using the Piecewise or UnitStep functions. ?? For some reason, the derived function is not considering the x distance point as a variable ? My ultimate goal is to get shear, deflection, slope plots and evaluate the values as specific points.

(*3 supports A(left end- fixed), B(middle-pinned), C(right end-pinned), uniform and concentrated load *)



qc = 1;
qu = 1;
lab = 10;
lbc = 10;
deltaab = 0;
deltabc = 0;
anglea = 0;
ei = 1;
momableft = (2.*ei)/lab*(2*anglea + angleb - (3*deltaab)/lab) + (-qu*
lab^2)/12;
mombaright = (2.*ei)/lab*(anglea + 2*angleb - (3*deltaab)/lab) + (+qu*
lab^2)/12;
mombcleft = (2.*ei)/lbc*(2*angleb + anglec - (3*deltabc)/lbc) + (-qc*
lbc)/8;
momcbright = (2.*ei)/lbc*(angleb + 2*anglec - (3*deltabc)/lbc) + (+qc*
lbc)/8;

ssol1 = Solve[
(mombaright + mombcleft) == 0  &&
(momcbright) == 0 ,
{angleb, angleb, anglec} ];

anglebValue = angleb /. ssol1[[1]];
anglecValue = anglec /. ssol1[[1]];

angleaValue = anglea;

Print[" The angle for support A: (ei=1) ", angleaValue]
Print[" The angle for support B: (ei=1) ", anglebValue]
Print[" The angle for support C: (ei=1) ", anglecValue]
momableft = (2.*ei)/
lab*(2*angleaValue + anglebValue - (3*deltaab)/lab) + (-qu*lab^2)/
12;
mombaright = (2.*ei)/
lab*(angleaValue + 2*anglebValue - (3*deltaab)/lab) + (+qu*lab^2)/
12;
mombcleft = (2.*ei)/
lab*(2*anglebValue + anglecValue - (3*deltabc)/lbc) + (-qc*lab)/
8;
momcbright = (2.*ei)/
lab*(anglebValue + 2*anglecValue - (3*deltabc)/lbc) + (+qc*lab)/8;

Print[" ********* Moments beam A-B & B-C *********  "];
Print[" Moment AB: ", momableft, " ft-kips"];
Print[" Moment BA: ", mombaright, " ft-kips"];
Print[" Moment BC: ", mombcleft, " ft-kips"];
Print[" Moment CB: ", momcbright, " ft-kips"];

********* Moments beam A-B & B-C *********

Moment AB: -10.1786 ft-kips

Moment BA: 4.64286 ft-kips

Moment BC: -4.64286 ft-kips

Moment CB: 0. ft-kips

(*  now for each beam in the model and using the end moments just \

solutionA =
Solve[(momableft + mombaright - (qu*lab^2)/2 +
rAbeamabqusolve*lab) == 0,
{rAbeamabqusolve}]; (*moment about point B *)
rAbeamabqusolveValue = rAbeamabqusolve /. solutionA[[1]];

solutionB =
Solve[(momableft + mombaright + (qu*lab^2)/2 -
rBbeamabqusolve*lab) == 0,
{rBbeamabqusolve}]; (*moment about point A *)
rBbeamabqusolveValue = rBbeamabqusolve /. solutionB[[1]];

Print["   "];
Print[" ********* Reaction beam A-B *********  "];
Print[" Beam A-B, REACTION A: ", rAbeamabqusolveValue, " kips"];
Print[" Beam A-B, REACTION B: ", rBbeamabqusolveValue, " kips"];
Print[" ZERO FORCE BALANCE CHECK ",
rAbeamabqusolveValue + rBbeamabqusolveValue - qu*lab, " kips"];

solutionBbc =
Solve[(mombcleft + momcbright - (qc*lbc)/2 + rBbeambcqcsolve*lbc) ==
0,
{rBbeambcqcsolve}]; (*moment about point C *)
rBbeambcqcsolveValue = rBbeambcqcsolve /. solutionBbc[[1]];

solutionCbc =
Solve[(mombcleft + momcbright + (qc*lbc)/2 - rCbeambcqcsolve*lbc) ==
0,
{rCbeambcqcsolve}]; (*moment about point B *)
rCbeambcqcsolveValue = rCbeambcqcsolve /. solutionCbc[[1]];

Print["   "];
Print[" ********* Reaction beam B-C *********  "];
Print[" Beam B-C, REACTION B: ", rBbeambcqcsolveValue, " kips"];
Print[" Beam B-C, REACTION C: ", rCbeambcqcsolveValue, " kips"];
Print[" ZERO FORCE BALANCE CHECK ",
rBbeambcqcsolveValue + rCbeambcqcsolveValue - qc , " kips"];

********* Reaction beam A-B *********

Beam A-B, REACTION A: 5.55357 kips

Beam A-B, REACTION B: 4.44643 kips

ZERO FORCE BALANCE CHECK 0. kips

********* Reaction beam B-C *********

Beam B-C, REACTION B: 0.964286 kips

Beam B-C, REACTION C: 0.0357143 kips

ZERO FORCE BALANCE CHECK 0. kips

(*     TOTAL MODEL REACTIONS   *)

Ra = rAbeamabqusolveValue;
Rb = rBbeamabqusolveValue + rBbeambcqcsolveValue;
Rc = rCbeambcqcsolveValue;

Print[" TOTAL MODEL REACTIONS A: ", Ra, " kips"]
Print[" TOTAL MODEL REACTIONS B: ", Rb, " kips"]
Print[" TOTAL MODEL REACTIONS C: ", Rc, " kips"]
Print[" TOTAL MODEL REACTION check: ",
Ra + Rb + Rc - qc - qu*lab, " kips"]

TOTAL MODEL REACTIONS A: 5.55357 kips

TOTAL MODEL REACTIONS B: 5.41071 kips

TOTAL MODEL REACTIONS C: 0.0357143 kips

TOTAL MODEL REACTION check: 1.77636*10^-15 kips

(*Define the total beam  model functions  *)

momABCtotal[x_] := Piecewise[{
{momableft - (qu*x^2)/2 + (Ra*x), x < lab},
{momableft + mombaright +
mombcleft - (qu*lab)*(x - lab/2) + (Ra*x) + Rb*(x - lab),
x == lab },
{momableft + mombaright +
mombcleft - (qu*lab)*(x - lab/2) + (Ra*x) + Rb*(x - lab),
x > lab && x < lab + (lbc/2)},
{momableft + mombaright + mombcleft +
momcbright - (qu*lab)*(x - lab/2) + (Ra*x) + Rb*(x - lab) +
Rc*(x), x == lbc + lab}
}]

momABCtotalUSTEP3[x_] :=
(momableft - (qu*x^2)/2 + (Ra*x))*UnitStep[lab - x] +
(momableft - (qu*lab)*(x - lab/2) + (Ra*x) + Rb*(x - lab))*
UnitStep[x - lab]*UnitStep[(lab + lbc/2) - x] +
(momableft - (qu*lab)*(x - lab/2) - qc*(x - (lab + lbc/2)) + Ra*x +
Rb*(x - lab))*UnitStep[x - (lab + (lbc/2))]  +
UnitStep[(lab + lbc) - Infinity]

momABCtotal[3] (*     see if function can evaluate value at x = 3 \
point *)
momABCtotalUSTEP3[3]

Plot[momABCtotal[x], {x, 0, lab + lbc}, PlotRange -> All,
AxesLabel -> {x, Bending_Moment (k - ft)[x]}]

Plot[momABCtotalUSTEP3[x], {x, 0, lab + lbc}, PlotRange -> All,
AxesLabel -> {x, USTEP3_Bending _Moment (k - ft)[x]}]

Plot[{momABCtotal[x], momABCtotalUSTEP3[x]}, {x, 0, lab + lbc},
PlotLegends -> {"Piecewise", "UnitStep"}, PlotRange -> All]

shearmomABCtotal[x_][x_] := D[momABCtotal[x], x]
shearABCtotalUSTEP3[x_] := D[momABCtotalUSTEP3[x], x]

shearmomABCtotal[3] (*     see if function can evaluate value at x = \
3 point *)
shearABCtotalUSTEP3[3]

simplifiedshearmomABCtotal[x_] := Simplify[shearmomABCtotal[x]]
simplifiedShearABCtotalUSTEP3[x_] := Simplify[shearABCtotalUSTEP3[x]]

simplifiedshearmomABCtotal[3]
simplifiedShearABCtotalUSTEP3[3]

(*Evaluate the shear force function at a specific point*)
shearValueAt1 = simplifiedShearABCtotalUSTEP3[1]
shearValueAt9 = simplifiedShearABCtotalUSTEP3[9]

(*Output the value*)
Print["shearValueAt1 = ", shearValueAt1]
Print["shearValueAt9 = ", shearValueAt9]

(*Display the shear force function*)
simplifiedShearABCtotalUSTEP3[x]

(*Evaluate the shear force function at specific points*)
shearValueAt9 = simplifiedShearABCtotalUSTEP3[9]
shearValueAt4 = simplifiedShearABCtotalUSTEP3[4]
(*Plot the shear force function*)

Plot[simplifiedShearABCtotalUSTEP3[x], {x, 0, 20},
PlotLabel -> "Shear Force Function", AxesLabel -> {"x", "V(x)"}]

Plot[simplifiedShearABCtotalUSTEP3[x_], {x, 0, 20},
PlotLabel -> "Shear Force Function", AxesLabel -> {"x", "V(x)"}]

Plot[shearForce[x], {x, 0, lab + lbc}, PlotRange -> All,
AxesLabel -> {x, Shear (k)[x]}]
Plot[shearABCtotalUSTEP3[x], {x, 0, lab + lbc}, PlotRange -> All,
AxesLabel -> {x, Shear (k)[x]}]
Plot[simplifiedShearABCtotalUSTEP3[x], {x, 0, lab + lbc},
PlotRange -> All, AxesLabel -> {x, Shear (k)[x]}]



moment plots are correct and match hand calcs.

• This is not a minimum working example of your problem. You should pay attention to the error messages. You are not defining your derivatives correctly. You are using the form g[x_] := D[f[x], x]; however, when called x has a value and things like D[f[3], 3] make no sense. Instead use the form g[x_] := Derivative[1][f][x] so that the derivative is taken before the value is assigned to the variable. Alternatively, use Set instead of SetDelayed, e.g., g[x_] = D[f[x], x] Commented Jun 16 at 18:45