Chapter 8
In the Moment of Heat

Code from the Appendix

Supporting Interval Functions

mid[X_] := (Min[X] + Max[X])/2;
extremes[X_] := {Min[X], Max[X]};
hull[X_] := Interval[extremes[X]];

Count of Digit Agreement and Form for Pretty Interval Output

DigitsAgreeCount[a_, 
          b_] := (prec = 
            Ceiling@Min[Precision /@ {a, b}]; {{ad, ae}, {bd, 
                be}} = RealDigits[#, 10, prec] & /@ {a, b}; 
          If[ae ≠ be ∨ a b ≤ 0, Return[0]]; 
          If[ad == bd, Return@Length[ad]]; {{com}} = 
            Position[MapThread[Equal, {ad, bd}], False, 1, 1] - 1; 
          com);
DigitsAgreeCount[Interval[{a_, b_}]] := DigitsAgreeCount[a, b]; 
    IntervalForm[
        Interval[{a_, 
            b_}]] := (If[(com = DigitsAgreeCount[a, b]) == 0, 
          Return@Interval[{a, b}]]; 
        start = Sign[a] N[
              FromDigits@{ad[[
                    Range@com]], 1}, com]; {low, 
            up} = SequenceForm @@ Take[#, {com + 1, prec}] & /@ {ad, bd}; 
        If[ae == 0, start /= 10; ae++]; 
        SequenceForm[
          DisplayForm@SubsuperscriptBox[NumberForm@start, low, up], 
          If[ae ≠ 1, 
            Sequence @@ {" × ", 
                DisplayForm@
                  SuperscriptBox[10, ae - 1]}, ""]])

Interval Newton Iteration

IntervalNewton[f_, {a_, b_}]:=Block[{Experimental`$EqualTolerance = 0},X1 = Interval[{a, b}];X = Interval[a];While[X =!= X1,X = X1;x = Interval[mid[X]];X1=IntervalIntersection[x-((f[x])/(f'[X])),X]];X]

§8.3 Solving it Analytically: Fourier Series

§8.3.1 Solving the Transcendental Equation

θ[q_Real] := 
      2  (q^(1/4)) 
        NSum[((((-1))^k)/(2  k + 1)) (q^(k (k + 1))), {k, 
            0, ∞}]
\(\(t\_*\) == \((t /. 
          FindRoot[θ[E\^\(\(-π\^2\) t\)] == π\/\(2 \@ 5\), {t, 0.4, 0.5}, 
            AccuracyGoal -> MachinePrecision])\)\)

§8.3.2 Using Interval Arithmetic

SetAttributes[{θ0, θ1}, 
        Listable];
θ0[q_, K_] := 
      2  (q^(1/
              4)) (∑_(k = 0)^K)((((-1))^k)/(2  k + 1)) 
            (q^(k (k + 1))); 
θ1[q_, 
        K_] := ((q^(-3/4))/2) (∑_(k = 0)^K)(((-1))^k) (2  k + 1) (q^(k (k + 1))); θ[q_Interval, K_Integer] := 
      hull[θ0[q, {K - 1, K}]]; 
Derivative[1, 0][θ][q_Interval, K_Integer]:=hull[θ1[q, {K - 1, K}]];
\(Table[{k, \(t\_*\) == 
            IntervalForm@
              IntervalNewton[\((θ[Exp[\(-π\^2\) #], 
                        k] - π\/\(2 \@ 5\))\) &, {0.4, 0.5}]}, {k, 4}] // 
      TableForm\)