(*はりの大たわみ問題を微分方程式の数値解法（Ruge-Kutta法）で解く
　gz=s/L, q2= PL^2/EI *)

(*微分方程式 d^2 psi(gz)/d gz^2+q2*Cos[psi(gz)]==0 ，式（2．10）　*)

(*与えられた荷重 q のもとで，gz=0 でpsi[0]=0, psi'[0]=a と適当に初期値 a を仮定して, a を変更しながら psi'[1]=0 （先端でモーメントゼロ）となる解を探す*)

(*はじめに，適当なq2（荷重），psi1(先端たわみ角)のもとで，psi'[a]が符号を変える点を見つける*)

q2=1.69;psi1=0.1;
Do[
ans1=NDSolve[{psi''[gz]+q2*Cos[psi[gz]]==0,psi[0]==0,psi'[0]==psi1},psi,{gz,0,1}][[1]];
psi2=psi1+0.1;
ans2=NDSolve[{psi''[gz]+q2*Cos[psi[gz]]==0,psi[0]==0,psi'[0]==psi2},psi,{gz,0,1}][[1]];
graph=Plot[Evaluate[D[psi[gz]/.ans1,gz]],{gz,0,1},PlotRange->All,Frame->True,GridLines->Automatic,PlotStyle->Thick];
psid1=Evaluate[D[psi[gz],gz]/.ans1]/.gz->1;
psid2=Evaluate[D[psi[gz],gz]/.ans2]/.gz->1;
Print["psi'[0]=",Evaluate[D[psi[gz],gz]/.ans1]/.gz->0," , ","psi'[1]=",psi'[1]/.ans1," , ","psi[1]=",psi[1]/.ans1];psi1=psi2;
If[psid1*psid2<0,Break[]],{i,1,50,1}];Print[psi1-0.1," , ",psi2]
graph

(*x1, x2 はpsi'[1] が反転する初期値*)

(*次の文献の挟み撃ち法を使って初期値psi'[0]を求める－＞コンテ，ドボアー著，「数値解析と算法入門」，ブレイン図書出版*)

q2=1.69;x1=1.4;x2=1.5;
Do[
ans1=NDSolve[{psi''[gz]+q2*Cos[psi[gz]]==0,psi[0]==0,psi'[0]==x1},psi,{gz,0,1}];ans2=NDSolve[{psi''[gz]+q2*Cos[psi[gz]]==0,psi[0]==0,psi'[0]==x2},psi,{gz,0,1}];
f1=(Evaluate[D[psi[gz]/.ans1,gz]]/.gz->1)[[1]];
f2=(Evaluate[D[psi[gz]/.ans2,gz]]/.gz->1)[[1]];
x12=(x1*f2-f1*x2)/(f2-f1);
ans12=NDSolve[{psi''[gz]+q2*Cos[psi[gz]]==0,psi[0]==0,psi'[0]==x12},psi,{gz,0,1}];
f12=(Evaluate[D[psi[gz]/.ans12,gz]]/.gz->1)[[1]];
If[f1*f12<0,{x1=x1;x2=x12;f1=f1;f2=f12;},{x1=x12;x2=x2;f1=f12;f2=f2;} ];
If[Abs[x2-x1]<=1E-10,Break[]];
,{i,1,10,1}];Print["x1=",x1," , ","x2=",x2," , ","x12=",x12];
graph=Plot[Evaluate[D[psi[gz]/.ans1,gz]],{gz,0,1},PlotRange->All,Frame->True,GridLines->Automatic,PlotStyle->Thick];
psid1=Evaluate[D[psi[gz],gz]/.ans1]/.gz->1;
x=NIntegrate[Cos[Evaluate[psi[gz]/.ans1]],{gz,0,1}];y=NIntegrate[Sin[Evaluate[psi[gz]/.ans1]],{gz,0,1}];
Print["psi[1]=",(psi[1]/.ans1)[[1]]*180/Pi," , ","psi'[1]=", psid1[[1]]]
graph

(*以上を一つのプログラムにまとめる*)

(*1. q を与える（少しずつ増加させる）． 2. 根の存在範囲決定，3．その後に挟み撃ち法  4.結果の出力*)

psi1=0.1;
Print[" q^2   psi'[0]  psi[1](deg)   psi'[1]       u/L         w/L  "];
Do[

(*はじめに，psi'[1] の符号変化を捉える*)

Do[
ans1=NDSolve[{psi''[gz]+q2*Cos[psi[gz]]==0,psi[0]==0,psi'[0]==psi1},psi,{gz,0,1}][[1]];
psi2=psi1+0.1;
ans2=NDSolve[{psi''[gz]+q2*Cos[psi[gz]]==0,psi[0]==0,psi'[0]==psi2},psi,{gz,0,1}][[1]];
psid1=Evaluate[D[psi[gz],gz]/.ans1]/.gz->1;
psid2=Evaluate[D[psi[gz],gz]/.ans2]/.gz->1;
(*Print["psi'[0]=",Evaluate[D[psi[gz],gz]/.ans1]/.gz->0," , ","psi'[1]=",psi'[1]/.ans1," , ","psi[1]=",psi[1]/.ans1];*)
psi1=psi2;
If[psid1*psid2<0,Break[]]
,{i,1,50,1}];
(*Print[psi1-0.1," , ",psi2]*)

(*これで根の存在範囲が求まった．次は挟み撃ち法*)

x1=psi1-0.1; x2=psi2;

Do[
ans1=NDSolve[{psi''[gz]+q2*Cos[psi[gz]]==0,psi[0]==0,psi'[0]==x1},psi,{gz,0,1}];ans2=NDSolve[{psi''[gz]+q2*Cos[psi[gz]]==0,psi[0]==0,psi'[0]==x2},psi,{gz,0,1}];
f1=(Evaluate[D[psi[gz]/.ans1,gz]]/.gz->1)[[1]];
f2=(Evaluate[D[psi[gz]/.ans2,gz]]/.gz->1)[[1]];
x12=(x1*f2-f1*x2)/(f2-f1);
ans12=NDSolve[{psi''[gz]+q2*Cos[psi[gz]]==0,psi[0]==0,psi'[0]==x12},psi,{gz,0,1}];
f12=(Evaluate[D[psi[gz]/.ans12,gz]]/.gz->1)[[1]];
If[f1*f12<0,{x1=x1;x2=x12;f1=f1;f2=f12;},{x1=x12;x2=x2;f1=f12;f2=f2;} ];
If[Abs[x2-x1]<=1E-10,Break[]];
,{i,1,10,1}];
(*Print["x1=",x1," , ","x2=",x2," , ","x12=",x12];*)

graph=Plot[Evaluate[D[psi[gz]/.ans1,gz]],{gz,0,1},PlotRange->All,Frame->True,GridLines->Automatic,PlotStyle->Thick];
psid1=Evaluate[D[psi[gz],gz]/.ans1]/.gz->1;
x=NIntegrate[Cos[Evaluate[psi[gz]/.ans1]],{gz,0,1}][[1]];y=NIntegrate[Sin[Evaluate[psi[gz]/.ans1]],{gz,0,1}][[1]];
Print[q2," , ",x12," , ",(psi[1]/.ans1)[[1]]*180/Pi," , ", psid1[[1]]," , ",x," , ",y]
,{q2,0.5,10,0.5}];(*graph*)

(*以下に楕円積分解を示すが，Runge-Kutta 解はほとんど楕円積分解に一致*）

(*
psi0 p q^2 delta/L (L-Delta)/L) sigma0
0.,0.707107, 0,        0        , 0.     ,0.
2.5,0.722364,0.0873274,0.0290838,0.999492,0.0872831
5.,0.737277,0.175021,0.0581376,0.99797,0.174666
7.5,0.75184,0.263453,0.0871313,0.995433,0.26225
10.,0.766044,0.353003,0.116035,0.991884,0.350138
12.5,0.779884,0.444065,0.14482,0.987325,0.438437
15.,0.793353,0.537053,0.173455,0.981758,0.527256
17.5,0.806445,0.632405,0.201913,0.975188,0.616714
20.,0.819152,0.730592,0.230164,0.967617,0.706933
22.5,0.83147,0.832121,0.25818,0.959051,0.798046
25.,0.843391,0.937551,0.285935,0.949493,0.890198
27.5,0.854912,1.0475,0.3134,0.938948,0.983544
30.,0.866025,1.16264,0.340551,0.927422,1.07826
32.5,0.876727,1.28375,0.367362,0.914919,1.17453
35.,0.887011,1.41171,0.393808,0.901442,1.27257
37.5,0.896873,1.54751,0.419867,0.886997,1.37263
40.,0.906308,1.6923,0.445517,0.871586,1.47498
42.5,0.915311,1.84742,0.470738,0.855211,1.57994
45.,0.92388,2.01447,0.495511,0.837871,1.68786
47.5,0.932008,2.1953,0.51982,0.819565,1.79919
50.,0.939693,2.39217,0.543651,0.800287,1.91442
52.5,0.94693,2.60781,0.566993,0.780028,2.03417
55.,0.953717,2.84559,0.589836,0.758772,2.15915
57.5,0.96005,3.10968,0.612178,0.736498,2.29027
60.,0.965926,3.40541,0.63402,0.713174,2.42865
62.5,0.971342,3.73963,0.655369,0.688756,2.57569
65.,0.976296,4.12138,0.676242,0.66318,2.73322
67.5,0.980785,4.56289,0.696667,0.63636,2.90364
70.,0.984808,5.0812,0.716688,0.60817,3.09023
72.5,0.988362,5.70096,0.736375,0.57843,3.29761
75.,0.991445,6.45953,0.755834,0.546873,3.53254
77.5,0.994056,7.41706,0.775232,0.513086,3.80559
80.,0.996195,8.67877,0.794849,0.476389,4.13447
82.5,0.997859,10.4521,0.815193,0.435559,4.55252
85.,0.999048,13.2333,0.837329,0.38802,5.13477
*)