(*片持ちはりのべき級数法によるたわみ計算*)

(*5章1節　等分布荷重を受ける片持ちはりの解析*)

Clear["Global`*"]

(*級数項の上限, この値を適当に仮定して解析する．*)

nn=12

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

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

psi

(*s ははり先端から測る, 境界条件*)

eqbc1=(psi/.s->0)-psi0

eqbc2=D[psi,{s,1}]/.s->0

ans1=Solve[{eqbc1==0,eqbc2==0},{a[0],a[1]}][[1]]

(*境界条件をpsiに代入し, psiを書き換える．*)

psi=psi/.ans1

(*Governing equationへ代入．ここで, k^2=P/EI*)

eq1=D[psi,{s,2}]+k^2*Cos[psi]

(*Cos[] の項をsについて級数展開し，式全体をｓのべき級数で表現*)

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

(*この級数の s の係数を見てみる*)

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

(*a_2,..., a_nn を求めるために連立方程式を解く*)

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

psians=FullSimplify[psi/.ans1]

(*無次元変数 s1=s/L, s=L でたわみ角ゼロ*)

eq2=(psians/.{s->s1*L})/.s1->1

(*k^2=P/EI => k^2L^2=PL^2/EI => k^2L^2->P1*)

(*P1-> 無次元荷重*)

nonlin1=eq2/.{k->P1^(1/2)/L}

(*psi0 と P1の関係を求める*)

Print["psi0   P1"];
Do[
anspsi2=FindRoot[nonlin1/.{psi0->i*Pi/180},{P1,1}][[1]][[2]];
Print[i," , ",anspsi2],{i,20,80,20}]

(*参考のために楕円積分解も表示*)

(*psi0   p         q^2      delta/L   ( L-Delta)/L)
0,0.707107,8.88178*10^-16,7.45058*10^-9,0.
5,0.737277,0.175021,0.0581376,0.99797
10,0.766044,0.353003,0.116035,0.991884
15,0.793353,0.537053,0.173455,0.981758
20,0.819152,0.730592,0.230164,0.967617
25,0.843391,0.937551,0.285935,0.949493
30,0.866025,1.16264,0.340551,0.927422
35,0.887011,1.41171,0.393808,0.901442
40,0.906308,1.6923,0.445517,0.871586
45,0.92388,2.01447,0.495511,0.837871
50,0.939693,2.39217,0.543651,0.800287
55,0.953717,2.84559,0.589836,0.758772
60,0.965926,3.40541,0.63402,0.713174
65,0.976296,4.12138,0.676242,0.66318
70,0.984808,5.0812,0.716688,0.60817
75,0.991445,6.45953,0.755834,0.546873
80,0.996195,8.67877,0.794849,0.476389
85,0.999048,13.2333,0.837329,0.38802   *)

(*たわみ計算*)

(* y=Sin[psi_0]*Integrate[Cos[beta],{s,0,s}]+Cos[psi_0]*Integrate[Sin[beta],{s,0,s}]*)

(* x=Cos[psi_0]*Integrate[Cos[beta],{s,0,s}]-Sin[psi_0]*Integrate[Sin[beta],{s,0,s}]*)

(*Here, beta=a_2*s^2+a3*s^3+....*)

beta1=psi-psi0

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

(*級数展開して，項別に積分する*)

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

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

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

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

(*たわみの計算*)

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

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

(*たわみの無次元化*)

defy2=FullSimplify[((defy1/.ans1)/.{s->L*s1})/L]

defy3=defy2/.{k->P1^(1/2)/L}

defx2=FullSimplify[((defx1/.ans1)/.{s->L*s1})/L]

defx3=defx2/.{k->P1^(1/2)/L}

(* Deflections of the beam *)

(*Preparation of relation between psi0 and P1 *)

nonlin1

Plot[nonlin1/.P1->120*Pi/180,{psi0,-2,5}]

FindRoot[nonlin1/.psi0->20*Pi/180,{P1,1}][[1]][[2]]

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

Print["psi0    P1      y/L        x/L "];
deflstx={};
deflsty={};
Do[
  psi00=i*10*Pi/180;
  ansP1=FindRoot[(nonlin1/.{ psi0->psi00})==0,{P1,1}];(*Print[];
  Print["alpha=",i*10,", w1=", ansP1[[1]][[2]]]*);
     defx4=(defx3/.ansP1)/.{s1->1, psi0->psi00}; 
    defy4=(defy3/.ansP1)/.{s1->1, psi0->psi00};
    deflstx=Append[deflstx,{defx4,ansP1[[1]][[2]]}];
    deflsty=Append[deflsty,{defy4,ansP1[[1]][[2]]}];
Print[i*10," , ",ansP1[[1]][[2]]," , ",defy4," , ",defx4];
,{i,1,8}]
(*結果のプロット*)
ListPlot[{deflstx,deflsty},Frame->True,GridLines->Automatic,Joined->False,PlotStyle->PointSize[0.02],FrameLabel->{"Deflections","Load (Pl^2/EI)"}]