Graphics Settings

<< Miscellaneous`Units`;
Remove[Gray];
<< Graphics`Master`;
Off[General::spell1];
fs = 7;
fsbig = 9;
width  =  Convert[8.64 Centi Meter,  Point]/Point;
th  =  0.0008;

Figure 9.1

Clear[p, f, q];
p[α_]  :=  
        2  +  Sin[10 α];
f[x_,  α_]  :=  (x^α) 
          Sin[α/(2 - x)];
q[α_,  
        n_] := (1/n) (∑_(k = 1)^(2 n))f[((k - 0.5)/n), α]
II[α_]  :=  p[α] q[α];

dataq  =  
        Table[{α, q[α, 20001]},  {α, 0,  5,  
            0.02}];

datap  =  
        Table[{α, p[α]},  {α, 0,  5,  
            0.02}];
dataI  = 
        Transpose[{ 
            First  /@  
              dataq, ( Last  /@  dataq) ( 
                Last  /@  datap)}];

p1  =  
        ListPlot[dataq, DisplayFunction  -> Identity, 
          PlotJoined -> 
            True, 
PlotStyle -> {Thickness[th], 
              Dashing[{ .007}]}];
p2  =  
        Plot[p[α],  {α,  0,  5}, 
          DisplayFunction  -> Identity, MaxBend -> 1, 
          PlotPoints -> 
            100, 
PlotStyle -> {Thickness[th], 
              Dashing[{ .007}]},  
          DisplayFunction  ->  Identity];

p3  =  
        ListPlot[dataI, DisplayFunction  -> Identity, 
          PlotJoined -> 
            True, 
PlotStyle -> {Thickness[
                th]}];

Show[p1, p2,p3, Frame -> True, 
Axes -> False, PlotRange -> {{- .001, 5},  {-1.001,  
                3.45}}, Epilog -> {Text["p(α)",{2.26,  3.23},TextStyle ->  {FontFamily -> "Times", 
                        FontSize -> fsbig} ],
Text["I(α)",{2.35,   .6} ],
\(Text["α", {4.95, \ \(-1\)}, \ {1, \(-1\)}\ ]\),
Text["q(α)",{1.3,   .5} ]},
FrameTicks -> {Range[0, 5],  
              Range[-1,  3.5,  1],  {{π/4, "π/4"}},  
              None},
DisplayFunction  -> $DisplayFunction, TextStyle -> {FontFamily -> "Times", 
              FontSize -> fs},ImageSize -> width];

Figure 9.2a

f[x_,  α_]  :=  (x^α) 
          Sin[α/(2 - x)];

    
Plot[f[x,  π/4],  {x, 0, 2},  Frame -> True,  
        Axes -> False,  
        PlotRange -> {{- .001, 2},  {-2.001, 2.001}}, 
        FrameTicks ->   {Range[0,  2,   .5], {-1, 0, 1, 2}, 
            None, None},  PlotPoints -> 1000, MaxBend -> 1, 
        PlotStyle -> Thickness[th], 
        TextStyle -> {FontFamily -> "Times", FontSize -> fs}, 
        ImageSize -> width];

Figure 9.2b

Clear[f1, f2, z, dz, t];
f1[x_Real,  α_]  :=  (x^α) 
        Sin[(α/(2 - x))]
f1[z_, α_] := - I (E^(α (Log[z] + (I/(2 - z)))))
z[t_] := t + I t (2 - t);
dz[t_] = 
        1 + 2 (1 - t) I;
f2[t_, α_] := Re[ f1[z[t], α]   dz[t]]

Plot[f2[t, π/4],  {t, 0, 2},  Frame -> True,  
        Axes -> None,  MaxBend -> 1,  
        FrameTicks -> {Range[0,  2,   .5],  Range[0,  1,   .25],  
            None,  None},  PlotStyle -> Thickness[th], 
        TextStyle -> {FontFamily -> "Times", FontSize -> fs}, 
        ImageSize -> width];

Figure 9..3

Clear[II];
II[α_]  :=  II[α] = 
        N[NIntegrate[(2 + Sin[10 α]) Re[
                f2[t,  α]],  {t, 0, 2},  
            PrecisionGoal -> 13,  AccuracyGoal -> 13],  
          MachinePrecision]

