Mathematica Solutions to the 
ISSAC System Challenge 
1997

Wolfram Research, Inc.

Problem 1

Method 1: Calculate the ratio of the eigenvalues explicitly, using brute force (the straightforward way)

The condition number for the matrix under consideration is the ratio of the largest to the smallest eigenvalue.

The condition number of a Hilbert matrix of dimension n is roughly equal to (e^(3.5 n))(see [1]); for n = 256 this is approximately (10^389). We expect to multiply and subtract expressions in intermediate stages, so to be on the safe side we'll use about twice as many digits.

Hilbert[n_] := Table[(1/(i + j - 1)), {i, n}, {j, n}]; 
    evals = Eigenvalues[N[Hilbert[256], 800]];

Here is a plot of the 256 eigenvalues.

ListPlot[Log[10, N[evals]], Frame -> True];

This is the answer.

N[First[evals]/Last[evals], 10]

Method 2: Use the power method for the matrix and its inverse (the fast way)

Fortunately, there is a closed-form expression for the matrix elements of the inverse of a Hilbert matrix (see [2]). This allows us to iterate toward the smallest and the largest eigenvalues using the power method.

Hilbert[n_] := Table[(1/(i + j - 1)), {i, n}, {j, n}]

\(InverseHilbert[n_] := 
	
      Table[\((\(-1\))\)\^\(i + j\)\ \((i + j - 1)\)\ 
          Binomial[n + i - 1, n - j]\ Binomial[n + j - 1, n - i]\ 
          Binomial[i + j - 2, i - 1]\^2, {i, n}, {j, n}]\)

As an aside, we note that the elements of this matrix have an interesting distribution.

