Quit

SetDirectory["~/Documents/dp3metric/webpage"]

discretetest[i_,  j_,  tf_ ,  de_, shift_]  =  
      Log[Abs[(( (tf[i + 1, j] - tf[i - 1, j])/
                          2 + (i - 1) (tf[i + 1, j] - 2 tf[i, j] + 
                              tf[i - 1, j])) ((tf[i, j + 1] - 
                              tf[i, j - 1])/
                          2 + (j - 1) (tf[i, j + 1] - 2  tf[i, j] + 
                              tf[i, j - 1]))  -  
                  1/16 (i - 1) (j - 
                        1) (tf[i + 1, j + 1] - tf[i - 1, j + 1] - 
                          tf[i + 1, j - 1] + tf[i - 1, j - 1])^2)/
              de^2 ]]  + 2 (tf[i, j] + shift)  -  
        Log[(i - 1) (j - 1) de^2]

eqntest[x_,  y_, tf_,  shift_] := ((( Derivative[1, 0][tf][x, y]+x Derivative[2, 0][tf][x, y])(  Derivative[0, 1][tf][x, y]+y Derivative[0, 2][tf][x, y] )- x y Derivative[1, 1][tf][x, y]^2) ) -Exp[-2 tf[x, y]  -  2 shift] x y

logeqntest[x_,  y_, tf_] := Log[Abs[((( Derivative[1, 0][tf][x, y]/x+ Derivative[2, 0][tf][x, y])(  Derivative[0, 1][tf][x, y]/y+Derivative[0, 2][tf][x, y] )- Derivative[1, 1][tf][x, y]^2) )]]+ 2 tf[x, y]

