This Mathematica notebook gives a very elementary way to think about the Gradient Newton Galerkin Algorithm (GNGA).  This hopefully will illuminate the “polynomial f” section in the original paper by Neuberger and Swift.

Define ψ1 to be the first eigenfunction on the Laplacian on the interval [0, π], normalized in L2.

ψ1[x_] = (Sqrt[(2/π)])Sin[x]

(∫_0^π)((ψ1[x])^2)\[DifferentialD]x

ψ2[x_] = (Sqrt[(2/π)])Sin[2x]
ψ3[x_] = (Sqrt[(2/π)])Sin[3x]

Define u[x] to be the general function in the 3-dimensional Galerkin space spanned by the first 3 eigenfunctions.  The Fourier coefficients a1, a2 and a3 are coordinates on the Galerkin space.

u[x_] = a1 ψ1[x] + a2 ψ2[x]+a3 ψ3[x]

The eigenvalues of the Laplacian are {1, 4, 9, ⋯}

u''[x]

We want  the PDE evaluated at a function in the Galerkin space, which is a function of x, to be identically 0.

pde[x_] = u''[x] + s u[x] + ((u[x])^3)

But it can’t be 0 for a function in the Galerkin space.  The best we can do is that the projection of the PDE onto the Galerkin space is 0.

Note that g1, g2, and g3 are each a cubic polyomial in the variables a1, a2 and a3.  This is because the nonlinearity in the PDE is cubic.

g1[s_, a1_, a2_, a3_] = (∫_0^π)pde[x] ψ1[x]\[DifferentialD]x
g2[s_, a1_, a2_, a3_] = (∫_0^π)pde[x] ψ2[x]\[DifferentialD]x
g3[s_, a1_, a2_, a3_] = (∫_0^π)pde[x] ψ3[x]\[DifferentialD]x

To approximate a solution, fix a value of s, and use the FindRoot[ ] function with a guess for the Fourier coefficients.  In the saved notebook file the guess is {1, 0, 0}, but you can change sval and the guess.

sval = 0;soln = FindRoot[{g1[sval, a1, a2, a3] == 0, g2[sval, a1, a2, a3] == 0, g3[sval, a1, a2, a3] == 0} , {{a1,1}, {a2,0}, {a3, 0}}]

Note that a1, a2, and a3 continue to be unassigned variable.  Use the “subtitution” feature “/.” to plot the solution.

Plot[u[x]/. soln, {x, 0, π}]

Plot[pde[x]/.soln/.s-> sval, {x, 0, π}]

Plot[{u[x], pde[x]}/. soln/.s-> sval, {x, 0, π}]