Problem 3.2

"Length of the first beam segment";
a=1;
"Length of the second beam segment";
b=0.8;
"Diameter of the first cross section";
d1=0.06;
"Area moment of inertia of the first cross section";
Iz1=(d1^4) π/64;
"Diameter of the second cross section";
d2=0.046;
"Area moment of inertia of the second cross section";
Iz2=(d2^4) π/64;
"Young modulus of the first beam segment";
E1=210 (10^9);
"Young modulus of the second beam segment";
E2=190  (10^9);
"Stiffness of the spring";
s=2 (10^5);
"Value of the distributed force";
p=800;
"Value of the concentrated force";
F=1000;
"Calculation of the exact solution";
vs=DSolve[{{v1''[x]==((Bf(a+b-x))/(Iz2 E2)),v2''[x]==((Bf(a+b-x))/(Iz2 E2))-((F (a+b/2-x))/(Iz2 E2)),v3''[x]==((Bf(a+b-x))/(Iz1 E1))-((F (a+b/2-x))/(Iz1 E1))-((p (((a-x))^2))/(2 Iz1 E1))},{v3[0]==0,v3'[0]==0,v2[a]==v3[a],v2'[a]==v3'[a],v1[a+b/2]==v2[a+b/2],v1'[a+b/2]==v2'[a+b/2]}},{v1,v2,v3},x][[1]];
BF=Solve[-s (v1[a+b]/.vs)==Bf,Bf][[1,1]];
{vE1,vE2,vE3}={v1,v2,v3}/.vs/.BF;
"Plotting the exact results";
PvE=Show[Plot[vE1[x],{x,a+b/2,a+b},PlotStyle->Red,PlotRange->{{0,a+b},{0.0002,-0.003}}],Plot[vE2[x],{x,a,a+b/2},PlotStyle->Red,PlotRange->{{0,a+b},{0.0002,-0.003}}],Plot[vE3[x],{x,0,a},PlotStyle->Red,PlotRange->{{0,a+b},{0.0002,-0.003}}]];
{φE1,φE2,φE3}=D[{vE1[x],vE2[x],vE3[x]},x];
PφE=Show[Plot[φE1,{x,a+b/2,a+b},PlotStyle->Red,PlotRange->{{0,a+b},{0.0002,-0.003}}],Plot[φE2,{x,a,a+b/2},PlotStyle->Red,PlotRange->{{0,a+b},{0.0002,-0.003}}],Plot[φE3,{x,0,a},PlotStyle->Red,PlotRange->{{0,a+b},{0.0002,-0.003}}]];
{Mh1,Mh2,Mh3}={-Iz2 E2,-Iz2 E2,-Iz1 E1} *D[{vE1[x],vE2[x],vE3[x]},{x,2}];
PMhE=Show[Plot[Mh1,{x,a+b/2,a+b},PlotStyle->Red,PlotRange->{{0,a+b},{-250,900}}],Plot[Mh2,{x,a,a+b/2},PlotStyle->Red,PlotRange->{{0,a+b},{-250,900}}],Plot[Mh3,{x,0,a},PlotStyle->Red,PlotRange->{{0,a+b},{-250,900}}]];
{V1,V2,V3}=-D[{Mh1,Mh2,Mh3},x];
PVE=Show[Plot[V1,{x,a+b/2,a+b},PlotStyle->Red,PlotRange->{{0,a+b},{-600,1300}}],Plot[V2,{x,a,a+b/2},PlotStyle->Red,PlotRange->{{0,a+b},{-600,1300}}],Plot[V3,{x,0,a},PlotStyle->Red,PlotRange->{{0,a+b},{-600,1300}}]];
"Calculation of the approximate solution using the method of Ritz";
"Shape function";
ω=(x^2);
"Lowest order of approximation";
Nl=0;
"Highest order of approximation";
Nh=8;
Nappr=Nl;
vp=ConstantArray[0,{Nh-Nl+1}];
φp=vp;
Mp1=vp;
Mp2=vp;
Vp1=vp;
Vp2=vp;
vB=vp;
φB=vp;
Do[
B=ω Table[{(x^(i-1))},{i,1,Nappr+1}];
Bt=Transpose[B];
dB=D[B,x];
ddB=D[B,{x,2}];
ddBt=Transpose[ddB];
S=Iz1 E1 Integrate[ddB.ddBt,{x,0,a}]+Iz2 E2 Integrate[ddB.ddBt,{x,a,a+b}]+s ((B.Bt)/.x->a+b);
Q=-p Integrate[B,{x,0,a}]-F (B/.x->a+b/2);
A=Inverse[S].Q;
v=Transpose[B].A;
φ=D[v,x];
Mh1=-Iz1 E1 D[v,{x,2}];
Mh2=-Iz2 E2 D[v,{x,2}];
V1=-D[Mh1,x];
V2=-D[Mh2,x];
vp[[j]]=v;
φp[[j]]=φ;
Mp1[[j]]=Mh1;
Mp2[[j]]=Mh2;
Vp1[[j]]=V1;
Vp2[[j]]=V2;
vB[[j]]=(v/.x->a+b)[[1,1]];
φB[[j]]=(φ/.x->a+b)[[1,1]];
Nappr=Nappr+1;
,{j,1,Nh-Nl+1}]
Do[
vB[[i]]={i+Nl-1,vB[[i]]};
φB[[i]]={i+Nl-1,φB[[i]]};
,{i,1,Nh-Nl+1}]
"Animation of the results";
Animate[Grid[{{Style[StringJoin[{"N=",ToString[k-1+Nl]}],20],SpanFromLeft},{Row[{Style["Exact ",16],Graphics[{Red,Rectangle[{0,0},{1,0.03}]}],Style["  Ritz's method ",16],Graphics[{Blue,Rectangle[{0,0},{1,0.03}]}]}],SpanFromLeft},{Show[PvE,Plot[vp[[k]],{x,0,a+b}],AxesStyle->10,PlotRange->{0.0002,-0.003},AxesLabel->{Style["x [m]",14],Style["v(x) [m]",14]},ImageSize->{400,200}],Show[PφE,Plot[φp[[k]],{x,0,a+b}],AxesStyle->10,PlotRange->{0.0002,-0.003},AxesOrigin->{0,0},AxesLabel->{Style["x [m]",14],Style["φ(x) [rad]",14]},ImageSize->{400,200}]},{Show[{PMhE,Plot[Mp1[[k]],{x,0,a}],Plot[Mp2[[k]],{x,a,a+b}]},AxesStyle->10,PlotRange->{-250,900},AxesLabel->{Style["x [m]",14],Style["Mh(x) [Nm]",14]},ImageSize->{400,200}],ListPlot[vB[[1;;k]],AxesStyle->10,PlotRange->{{0,Nh+1},{0.0002,-0.003}},AxesLabel->{Style["N",14],Style["(Subscript[v, B]) [m]",14]},PlotMarkers->{Automatic,Medium},ImageSize->{400,200}]},{Show[{PVE,Plot[Vp1[[k]],{x,0,a}],Plot[Vp2[[k]],{x,a,a+b}]},AxesStyle->10,PlotRange->{-800,1500},AxesLabel->{Style["x [m]",14],Style["V(x) [N]",14]},ImageSize->{400,200}],ListPlot[φB[[1;;k]],AxesStyle->10,PlotRange->{{0,Nh+1},{0.0002,-0.003}},AxesLabel->{Style["N",14],Style["(Subscript[φ, B]) [rad]",14]},PlotMarkers->{Automatic,Medium},ImageSize->{400,200}]}},Dividers->{{},{2->Thick,3->Thick}},Spacings->{{},{2,2,2}}],{k,1,Nh-Nl+1,1},AnimationRate->0.5]