\(mygraphicaldisplay[WIDTH_, \ 
        shift_]\  := \ {
NMAX = Length[array] - 2; \ 
        Delta\  = \ WIDTH/\((NMAX - 1)\); 
For[i = 0, 
          i ≤ \ NMAX + 1, \ \(i++\), \ 
          For[j = 0, \ 
            j ≤ 
              NMAX + 1, \ \(j++\), \ {\(new[i, j] = 
                  array[\([i + 1, 
                      j + 1]\)];\)}]]; 
 (*\(newsol\  = \ 
              Interpolation[
                Table[{\((i - \((NMAX)\)\ Floor[\((i - 1)\)/NMAX] - 
                          1)\)\ Delta, \((\ Floor[\((i - 1)\)/NMAX])\) 
                      Delta, \ 
                    new[i - \((NMAX)\) Floor[\((i - 1)\)/NMAX], 
                      Floor[\((i - 1)\)/NMAX] + 1]}, \ {i, 1, NMAX^2}], \ 
                InterpolationOrder -> 
                  3];\)*) 
newsol\  = \ 
          ListInterpolation[
            array, \ {{0 - Delta, WIDTH + Delta}, \ {0 - Delta, 
                WIDTH + Delta}}, \ 
            InterpolationOrder -> 4]; 
Plot3D[
          newsol[x, y], \ {x, 0, WIDTH}, \ {y, 0, 
            WIDTH}]; 


Print["continuous error, epsilon in paper"]; 
 (*\(Print[
              Table[eqntest[i\ , \ j\ , newsol, shift], \ {i, 0.01, WIDTH, \ 
                  WIDTH/6.1}, \ {j, 0.01, WIDTH, \ 
                  WIDTH/6.1}]];\)*) 
cerror3d[NMAX]\  = \ 
          Plot3D[eqntest[x, y, newsol, shift], \ {x, 0.0, WIDTH}, {y, 0.0, 
              WIDTH}]; 
acerror[NMAX]\  = \ 
          NIntegrate[
            eqntest[x, y, newsol, shift]^2/
                0.5\ Exp[\(-2\) newsol[x, y]], \ {x, 0, 1}, {y, 0, 1}, \ 
            PrecisionGoal -> 
              2]; 
Print["average value of error squared"]; 
Print[
          acerror[NMAX]]; 

Print["log continuous error, hdot"]; 
 (*\(Print[
              Table[logeqntest[i\ , \ j\ , newsol] + 2\ shift, \ {i, 0.01, 
                  WIDTH, \ WIDTH/6.1}, \ {j, 0.01, WIDTH, \ 
                  WIDTH/6.1}]];\)*) 
logcerror3d[
            NMAX]\  = \ 
          Plot3D[logeqntest[x, y, newsol] + 2\ shift, \ {x, 0.0001, 
              WIDTH}, {y, 0.0001, WIDTH}]; 
logacerror[
            NMAX]\  = \ 
          NIntegrate[\((logeqntest[x, y, newsol] + 2\ shift)\)^2/
                0.5\ Exp[\(-2\) newsol[x, y]], \ {x, 0, 1}, {y, 0, 1}, \ 
            PrecisionGoal -> 
              2]; 
Print["average value of error squared"]; 
Print[
          logacerror[
            NMAX]]; 


 (*Print["discrete error"]; 

          Print[Table[
              discretetest[i\ , \ \ j, new, \ WIDTH/\((NMAX - 1)\), 
                shift], \ {i, Floor[NMAX/8], NMAX, \ Floor[NMAX/8]}, \ {j, 
                Floor[NMAX/8], NMAX, \ 
                Floor[NMAX/
                    8]}]]; 
Print["Average absolute value of discrete error"]; 
Print[
            Sum[Abs[discretetest[i, j, new, Delta, shift]], \ {i, 2, 
                  NMAX}, {j, 2, 
                  NMAX}]/\((NMAX)\)^2];*) 


Print["plot of discrete error"]; 
derror3d[NMAX]\  = \ 
          ListPlot3D[
            Table[discretetest[i, j, new, Delta, shift], \ {i, 2, NMAX, 
                10}, \ {j, 2, NMAX, \ 
                10}]]; 
 (*Print["divisor volume"]; 
divvol[NMAX]\  = \ 
            NIntegrate[
              D[newsol[0, y], y, y]\ y\  + \ D[newsol[0, y], y], \ {y, 0, 
                1}]; 
Print[
            divvol[NMAX]];*) 

Print["divisor volume using Riemann sum"]; 
bryi\  = 
          Floor[\((NMAX - 1)\)/WIDTH + 1/2]; 
divvol[
            NMAX]\  = \ 
          Sum[\((\((\ 
                    new[0, i - 1] - 2  new[0, i] + new[0, i + 1])\) \((i - 
                      1)\)\  + \((\ new[0, i + 1] - new[0, i - 1])\)/
                  2)\), \ {i, 1, bryi}]; 
divvol[
            NMAX]\  += \ \((\((\ 
                    new[0, bryi] - 2  new[0, bryi + 1] + 
                      new[0, bryi + 2])\) \((bryi)\)\  + \((\ 
                    new[0, bryi + 2] - new[0, bryi])\)/2)\) \((1.0 + 
                  Delta/2 - 
\((bryi)\) Delta)\)/
              Delta; 
Print[
          divvol[NMAX]]; 

Print["total volume, using interpolation function"]; 
cvol[NMAX]\  = \ 
          NIntegrate[
            x\ y\ Exp[\(-2\)\ newsol[x, y] - 2\ shift], \ {x, 0, 1}, \ {y, 0, 
              1}]; 
Print[
          cvol[NMAX]]; 

Print["total volume using Riemann sum"]; 

vol[
            NMAX]\  = \ 
          Sum[\((i - 1)\)\ \((j - 1)\) 
              Delta^4\ Exp[\(-2\) new[i, j] - 2\ shift], \ {i, 2, 
              bryi}, \ {j, 2, 
              bryi}\ ]; 

For[i = 1, 
          i <= bryi, \ \(i++\), {
\(vol[
                  NMAX]\  += \ \ 2 \((i - 
                      1)\) \((bryi)\)\ Delta^4\ Exp[\(-2\)\ new[i, 
                          bryi + 1] - 
                      2\ shift]\ \((1.0 + Delta/2 - \((bryi)\) Delta)\)/
                    Delta;\)
}]; 

vol[
            NMAX] += \ \((bryi)\)^2\ Delta^2\ Exp[\(-2\) 
                  new[bryi + 1, bryi + 1] - 
                2\ shift] \((1.0 + 
                  Delta/2 - \((bryi)\) 
                    Delta)\)^2; 

Print[
          vol[NMAX]]; 

Print["slice of continuous error"]; 
edgeslice[NMAX]\  = \ 
          Plot[eqntest[x, WIDTH, newsol, \ shift], \ {x, 0, WIDTH}, \ 
            PlotRange -> All]; 
low[NMAX]\  = \ 
          new[1, 1] + shift; 
high[NMAX]\  = \ 
          new[NMAX, \ NMAX] + shift; 
Aval[NMAX]\  = \ 
          newsol[1, 1] + 
            shift; 
Print["value of potential at (1,1) using interpolation function"]; 
Print[
          Aval[NMAX]]; 
 (*bry\  = \ 
            Floor[\((NMAX - 1)\)/1.2] + 
              1; 
xmin\  = \ \((bry - 
                  1)\)\ Delta; 
Aval1[
              NMAX]\  = \((\ 
                new[bry, 
                      bry] \((1 + \((\(-1\) + xmin)\)/Delta)\)^2\  + \ \((new[
                          bry, bry + 1] + 
                        new[bry + 1, bry])\) \((1 + \((\(-1\) + xmin)\)/
                          Delta)\) \((1 - xmin)\)/Delta
 + 
                  new[bry + 1, bry + 1] \((1 - xmin)\)^2/Delta^2)\)\  + \ 
              shift; 
Print[
            Aval1[NMAX]];*) 
}\)