\(\(ListDensityPlot[InverseHilbert[128], 

      	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
      Mesh -> False, \ ColorFunction -> \((Hue[0.7  #]&)\)]; \)\)

Again, about 800 digits should be enough.

{h, ih} = N[{Hilbert[256], InverseHilbert[256]}, 800];

We define two functions, 𝒩 to normalize an eigenvector and RayleighQuotient to calculate an eigenvalue.

𝒩[v_List] := (v/(Sqrt[v . v]))

RayleighQuotient[u_List, v_List] := 
      (((u . h . u) (v . ih . v))/((u . u) (v . v)))

We iterate until the ratio of the two eigenvalues remains constant up to the required precision.

fp = 
      FixedPoint[
        Function[x, 
          {{#1, #2}, RayleighQuotient[#1, #2]}&@@{
              𝒩[
                h . x[[1, 1]]], 
              𝒩[
                ih . x[[1, 2]]]}], {
          Table[1, {2}, {256}], 2}, 
        SameTest -> 
          (Abs[1 - #1[[2]]/#2
                        [[2]]] < 
                (10^(-6))&)];

Here are plots of the components of these two eigenvectors.

Show[
      GraphicsArray[
        ListPlot[#, PlotRange -> All, PlotJoined -> True, 
              Frame -> True, Axes -> False, 
              DisplayFunction -> Identity]&/@
          fp[[1]]]];

This is the ratio of the two eigenvalues.

N[fp[[2]], 6]

Method 3: Use LUDecomposition (the general way)

This method works for arbitrary matrices; we don't have to know the inverse explicitly.

Hilbert[n_] := Table[(1/(i + j - 1)), {i, n}, {j, n}]

An LU-decomposition is more stable than an eigenvalue calculation, so we don't have to use 800 digits. It is enough to have somewhat more digits than the expected size of the condition number.

lud = LUDecomposition[h = N[Hilbert[256], 500]];

\(estimateNorm2[mat_?MatrixQ] := 
	
      Module[{l = Length[mat], prec = Precision[mat], x, y}, \ \ \ \ \ \ 

        		\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
        y = Table[Random[Real, {1, 2}, prec + 5], {l}]; 

        \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 		
        For[j = 1, j ≤ 8, \(j++\), x = y; y = mat . x]; 

        \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 	
        \@\(y . Conjugate[y]\/x . Conjugate[x]\)]\)

\(estimateInverseNorm2[{mat_?MatrixQ, perm_, estcondnum_}] := 
      Module[{l = Length[mat], prec = Precision[mat], x, y}, \ \ \ \ 

        		\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
        x = Table[Random[Real, {1, 2}, prec], {l}]; 

        \ \ \ \ \ \ \ \ \ \ \ \ \ \ 		
        For[j = 1, j ≤ 8, \(j++\), y = x; \ \ \ \ \ \ 

          			\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
          x = LUBackSubstitution[{mat, perm, estcondnum}, y]]; 

        	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
        \@\(x . Conjugate[x]\/y . Conjugate[y]\)]\)

N[estimateNorm2[h] estimateInverseNorm2[lud], 5]

Problem 2

Method 1: Use a built-in function

The most natural approach gives the correct result but with a warning message.

NIntegrate[(x^((x^x))), {x, 1, 6}]

To guarantee the required digits and to avoid the message we use arbitray-precision arithmetic, allow more recursive subdivisions at the right endpoint, and use a more appropriate integration method.

NIntegrate[(x^((x^x))), {x, 1, 6}, MaxRecursion -> 20, 
      WorkingPrecision -> 30,  PrecisionGoal -> 7, 
      Method -> DoubleExponential]

Method 2: Substitute to get a nicer integrand

We make the substitution (x^x)⟶z. Using Solve[(x^x) == z, x], this is equivalent to x⟶log(z)/ProductLog(log(z)). This makes the integrand much better behaved, getting rid of any subdivision problems.

NIntegrate[
      Evaluate[(((Log[z]/ProductLog[Log[z]]))^z) 
          (∂_z) Log[z]/ProductLog[Log[z]]], {z, 1, (6^6)}]

Problem 3

Method: Use a built-in function

The most natural approach gives the first seven digits but also generates messages that tell us that the tail estimate failed.

\(NSum[1\/\@\(n\^π + n\^2 + n\^\(\@2\) + 1\)\%3, {n, 1, ∞}]
        \  // N[#, 14]&\)

To guarantee 14 digits we use bignum arithmetic (as recommended in the error message). Because we know that the sum converges, we don't need to verify the convergence numerically.

\(NSum[1\/\@\(n\^π + n\^2 + n\^\(\@2\) + 1\)\%3, {n, 1, ∞}, 
      VerifyConvergence -> False, WorkingPrecision -> 120, 
      PrecisionGoal -> 16, Method -> NIntegrate]\)

Problem 4

Method 1: Calculate the residue numerically

The coefficient of (z^n)in the polynomial p(z) can be obtained as the residue of p(z) (z^(-n - 1))at z = 0. (This old and widely used technique is used to calculate large coefficients in enumeration theory [3].)

The default precision goal of NIntegrate is only 6, so we cannot be sure of the last digits.

\(N[Im[
          NIntegrate[
            Evaluate[
				
              \((\(\((x + 1)\)\^2000\ \((x\^2 + x + 1)\)\^1000\ 
                        \((x\^4 + x\^3 + x\^2 + x + 1)\)\^500\)\/x\^3001 /. 
                    x -> E\^\(I\ φ\))
                  \)\ I\ 
                E\^\(I\ φ\)], 

            		{φ, 0, 2\ π}]]/\((2  π)\), 15]\)

We use bignum arithmetic to make sure to get 13 correct digits. We also use the trapezoidal rule because the integrand is periodic.

\(N[Im[
          NIntegrate[
            Evaluate[
				
              \((\(\((x + 1)\)\^2000\ \((x\^2 + x + 1)\)\^1000\ 
                        \((x\^4 + x\^3 + x\^2 + x + 1)\)\^500\)\/x\^3001 /. 
                    x -> E\^\(I\ φ\))
                  \)\ I\ 
                E\^\(I\ φ\)], 

            		{φ, 0, 2\ π}, 	Method -> Trapezoidal, 
            WorkingPrecision -> 400, PrecisionGoal -> 16, 
            MaxRecursion -> 10]]/\((2  π)\), 15]\)

Method 2: Calculate the residue symbolically

Residue[
      (((x + 1))^2000) ((((x^2) + x + 1))^1000) 
        ((((x^4) + (x^3) + (x^2) + x + 1))^500)/(x^3001), {x, 0}]

N[%, 14]

Method 3: Calculate the series up to order 3000

ser = 
      Series[(((x + 1))^2000) ((((x^2) + x + 1))^1000) 
          ((((x^4) + (x^3) + (x^2) + x + 1))^500), {x, 0, 3000}];

SeriesCoefficient[ser, 3000]

N[%, 14]

Method 4: Expand the polynomial

Since the polynomial is a square, we find the appropriate coefficient in the square of the expansion of the square root.

Sum[#[[i]] #
              [[Length[#] + 1 - i]], {i, 
            1, Length[#]}]&[
      CoefficientList[
        Expand[(((x + 1))^1000) ((((x^2) + x + 1))^500) 
            ((((x^4) + (x^3) + (x^2) + x + 1))^250)], x]]

N[%, 14]

Problem 5

Method 1: Find the root numerically

A well-known upper bound for the largest root (Subscript[z, ν]^((ν)))of the (v^th) Laguerre polynomial is 4ν+2 (see [4, 5]). The Laguerre polynomials satisfy the well-known differential equation x (y^′)′(x)+(1 - x)y'(x)+n y(x)=0.

\(\(\(DSolve[
            x\ \(\(y'\)'\)[x]\  + \ \((1 - x)\) \(y'\)[x]\  + \ ν\ y[x] 
              == 0, \ y[x], x] // 
	FunctionExpand\) // Simplify\) // 
      ExpandAll\)

Transforming the equation to Liouville normal form using \(TraditionalForm
      \`\(y(x)\)⟶\(\@x\) \(e\^\(x/2\)\) 
        \(u(x)\)\) we get ((u^′))'(x)+(-(x^2) + 4 ν x + 2 x + 1) /(4  (x^2)) 
              u(x)=0. When -(x^2) + 4 ν x + 2 x + 1 > 0, the solution u(x) will oscillate and will therefore have zeros. When -(x^2) + 4 ν x + 2 x + 1 < 0, the solution u(x) will increase or decrease monotonically. The larger zero, 
      2 ν + (Sqrt[4 ν (ν + 1) + 2]) + 1 = 
        4 ν + 2 + O((1/ν)), approximates the largest zero of (Subscript[L, 1000])(x).

This gives us a starting point ((Subscript[z, start]) = 4 ν + 2)  for the numerical root finding procedure. We slide down the curve on the strictly concave side toward the root.

We define the function f to be the Laguerre polynomial and the gradient g of f. For faster and more reliable numerical calculation we avoid evaluating these functions as explicit polynomials.

f[x_?NumberQ] := LaguerreL[1000, SetPrecision[x, 100]]g[x_?NumberQ] := -LaguerreL[999, 1, SetPrecision[x, 100]]

FindRoot[Evaluate[f[x]], {x, 4002}, WorkingPrecision -> 100, 
        AccuracyGoal -> 1000, MaxIterations -> 50, 
        Jacobian -> g[x]] // N[#, 12]&

The secant method gives the same result.

FindRoot[Evaluate[f[x]], {x, 4002, 4003}, WorkingPrecision -> 100, 
        AccuracyGoal -> 1000, MaxIterations -> 100] // 
      N[#, 12]&

We check numerically to see that this is really the largest root.

\(RootCount[ρ_] := 
	
      Round[Im[NIntegrate[
              \(-\(\(LaguerreL[999, 1, 
                        ρ\ 
                          E\^\(I\ φ\)]\ 
                      \)\/LaguerreL[1000, 
                        ρ\ 
                          E\^\(I\ φ\)]
                        \)\) I\ ρ\ 
                E\^\(I\ φ\), {
                φ, 0, π, 2\ π}, 
              WorkingPrecision -> 100, PrecisionGoal -> 4]]/
          \((2  π)\)]\)

{RootCount[3940], RootCount[3950]}

Using real root isolation techniques from the package Algebra`RootIsolation` we can prove that this is the largest root. We can use the same package to find an interval with rational endpoints containing the root to the required precision.

Method 2: Calculate all roots one after another

We start at the origin and step toward the largest root by explicitly calculating all the roots on the way.

\(f[x_?NumberQ] := LaguerreL[1000, SetPrecision[x, 100]]; 

    g[x_?NumberQ] := \(-LaguerreL[999, 1, SetPrecision[x, 100]]\); \)

{rootf[1], rootf[2]} = 
      (x /. 
              FindRoot[f[x], {x, #}, WorkingPrecision -> 100, 
                AccuracyGoal -> 1000, MaxIterations -> 50, 
                Jacobian -> g[x]])&/@{0, 1/100};

\(Do[rootf[i] = 
        x /. FindRoot[f[x], 

            				{x, rootf[i - 1] + 1/2 \((rootf[i - 1] - rootf[i - 2])\), 
              
		\ \ \ \ \ \ \ \ \ \ \ \ \ 
              rootf[i - 1] + 3/2 \((rootf[i - 1] - rootf[i - 2])\)}, 
            WorkingPrecision -> 100, AccuracyGoal -> 1000, 
            MaxIterations -> 50], 
	\ \ \ \ {i, 3, 1000}]\)

This is the largest root.

N[rootf[1000], 14]

Here is the cumulative root density as well as the spacing between the roots.

Show[
      GraphicsArray[{
          Graphics[{Table[Line[{{i, 0}, {i, rootf[i]}}], {i, 1000}]}, 
            Frame -> True], 
          Graphics[{
              MapIndexed[
                Point[{#2[[1]], #1}]&, 
                -Apply[Subtract, 
                    Partition[Table[rootf[i], {i, 1000}], 2, 1], {1}]]}, 
            Frame -> True]}]];

Problem 6

Method: Use a built-in function

The default term ordering for GroebnerBasis is the lexicographic ordering, so the input is straightforward.

\(\(gb = 
      GroebnerBasis[

        		{8\ w\^2 + 5\ w\ x - 4\ w\ y + 2\ w\ z + 3\ w + 5\ x\^2 + 
            2\ x\ y - 7\ x\ z - 7\ x + 7\ y\^2 - 8\ y\ z - 7\ y + 7\ z\^2 - 
            8\ z + 8, 
			
          3\ w\^2 - 5\ w\ x - 3\ w\ y - 6\ w\ z + 9\ w + 4\ x\^2 + 2\ x\ y - 
            2\ x\ z + 7\ x + 9\ y\^2 + 6\ y\ z + 5\ y + 7\ z\^2 + 7\ z + 5, 

          			\(-2\)\ w\^2 + 9\ w\ x + 9\ w\ y - 7\ w\ z - 4\ w + 8\ x\^2 + 
            9\ x\ y - 3\ x\ z + 8\ x + 6\ y\^2 - 7\ y\ z + 4\ y - 6\ z\^2 + 
            8\ z + 2, 
			
          7\ w\^2 + 5\ w\ x + 3\ w\ y - 5\ w\ z - 5\ w + 2\ x\^2 + 9\ x\ y - 
            7\ x\ z + 4\ x - 4\ y\^2 - 5\ y\ z + 6\ y - 4\ z\^2 - 9\ z + 2}, {
          x, y, z, w}]; \)\)

gb

The basis consists of a univariate polynomial in w and three linear polynomials in x, y, and z.

Exponent[#, {x, y, z, w}]& /@ gb

A numerical version of the Gröbner basis can be calculated more than 150 times faster.

\(\(gbN = 
      GroebnerBasis[

        		{8\ w\^2 + 5\ w\ x - 4\ w\ y + 2\ w\ z + 3\ w + 5\ x\^2 + 
            2\ x\ y - 7\ x\ z - 7\ x + 7\ y\^2 - 8\ y\ z - 7\ y + 7\ z\^2 - 
            8\ z + 8, 
			
          3\ w\^2 - 5\ w\ x - 3\ w\ y - 6\ w\ z + 9\ w + 4\ x\^2 + 2\ x\ y - 
            2\ x\ z + 7\ x + 9\ y\^2 + 6\ y\ z + 5\ y + 7\ z\^2 + 7\ z + 5, 

          			\(-2\)\ w\^2 + 9\ w\ x + 9\ w\ y - 7\ w\ z - 4\ w + 8\ x\^2 + 
            9\ x\ y - 3\ x\ z + 8\ x + 6\ y\^2 - 7\ y\ z + 4\ y - 6\ z\^2 + 
            8\ z + 2, 
			
          7\ w\^2 + 5\ w\ x + 3\ w\ y - 5\ w\ z - 5\ w + 2\ x\^2 + 9\ x\ y - 
            7\ x\ z + 4\ x - 4\ y\^2 - 5\ y\ z + 6\ y - 4\ z\^2 - 9\ z + 2}, {
          x, y, z, w}, CoefficientDomain -> InexactNumbers]; \)\)

Short[gbN, 12]

Exponent[#, {x, y, z, w}]& /@ gbN

This shows that the numerically calculated Gröbner basis and the exact one are the same.

Chop[
      Expand[(gb/(#[[1, 1]]&))/@
              MonomialList[gb, {x, y, z, w}]] - gbN, (1/(10^30))]

Problem 7

Method 1: Calculate the first nonvanishing term in the series expansion

\(f[x_] = 
	
      Tan[Tanh[Sin[x]]] + Tanh[Sin[Tan[x]]] + Sin[Tan[Tanh[x]]] - 
        Tan[Sin[Tanh[x]]] - Sin[Tanh[Tan[x]]] - Tanh[Tan[Sin[x]]] - 
        Tan[Sinh[Tanh[x]]] - Sinh[Tanh[Tan[x]]] - Tanh[Tan[Sinh[x]]] + 
        Tan[Tanh[Sinh[x]]] + Tanh[Sinh[Tan[x]]] + Sinh[Tan[Tanh[x]]]; 


    g[x_] = Sinh[Tanh[Sin[x]]] + Tanh[Sin[Sinh[x]]] + Sin[Sinh[Tanh[x]]] - 
        Sinh[Sin[Tanh[x]]] - Sin[Tanh[Sinh[x]]] - Tanh[Sinh[Sin[x]]] - 
        Tan[Sinh[Sin[x]]] - Sinh[Sin[Tan[x]]] - Sin[Tan[Sinh[x]]] + 
        Tan[Sin[Sinh[x]]] + Sin[Sinh[Tan[x]]] + Sinh[Tan[Sin[x]]]; \)

We calculate the first nonvanishing terms in the series for f(x) and g(x) at x = 0.

fSeries = Series[f[x], {x, 0, 25}] // Normal

gSeries = Series[g[x], {x, 0, 25}] // Normal

The limit follows easily.

((fSeries[[1]] /. 
        x -> 
          gSeries[[1]])/(gSeries
          [[1]] /. 
        x -> fSeries[[1]]))

N[%, 12]

Method 2: Calculate the limit numerically

We approach the origin with the sequence (10^(-i)), stopping when we reach a stationary result.

\(f[x_] = 
	
      Tan[Tanh[Sin[x]]] + Tanh[Sin[Tan[x]]] + Sin[Tan[Tanh[x]]] - 
        Tan[Sin[Tanh[x]]] - Sin[Tanh[Tan[x]]] - Tanh[Tan[Sin[x]]] - 
        Tan[Sinh[Tanh[x]]] - Sinh[Tanh[Tan[x]]] - Tanh[Tan[Sinh[x]]] + 
        Tan[Tanh[Sinh[x]]] + Tanh[Sinh[Tan[x]]] + Sinh[Tan[Tanh[x]]]; 


    g[x_] = Sinh[Tanh[Sin[x]]] + Tanh[Sin[Sinh[x]]] + Sin[Sinh[Tanh[x]]] - 
        Sinh[Sin[Tanh[x]]] - Sin[Tanh[Sinh[x]]] - Tanh[Sinh[Sin[x]]] - 
        Tan[Sinh[Sin[x]]] - Sinh[Sin[Tan[x]]] - Sin[Tan[Sinh[x]]] + 
        Tan[Sin[Sinh[x]]] + Sin[Sinh[Tan[x]]] + Sinh[Tan[Sin[x]]]; \)

\(oldRatio = 10\^9; i = 2; 

    While[newRatio = 
        N[f[g[N[10\^\(-i\), 200\ i]]]\/g[f[N[10\^\(-i\), 200\ i]]], 12]; \ 

      	\ \ \ \ \ \ \ \ \ \ \ 
      Abs[\(newRatio - oldRatio\)\/newRatio] > 10\^\(-12\), \ 

      	\ \ \ \ \ \ \ \ \ \ i = i + 1; oldRatio = newRatio]; 
newRatio\)

Problem 8

Method: Do a change of variables

The current version cannot do the integral directly. (However, the next version can.)

(∫_0^1) (x^2) PolyLog[3, (1/(x + 1))] \[DifferentialD]x

In such a situation it is always a good idea to simplify the argument of the most complicated function. So we let 1/( x + 1)⟶z, and this works.

(∫_1^((1/2)))(((((1 - z)/z)))^2) PolyLog[3, z] 
        (-(1/(z^2))) \[DifferentialD]z

The substitution ( x + 1)⟶z also works.

(∫_1^2)(((z - 1))^2) PolyLog[3, (1/z)] 
        \[DifferentialD]z

Note 1: We check the integral numerically to make sure the symbolic result is correct.

\(N[{%, 
	
        NIntegrate[x\^2\ PolyLog[3, 1\/\(x + 1\)], {x, 0, 1}, 
          PrecisionGoal -> 12]}, 12]\)

Note 2: The corresponding indefinite integrals can be done, too.

∫(((((1 - z)/z)))^2) PolyLog[3, z] 
            (-(1/(z^2))) \[DifferentialD]z // FunctionExpand // 
      Simplify

∫(((z - 1))^2) PolyLog[3, (1/z)] \[DifferentialD]z // 
        FunctionExpand // Simplify

Problem 9

Method 1: Solve a corresponding discrete version

𝒦[x_, y_] := 
      (E^(x + y + (x^2) + (y^2) + x y + (x^2) (y^2)));

We get the abscissas and weights for the Gaussian quadrature of order 32 in the interval (0,1) to 20 digits. Any order between, say, 10 and 40 would also work. Going to 20 digits makes sure that arbitrary-precision arithmetic is used; this controls the errors in the calculations.

Needs["NumericalMath`GaussianQuadrature`"]

{𝓅, 𝓌} = 
      Transpose[Re[GaussianQuadratureWeights[(2^5), 0, 1, 20]]];

We iterate the eigenfunction expansion with an increasing number of subdivisions until we have the desired precision (actually i = 1 is already enough).

\(eMaxOld = 1; i = 0; 

    While[c = 2\^i; 
\ 
      eMaxNew = 
        \(Eigenvalues[
            \(Outer[𝒦, #, #]&\)[
                Flatten[Table[𝓅 + s, {s, 0, c - 1}]/c]] . 

              				\ \ \ \ \ \ \ \ \ \ \ \ \ 
              DiagonalMatrix[Flatten[Table[𝓌, {c}]/c]]]
            \)[[1]]; 
      Abs[\(eMaxNew - eMaxOld\)\/eMaxNew] > 10\^\(-14\), i = i + 1; 
      eMaxOld = eMaxNew]; \ N[eMaxNew, 14]\)

Method 2: Iterate toward the largest eigenvalue

For this kernel we can iterate the application of the integral operator starting with some vector that has a nonvanishing component in the direction of the largest eigenvalue (see [6]). With increasing iterations we increase the precision, too.

𝒦[x_, y_] := 
      (E^(x + y + (x^2) + (y^2) + x y + (x^2) (y^2)));

\(\(setInterpolationQuality[p_, io_, ag_, \ ip_] := 
	
      SetOptions[FunctionInterpolation, InterpolationPoints -> p, 
        InterpolationOrder -> io, AccuracyGoal -> ag, 
        InterpolationPrecision -> ip]; \)\)

i = 0;Off[NIntegrate::inum];
fpl=FixedPointList[
	      
		 Function[f,setInterpolationQuality[i = i + 1; 20 + 10 i, i + 3,  i, 
                    16];
		    
		  a = FunctionInterpolation[
                      NIntegrate[
                        Evaluate[
                          SetPrecision[
                            𝒦[x, y] 
                              f[[1]][
                                y], 25]], {y, 0, 1}, PrecisionGoal -> i, 
                        WorkingPrecision -> 25], {x, 0, 1}];
		    
		   b = FunctionInterpolation[a[x]/a[1/2], {x, 0, 1}];
			  
			{b, Plus@@Table[a[x]/b[x], {x, 1/100, 1, 1/100}]/100}
                    ],{1&, 1}, SameTest -> 
              (Abs[(#1[[2]] - #2
                              [[2]])/#2
                          [[2]]] < 
                    (10^(-14))&)];

This is the largest eigenvalue.

N[fpl[[-1, -1]], 14]

Here are the iterated eigenfunctions; only the first few are visibly different.

Plot[Evaluate[First[#][x]&/@fpl], {x, 0, 1}, 
      PlotStyle -> Table[Hue[h], {h, 0, 0.7, 0.7/14}]];

Here are the local logarithmic differences between the iterated eigenvectors and the last eigenvector. The dip in the middle comes from the normalization condition.

Off[Plot::plnr]; 
    Plot[Evaluate[
        Log[10, 
              Abs[First[#][x] - 
                  fpl[[-1, 1]][
                    x]]]&/@Rest[Drop[fpl, -1]]], {x, 0, 1}, 
      PlotStyle -> Table[Hue[h], {h, 0, 0.7, 0.7/12}], 
      Frame -> True];

This shows how the eigenvalues converge.

Last /@ fpl

This shows the convergence of the differences between the iterated eigenvalues and the last iteration.

(fpl[[-1, 2]] - #)&/@Drop[%, -1]

Problem 10

Method 1: Solve the ODE for (y^′)(x) and (x^′)(y)

Let's start with the straightforward numerical solution.

nsolNaive=NDSolve[{y'[x]==z[x],z'[x]==(x^3) + y[x]^3 + z[x]^3,y[0] == 0,z[0] == 0},{y[x], z[x]},{x, 0, 2}];

Let's look at the graphs of the function and its derivative. There seems to be a singularity of (y^′)(x) on the interval (0,2) at the point r ≈ 1.65.

Show[
      GraphicsArray[
        Apply[Plot[
              Evaluate[#1 /. nsolNaive], {x, 0, 
                nsolNaive[[1, 1, 2, 0, 1, 1, 2
                  ]]}, DisplayFunction -> Identity, 
              PlotLabel -> #2, PlotStyle -> {Thickness[0.012]}]&, {{
              y[x], "y[x]"}, {z[x], "y'[x]"}}, {1}]]];

If (y^′)(x) has a singularity at x = r, then near this singularity the differential equation reduces to \(TraditionalForm
      \`\(y\^′′\)(x) = \(\(y\^′\)(x)\)\^3\).

DSolve[y''[x] == y'[x]^3, y[x], x][[1
        ]] /. 
      f_ Power[s_, 1/2] :>  Sqrt[Expand[(f^2) s]]

The solution of this differential equation has the form y(x) = (Subscript[c, 1]) + (Sqrt[r - 2  x]). We conclude that y(x) is finite at x = r and that there is a square root singularity at x = r.

We first solve for y(x) up to some point to the left of the singularity and then switch to the differential equation for x(y). We can then solve across the singularity for y(x); it is not a singularity for x(y). We set the appropriate options in NDSolve to achieve the required precision goal.

nsol1=NDSolve[{y'[x]==z[x],z'[x]==(x^3) + y[x]^3 + z[x]^3,y[0] == 0,z[0] == 0},{y[x], z[x]},{x, 12/10, 12/10},StartingStepSize -> (10^(-3)),WorkingPrecision -> 40,PrecisionGoal -> 30,AccuracyGoal -> 20,MaxSteps -> 3000];

Now we change from y[x] to x[y] and continue solving.

{y0 = 
        y[x] /. nsol1[[1]] /. 
          x -> 12/10, 
      z0 = z[x] /. nsol1[[1]] /. 
          x -> 12/10};

nsol2=NDSolve[{-x''[y]==(x[y]^3 + (y^3)) ((x'[y])^3)+1,x[y0] == 12/10,x'[y0]==1/z0},x[y],{y, y0, 1},WorkingPrecision -> 40,PrecisionGoal -> 30,AccuracyGoal -> 20,MaxSteps -> 3000];

To get an approximation to y(r) we use the calculated solution.

y[x] /. nsolNaive[[1]] /. 
      x -> 
        nsolNaive[[1, 1, 2, 0, 1, 1, 2
          ]]

We use the resulting interpolating function to find the minimum of -x (y), which is the point r.

FindMinimum[
        Evaluate[
          -nsol2[[1, 1, 2]]], {y, 
          9/10, 1}, WorkingPrecision -> 40, PrecisionGoal -> 50, 
        AccuracyGoal -> 100]  //  N[#, 14]&

Here is a plot of the real part of y(x) near the square root branch point. The violet line is the path of integration.

Method 2: Change to the polar coordinate system centered at {0,1}

We change variables from x, y to ρ, φ.

x = ρ[φ] Cos[φ];
y=ρ[φ] Sin[φ] + 1;

y'=(∂_φ)y/(∂_φ)x; y''=(∂_φ)y'/(∂_φ)x;

This is the new differential equation in ρ(φ).

\(newOde = 
      \((Factor[Numerator[Together[\(y'\)' - \((x\^3 + y\^3 + y'\^3)\)]]] 
          == 0)\)\)

nsol=NDSolve[{newOde,ρ[-(π/2)] == 1,ρ'[-(π/2)]==0},ρ,{φ,-(π/2),(π/4)},WorkingPrecision -> 30,PrecisionGoal -> 30,AccuracyGoal -> 20,MaxSteps -> 5000];

Here is a plot of the solution.

ParametricPlot[Evaluate[{ρ[φ] Cos[φ],ρ[φ] Sin[φ] + 1}/.nsol],{φ,-(π/2),(π/4)},AxesLabel -> {"x", "y"}];

Now the solution always stays in the order 1, so there is not an accuracy problem.

The maximum of \(TraditionalForm
      \`\(-\(ρ(φ)\)\)\ \(cos(φ)\)\) determines r.

MapAt[N[#, 15]&,FindMinimum[Evaluate[-ρ[φ] Cos[φ]/.nsol[[1]]],{φ,-1/10,1/10},WorkingPrecision -> 30,AccuracyGoal -> 80],1]

Substituting back into the solution of the differential equation gives y(r).

ρ[φ] Sin[φ] + 1 /. %[[2
            ]] /. 
        nsol[[1]]  //  N[#, 15]&

All timings and memory measurements were done on a 200-MHz Pentium Pro running Linux with the standard shipping 3.0.1 Mathematica kernel.

({{problem, method, seconds, megabytes}, {1, 1, 33000, 48}, {1, 2, 1000, 63}, {1, 3, 1500, 33}, {2, 1, 3, 1}, {2, 2, 8, 1}, {3, 1, 32, 1}, {4, 1, 180, 1}, {4, 2, 270, 3}, {4, 3, 3200, 5}, {4, 4, 3200, 5}, {5, 1, 13, 1}, {5, 2, 2400, 2}, {6, 1, 160, 2}, {7, 1, 12, 1}, {7, 2, 14, 1}, {8, 1, 6, 3}, {9, 1, 40, 1}, {9, 2, 280, 1}, {10, 1, 280, 1}, {10, 2, 190, 2}})

[1] J. Todd, Appl. Math., Ser. 39, NBS 1954.

[2] N. J. Higham, Accuracy and Stability of Numerical Algorithms, SIAM 1996.

[3] A. Odlyzko in R. Graham, M. Grötschel, L. Lovász (eds.), Handbook of Combinatorics, Elsevier 1995.

[4] E. R. Neumann, Jahresber, DMV 1921, Teubner, Leipzig 1921.

[5] A. Zarzo, J. S. Dehesa, J. Comp. Appl. Math. 50, 613 (1994).

[6] R. P. Kanwal, Linear Integral Equations, Birkhäuser 1997.