α0 = (785933674/(10^9)); h = (1/(10^10)); 
    ymax = 3.03373258648549;
data = 
        Table[{n, II[α0 + h n]}, {n, -500, 500, 
            1}];

data[[-2]]

ListPlot[data1 = 
            data /. {n_, 
                  y_} :> 
{0.785933674 + 
                    n  10^-10, y - ymax},

PlotRange -> {{0.78593361, 
                0.78593374}, 
{-10^-12.85,  
                10^-13.6}},Frame -> True, Axes -> False,
\(TextStyle -> {FontFamily -> "Times", FontSize -> \ fs}\),FrameTicks->{{{0.78593362400001`,"0.785933624"},{0.785933674, "0.785933674"},
{0.785933724, "0.785933724"}},
Join[{{0, "(Subscript[(I), max])"}},  {-10^-ToExpression@#, "(Subscript[(I), max])-10"^StringForm["-``", #]}& /@ {13,  "13.5"}],None,None},
PlotStyle -> 
            PointSize[ .003;  .005],
GridLines -> {{}, {}},ImageSize ->    width, 
AspectRatio ->  .6,
Epilog -> 
{Text["(Subscript[(I), max])=3.03373258648549", {0.785933615,  10^-13.7},{-1, 1}, \(TextStyle -> \ {FontFamily -> "Times", FontSize -> \ \ \ fsbig - 1}\)],
\(Text["α", \ {0.785933732, \ \(-10^\(-\((12.88)\)\)\)}, {\(-1\), 0}, \ 
                  TextStyle -> \ {FontFamily -> "Times", 
                      FontSize -> \ \ fsbig}]\)}];

Figure 9.4

Clear[f1, f2, z, zder];
f1[x_Real, α_] := 
      (x^α) Sin[(α/(2 - x))]
f1[z_, α_] := -I (E^(α (Log[z] + (I/(2 - z)))))
z[t_] := t + I t (2 - t);
dz[t_] = 1 + 2 (1 - t) I;
f2[t_, α_] := Re[f1[z[t], α] dz[t]]
f2C[t_, α_] := 
      f1[z[t], α] dz[t]
f2der[t_, α_] := 
      Re[f1[z[t], α] (Log[
                z[t]] + (I/(2 - z[t]))) dz[
            t]]
f2derC[t_, α_] := 
      f1[z[t], α] (Log[z[t]] + (I/(2 - z[t]))) dz[
          t]

Plot[f2der[t, (π/4)], {t, 0, 2}, Frame -> True, 
        Axes -> False, GridLines -> {{}, {}}, 
        PlotRange -> {{- .0001, 2.001}, {-2.001, 1.5}}, 
        MaxBend ->  .5, PlotPoints -> 100, 
        PlotStyle -> Thickness[th], 
        FrameTicks -> {Range[0, 2,  .5], {-1, 0, 1, 2}, None, 
            None},  TextStyle -> {FontFamily -> "Times", 
            FontSize -> fs}, ImageSize -> width];

There is very fast exponential decay here. I believe that if t = 2 - (1/δ), then the value is about \(TraditionalForm\`e\^\(\(-\(2\/5\)\) α\ δ\)\). If δ=200, this is (10^(-28. ))

Plot[
        Log[10,  
          Abs[(2 + Sin[10  π/4]) 
                f2der[t, π/4] + 
  10 
                 Cos[10  π/4] f2[t, π/4]]],  {t, 0, 1.99},  
        Frame -> True,  Axes -> False,  PlotRange -> All,  
        PlotPoints -> 200,  MaxBend -> 2];

Clear[obj];
obj[α_Real] := 
      NIntegrate[(2 + Sin[10 α]) f2der[t, α] + 
          10 Cos[10 α] f2[t, α], {t, 0, 1.999}, 
        AccuracyGoal -> 10]

<<NumericalMath`GaussianQuadrature`

Clear[wts]
wts[i_] := wts[i]  =  
          Transpose[
            GaussianQuadratureWeights[200,   i  .4,  
              i  .4  +   .4]];
Do[wts[i],  {i, 0,  4}]
wts1  =  
        Transpose[
          GaussianQuadratureWeights[200,  0, 
            2]];

func1[
          t_,  α_]  := (2 + Sin[10 α]) f2der[
              t, α] + 
          10 Cos[10 α] f2[t, α];

    
obj1[α_Real] := 
      wts1[[2]] .  (func1[#, α] & /@  
            wts1[[1]])
obj2[α_Real] := 
      Sum[wts[
              i][[2]] .  (func1[#, α] & /@  wts[
                  i][[1]]),  {i, 0, 4}]

obj1[π/4. ; ans200]
obj2[Pi/4. ; ans200]

Plot[obj1[α],  {α,  0, 1}];

(αans = α /. 
                FindRoot[obj2[α],  {α, 7/10, 8/10}] ;) // 
        Timing // First
NumberForm[αans, 11]
N[ans200, 12]

(αans = α /. 
                FindRoot[obj[α],  {α, 7/10, 8/10}] ;) // 
        Timing // First
NumberForm[αans, 16]
ans200

Figure 9.5

x[t_] := ((2 t)/(1 + t));
dx[t_] := (2/(((1 + t))^2));
\(\(Plot[
        Evaluate[{x[t]\^α\ dx[t]\ Sin[
                  1\/2\ α\ \((1 + 
                        t)\)], 
1\/2\ \((1 + 
                      t)\)\ x[t]\^α\ dx[t]\ Cos[
                    1\/2\ α\ \((1 + t)\)] + 
Log[
                    x[t]]\ x[t]\^α\ dx[t]\ Sin[
                    1\/2\ α\ \((1 + 
                          t)\)]} /. α -> π\/4], {t, 0, 100}, PlotRange -> {{\(-2\), 100}, {\(- .22\),  .4}}, 
        Frame -> True, 
        Axes -> 
          False, 
FrameTicks -> {Range[0, 100, 20], \ 
            Range[\(- .2\), \  .4,  .2], \ None, \ None}, 
        PlotPoints -> 100, 
        PlotStyle -> {{Thickness[th]}, {Thickness[th], 
              Dashing[{ .01}]}}, \ 
        TextStyle -> {FontFamily -> "Times", FontSize -> fs}, 
        ImageSize -> width];\)\)

Figure 9.6

q::usage  =  "q[α] gives the q function as defined in the \
chapter. q = I / (2 + sin(10 α)).";

dI::usage  =  "dI[α] uses the term-by-term derivative of \
the G-series to give the derivative of I(α) reduced by a certain \
constant since we are interested in root-finding only. The formula used is \
efficient, but fails at integer values. Thus this is meant to be used on the \
interval (0,1).";

dItrue::usage  =  "dItrue fixes dI to give I'(α) by \
returning the constant.";

\(\(dqtrue::usage\  = \ "dqtrue uses dItrue to get q'(α).";\)
\)
qderRich::usage  =  "qderRich[α,h] uses one level of \
Richardson extrapolation on a central difference to estimate q' \
numerically.";


\(\(q[α_Real]\  := \ 
        Module[{prec = Precision[α], 
            z = α\^2/16}, 
s = \(t = 
              0\); 
stermpow = \(zpow = \(sterm = \(tprod = 
                  1\)\)\); 
tsum = \(tsumprod = 
              1\  - \ Log[α/2] - \ 2\ EulerGamma - 
                PolyGamma[α + 1]\); 
recip = 
            1/α; 
          n = 0; 
While[
            Abs[stermpow]\  > \ 10^\(-prec\)\  || \ 
              Abs[tsumprod]\  > \ 
                10^\(-prec\), \(n++\); 
αnprod2 = \((α - 2\ n + 2)\)\ \((α - 2\ n + 
                    1)\); αnprod1 = \((α - 
                    2\ n)\)\ \((α - 2\ n + 
                    1)\); 
αnprod21 = \(2  α - 4\ n + 3\)\/αnprod2\ ; 
s += \ stermpow; \ \ t += 
              zpow\ tsumprod; \ \ zpow\  *= \ \(-z\); 
sterm *= αnprod1\/\(n\ \((n + 1)\)\ \((2\ n + 1)\)\^2\); \ \ stermpow = 
              sterm\ zpow; \ \ tprod *= αnprod2\/\(n\^2\ \((4  n\^2 - 1)\)\); \ \ tsum += αnprod21 + 1\/n + 1\/\(2\ n - 1\) + 
                1\/\(2\ n + 1\); 
tsumprod = 
              tsum\ tprod]; 
\(2\^α\) α \((2\ π\ z\ s + t)\)];\)
\)
\(nterms[α_, prec_] := 
      If[prec ≤ 15, 6, 
        2 + Floor[
            n /. FindRoot[
                N[α]\^\(2  n\) == \(16\^n\) \(\((\(n!\))\)\^2\) 
                    10\^\(-prec\), {n, 2, 2000}]]]
\)
dI[α_Real] := 
        Module[{prec = Precision[α], 
            z = (α^2)/
                16}, 
s = t = stermderlog = tprodderlog = 0; sterm = 2 π z; tprod = 1; 
          recip = 1/α; 
tsum = tsumprod = 
              1  -  Log[α/2] -  2 EulerGamma - 
                PolyGamma[α + 
                    1]; 
der = tsumder = -recip - 
                PolyGamma[
                  1, α + 
                    1];  

Do[{αprod1, αprod2} = (α - 2 n + 1) (α - 
                    2 n + {0, 
                      2}); 
αprod21 = ((2  α - 4 n + 3)/α)prod2; 
            s +=  sterm;   t +=  
              tsumprod;  
sterm *= -(( αprod1  z     )/(n (n + 
                          1) (((2 n + 1))^2)));  
            stermderlog  +=  ((2  α - 4 n + 1)/α)prod1 ;  
tprod *= -((z αprod2)/((n^2) (4  (n^2) - 
                          1)));  
            tsum += αprod21 + (1/n) + (1/(2 n - 1)) + 
                (1/(2 n + 1)); 
tsumprod = tsum  tprod; 
            tsumder  -= (((((α - 2 n + 2))^2) + (((α - 2 n + 1))^2))/α)(prod2^2); 
tprodderlog  += αprod21; 
der +=  
              tsumprod (2 n recip  + tprodderlog + (tsumder/tsum)) + 
                sterm (2  n recip + 
                      stermderlog),  
{n, 
              nterms[α, 
                prec]}]; 
(2 + 
                  Sin[10 α]) (recip + 
                  Log[2] + ((2  s  recip + der)/(     s + t))) + 10 Cos[
                10  α] ];

 (* 
      dItrue is the correct derivative of I (α);  
      it can be used to get the correct derivative of q (α) *) 

dItrue[α_Real] := 
        dI[α] q[α];

dqtrue[α_Real] := ((dItrue[α] - 10 Cos[10 α] q[α])/(2 + Sin[10 α]));
qder[α_, 
        h_] := ((q[α + h] - q[α - h])/(2 h))
qderRich[α_, h_] := 
        (4/3) qder[α, h] - (1/3) qder[α, 2 h];

α = π/4. ; 
    data = Table[{-i, 
          Log[10, Abs[
              qderRich[α, (10^(-i))] - dqtrue[α]]]}, {i, 5, 
          2, - .1}];

\(\(ListPlot[data, PlotJoined -> True, \ 
        Epilog -> {PointSize[ .015], \ Point\  /@ \ data}, \ 
        Frame -> 
          True, \ 
FrameTicks -> {\({#, \ "10"\^StringForm["-``", \(-#\)]} &\)\  /@ \ 
              Range[\(-5\), \ \(-2\)], \ \({#, \ "10"\^StringForm["-``", \(-#\)]} &\)\  /@ 
              Range[\(-13\), \ \(-8\), 1], None, \ 
            None}, 
TextStyle -> {FontFamily -> "Times", \ FontSize -> \ \ fs}, \ 
        PlotRange -> {{\(-5.1\), \(-1.9\)}, \ {\(-13.4\), \ \(-7.99\)}}, \ ImageSize -> width];\)\)