<< array40;
mygraphicaldisplay[1.02, 0];

<< array80;
mygraphicaldisplay[1.02, 0];

<< array160;
mygraphicaldisplay[1.02, 0];

<< array250;
mygraphicaldisplay[1.02, 0];

conterror  =  {acerror[40],  acerror[80],  acerror[160],  
        acerror[250]}

logconterror  =  {logacerror[40],  logacerror[80], 
        logacerror[160],  logacerror[250]}

{conterror[[1]]/conterror[[2]],  
      conterror[[2]]/conterror[[3]]}

volumetable = {{ 1.2/39, 
          vol[40]}, {1.2/(80 - 1), vol[80]}, {1.2/(160 - 1),  
          vol[160]},  {1.2/249,  vol[250]}}

divisortable  =  {{1.2/39, divvol[40]}, {1.2/79,  
          divvol[80]}, {1.2/(159), divvol[160]},  {1.2/249,  
          divvol[250]}}

new11table  =  {{1.2/39, low[40]}, {1.2/(79), low[80]}, {1.2/159, 
          low[160]},  {1.2/249,  low[250]}}

newmaxtable  =  {{1.2/39,  high[40]}, {1.2/(79), 
          high[80]}, {1.2/159, high[160]},  {1.2/249,  
          high[250]}}

Avaltable  =  {{1.2/39, Aval[40]},  {1.2/79, Aval[80] }, {1.2/159, 
          Aval[160]},  {1.2/249,  Aval[250]}}

new11plot  =  
      ListPlot[new11table,  PlotStyle -> PointSize[0.03]]

new11fit  =  (a x^2  +  b  /.  
          FindFit[Join[{new11table[[3]]},  {new11table[[4]]}],  
            a x^2  +  b ,  {a, b},  x])

cornerplot  =  
      Show[{new11plot,  Plot[new11fit,  {x, 0, 0.032}]},  
        AxesLabel -> 
          TraditionalForm /@ {Δ, "f(0,0)"},  
        TextStyle -> {FontFamily -> "Times", FontSize -> 14},  
        AxesOrigin -> {0, -0.84655}]

