(*5章4節 図5．4の計算*)

(*Rohde の論文の計算*)

(*Rohde,F.V.,Large deflections of a cantilever beam with a uniformly
disributed load,Quart.,Appl.Math., 11(1953),p.337.*)

(*この結果をもとに変形図(Fig5-4)を描く*)

Clear["Global`*"]

(*級数項の上限*)

nn=12

(*\theta の級数和の表現式*)

tht=0;
Do[tht=tht+a[i]*s^i,{i,0,nn}]

tht

T=tht/.{a[0]->0,a[1]->0}

(*Governing equation*)

eq1=D[tht,{s,2}]+w*s/EI*(Cos[alp]*Cos[T]-Sin[alp]*Sin[T])

eq2=FullSimplify[Series[eq1,{s,0,nn}]]

eqs0=Coefficient[eq2,s,0]

eqs1=Coefficient[eq2,s,1]

eqs2=Coefficient[eq2,s,2]

eqs3=Coefficient[eq2,s,3]

eqs4=Coefficient[eq2,s,4]

eqs5=Coefficient[eq2,s,5]

eqs6=Coefficient[eq2,s,6]

Do[
eqs[i]=Coefficient[eq2,s,i];
Print[i," , ",eqs[i]],{i,0,nn-2}]

(*連立方程式を解くためのリスト作成*)

eqlst={};anlst={};
Do[eqlst=Append[eqlst,eqs[i]==0];
anlst=Append[anlst,a[i+2]];
,{i,0,nn-2}]

eqlst

anlst

ans1=Solve[eqlst,anlst][[1]]

tht1=(tht/.ans1)/.{a[0]->0,a[1]->0}

Series[Cos[tht1],{s,0,nn-1}]

(*級数展開*)

T1=FullSimplify[Integrate[Series[Cos[tht1],{s,0,nn}],{s,0,s}]]

Do[Print[i," , ",Coefficient[T1,s,i]],{i,1,nn}]

T2=FullSimplify[Integrate[Series[Sin[tht1],{s,0,nn}],{s,0,s}]]

Coefficient[T2,s,nn-2]

Do[Print[i," , ",Coefficient[T2,s,i]],{i,1,nn}]

defy1=T1*Sin[alp]+T2*Cos[alp]

s*Coefficient[defy1,s,1]+Coefficient[defy1,s,4]*s^4

defx1=T1*Cos[alp]-T2*Sin[alp]

(*無次元化*)

defy2=FullSimplify[((defy1)/.{w->EI*w1/L^3,s->L*s1})/L]

defx2=FullSimplify[((defx1)/.{w->EI*w1/L^3,s->L*s1})/L]

(* Deflections of the beam *)

(* Rohde's Paper equation (3)*)

alp1=0;Do[alp1=alp1-a[i]*L^i,{i,3,nn}]

alp1

alp2=FullSimplify[alp1/.ans1]

(*Preparation of relation between alpha and w *)

alp3=alp2-alp

nonline1=FullSimplify[(alp3)/.w->w1*EI/L^3]

Plot[nonline1/.alp->20*Pi/180,{w1,-2,5}]

FindRoot[nonline1/.alp->20*Pi/180,{w1,0}]

(*alpha を与えてw1（無次元荷重）を求める*)

Print["alpha     x/L        y/L "];
deflst[1]={};deflst[2]={};deflst[3]={};deflst[4]={};deflst[5]={};deflst[6]={};deflst[7]={};
Do[
  alp1=i*10*Pi/180;
  answ1=FindRoot[(nonline1/.{ alp->alp1})==0,{w1,3}];(*Print[];
  Print["alpha=",i*10,", w1=", answ1[[1]][[2]]]*);
    Do[
    s1=(j-1)/10;  
    defx3=(defx2/.answ1)/.{alp->alp1}; 
    defy3=(defy2/.answ1)/.{alp->alp1};
    deflst[i]=Append[deflst[i],{defx3,defy3}];
If[s1==1,xt[i]=defx3];
If[s1==1,yt[i]=defy3];
(*Print[defx3,", ",defy3,","]*),{j,1,11,1}
]
,{i,1,7}]
(*結果のプロット*)
ListPlot[{deflst[1],deflst[2],deflst[3],deflst[4],deflst[5],deflst[6],deflst[7]},Frame->True,GridLines->Automatic,Joined->True]

(*座標原点(s=0)をはり先端にとっているためにこのような図になる*)

(*以下のような座標変換をする*)

(*リストのx座標を，xt[i]-x, y座標を，-yt[i]+yと変換する*)

(*ここで，xt[i],yt[i]は，はり固定端の変位(i番目の負荷荷重)*)

Do[
deflst1[i]={};
      Do[deflst1[i]=Append[deflst1[i],{xt[i]-deflst[i][[j]][[1]],-yt[i]+deflst[i][[j]][[2]]}],{j,1,11}]
,{i,1,7}];
ListPlot[{deflst1[1],deflst1[2],deflst1[3],deflst1[4],deflst1[5],deflst1[6],deflst1[7]},Frame->True,GridLines->Automatic,Joined->True]

Do[
Do[
Print[deflst1[i][[j]][[1]]," , ",deflst1[i][[j]][[2]]," ,"],{j,1,11}]


,{i,1,7}]