newmaxplot  =  
      ListPlot[newmaxtable,  PlotStyle -> PointSize[0.03]]

Clear[a, b]

newmaxfit  =  (a x^2  +  b  /.  
          FindFit[Join[{newmaxtable[[3]]},  {newmaxtable[[4]]}],  
            a x^2  +  b ,  {a, b},  x])

Show[{newmaxplot,  Plot[newmaxfit,  {x, 0, 0.032}]},  
      AxesLabel -> TraditionalForm /@ {δ, "f(1.02,1.02)"},  
      TextStyle -> {FontFamily -> "Times", 
          FontSize -> 14}]

volfit  =  (a x^2  +  b  /.  
          FindFit[Join[{volumetable[[3]]},  {volumetable[[4]]}],  
            a x^2  +  b,  {a, b},  x])

volplot  =  
      Show[{ListPlot[volumetable,  PlotStyle -> PointSize[0.03]],  
          Plot[volfit,  {x, 0, 0.031}]}, 
        AxesLabel -> TraditionalForm /@ {Δ, "V"},  
        TextStyle -> {FontFamily -> "Times", FontSize -> 14},  
        AxesOrigin -> {0, 0.4997}]

divfit  =  
      a x^2  +  b  /.  
        FindFit[Join[{divisortable[[3]]},  {divisortable[[4]]}],  
          a x^2  +  b,  {a, b}, x]

\(divvolplot\  = \ 
      Show[{ListPlot[divisortable, PlotStyle -> \ PointSize[0.03]], \ 
          Plot[divfit, \ {x, 0, 0.031}]}, \ 
        AxesLabel -> 
          TraditionalForm /@ {Δ, \*""\!\(V\_d\)""}, \ 
        TextStyle -> {FontFamily -> "Times", 
            FontSize -> 14}]\)

Afit  =  
      a x^2  +  b  /.  
        FindFit[Join[{Avaltable[[3]]},  {Avaltable[[4]]}],  
          a x^2 + b,  {a, b}, x]

Aplot  =  
      Show[{ListPlot[Avaltable,  PlotStyle -> PointSize[0.03]],  
          Plot[Afit,  {x, 0, 0.031}]},  
        AxesLabel -> 
          TraditionalForm /@ {Δ, "f(1,1)"},  
        TextStyle -> {FontFamily -> "Times", 
            FontSize -> 14}]

edge80  =  
      Show[edgeslice[80],  
        AxesLabel -> TraditionalForm /@ {ξ, ""},  
        TextStyle -> {FontFamily -> "Times", 
            FontSize -> 14}]

edge160   = 
      Show[edgeslice[160],  
        AxesLabel -> TraditionalForm /@ {ξ, ""},  
        TextStyle -> {FontFamily -> "Times", 
            FontSize -> 14}]

edge250  = 
      Show[edgeslice[250],  
        AxesLabel -> TraditionalForm /@ {ξ, ""},  
        TextStyle -> {FontFamily -> "Times", 
            FontSize -> 14}]

fortyplot  =  
        Show[cerror3d[40], Mesh -> False,  
          AxesLabel -> TraditionalForm /@ {ξ,  η,  ""},  
          TextStyle -> {FontFamily -> "Times", 
              FontSize -> 14}] ;eightyplot  =  
        Show[cerror3d[80], Mesh -> False,  
          AxesLabel -> TraditionalForm /@ {ξ,  η,  ""},  
          TextStyle -> {FontFamily -> "Times", 
              FontSize -> 14}] ;onesixtyplot  =  
        Show[cerror3d[160], Mesh -> False,  
          AxesLabel -> TraditionalForm /@ {ξ,  η,  ""},  
          TextStyle -> {FontFamily -> "Times", 
              FontSize -> 14}] ;twofiftyplot  =  
        Show[cerror3d[250], Mesh -> False,  
          AxesLabel -> TraditionalForm /@ {ξ,  η,  ""},  
          TextStyle -> {FontFamily -> "Times", 
              FontSize -> 
                14}] ;