GNGA:  y'' + (y^3) = 
        0,   y(0) = 0 = y(1).

\(n\  = \ 100; \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ m\  = \ 
          20; \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 

dx\  = \ 
          1/n // N; \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \  (*\ delta\ x\ *) 

        x\  = \ Table[i/n, \ {i, \ 0, \ n}]\  // 
            N; \ \ \ \ \ \ \ \ \ \ \ \ \ \  (*\ 
          Grid\ *) 

 (*\ 
          nonlinearity\ *) 
\(f[u_]\ \ \  := 
            u^3;\)\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 

        \(fp[u_]\  := \ 3  u^2;\)

        \(F[u_]\ \ \  := u^4/4;\)

 (*\ Fourier\ Basis, \ 
          eigenfunctions\ of\ \  - y'' = \ λy, \ 
          normalized\ in\ L\^2\ *) 

        \(p[k_, \ x_] := \ Sqrt[2] Sin[k\ Pi\ x]\  // N;\)\ 

        \(lam\  = \ Table[k^2\ Pi^2, \ {k, \ 1, \ m}] // N;\)\ 

 (*\ 
          Variational\ Functional\ J\ *) 

        \(J[u_, a_] := 
            Sum[a[\([j]\)]^2  lam[\([j]\)]/2, {j, 1, m}] - 
              Sum[\(F[u]\)[\([i]\)], {i, 1, n + 1}]*
                dx;\)

 (*\ 
          gradient\ of\ J\ *) 

        \(g[u_, \ a_] := \ 
            Module[{v\  = \ f[u], \ j\  = \ 0}, 
Table[
                a[\([j]\)]\ lam[\([j]\)]\  - \ \((p[j, \ x]\  . \ v)\)*
                    dx, \ {j, \ 1, \ m}]];\)

 (*\ 
          Hessian\ of\ J\ *) 

        \(kron[b1_, b2_]\  := \ 
            If[b1 == b2, \ 1, \ 0];\)\ 

        \(h[u_, \ a_] := \ 
            Module[{v\  = \ fp[u], \ j\  = \ 0, \ 
                k\  = \ 0}, 
Table[
                kron[j, \ 
                      k]\ lam[\([j]\)]\  - \ \((\((p[k, \ x]*
                            p[j, \ x]\ )\) . \ v)\)*dx, \ {j, \ 1, \ 
                  m}, \ {k, \ 1, \ m}]];\)

 (*\ 
          Number\ of\ negative\ e - vals\ of\ Hessian\ *) 

        \(sig[elist_] := 
            Module[{cnt = 0, i = 0}, 
For[i = 1, 
                i ≤ 
                  Length[elist], \(i++\), 
\(If[
                    elist[\([i]\)] < \(- .01\), \(cnt++\), 
                    Null];\)
]; 
cnt
];\)\)

\(\(\( (*\ 
      initial\ guess\ of\ Fourier\ coeffs\ *) \)\(
\)\(\(a\  = \ 
        Table[0, \ {k, \ 1, \ m}];\)

    \(a[\([3]\)]\  = \ 7.5;\)

    \(u\  = \ 
        Sum[a[\([k]\)] p[k, \ x], \ {k, \ 1, \ m}]\  // 
          N;\)

 (*\ 
      Gradient\ and\ Hessian\ for\ first\ iteration\ *) 

    \(grad\  = \ g[u, \ a];\)

    \(hess\  = \ h[u, \ a];\)\)\)\)

\(go = \ 1; \ tol\  = \ 10^\((\(-10\))\); \ num\  = \ 0; \ 
    maxits\  = \ 100; 
    delta\  = \ 1.0;

 (*\ 
      do\ Newton\ iterations\ until\ grad\ small\ *) \)
\(\(While[
        go == 
          1, 
\(num++\); 

 (*\ 
          take\ step\ in\ search\ direction\ \ *) 
a\  = \ 
          a\  - \ delta\ *\ PseudoInverse[hess] . grad; 
	
\ \ \ \  (*\ 
          new\ Fourier\ expan\ with\ new\ coeffs\ *) 
u\  = \ Sum[a[\([k]\)] p[k, \ x], \ {k, \ 1, \ m}]\  // 
            N; 

\ \ \ \  (*\ 
          calc\ Gradient\ and\ Hessian, \ 
          evals\ of\ Hessian\ *) 
grad\  = \ 
          g[u, \ a]; 
hess\  = \ 
          h[u, \ a]; 
evals\  = \ 
          Eigenvalues[hess]; 
	
\ \ \ \  (*\ 
          calc\  || \n  abla\ J || _  2, \ test\ for\ convergence, \ 
          output\ info\ *) 
nrm\  = \ 
          Norm[grad]/
            n; 
If[\((nrm\  < \ 
                tol)\)\  || \ \((num\  ≥ \ maxits)\), 
          go = 0]; 
Print["iteration ", 
          num, \ \ \ \ \ \ \ \ ", Norm =  ", 
          nrm, 
\ \ \ \ \ \ \ \ \ \ \ ", J = ", \ \ \ \ \ \ \ \ \ J[u, a], \*"", ||u\!\(\( || \_∞\)\) = "", 
          Max[Abs[u]], ", MI = ", 
          sig[evals]];

];\)\)

a[[3]]  =  7.5;

ListPlot[u,  PlotJoined -> True]

Print[Chop[a, 10^-10]];

Plot[{J[tt*u, tt*a]}, {tt, -1.5, 1.5}]

p1 = p[1, x];
a  =  Table[0,  {k,  1,  m}];
a[[1]]  =  1;
Plot[{J0[tt*p1, tt*a]}, {tt, 0, 5.5}]

p1 = p[2, x];
a  =  Table[0,  {k,  1,  m}];
a[[2]]  =  1;
Plot[{J0[tt*p1, tt*a]}, {tt, 0, 10.5}]

gamma[u_, a_] := 
        Sum[a[[i]]^2  lam[[i]], {i, 1, m}] - (u^4) . one/
            n;

Sum[a[[i]]^2  lam[[i]], {i, 1, m}] - (u^4) . one/
        n

Plot[{gamma[tt*u, tt*a], J0[tt*u, tt*a]}, {tt, 0, 1.1}]

J0[u, a]

ListPlot[z[-10, 5], PlotJoined -> True];

Sort[Eigenvalues[h[v, a]]]

z[s_, t_] := s*p[1, x] + t*p[2, x];
plist = {};
For[s = -3, s ≤ 3, 
        s +=  .5, 
Print[
          s]; 
row = {}; 
For[t = -3, 
          t ≤ 3, 
          t +=  .5, 
v = z[s, t]; 
a = 
            Table[p[k, x] . v/n, {k, 1, m}]; 
eval = Sort[
                Eigenvalues[h[v, a]]][[1]]; 
If[
            eval > 0, 
row = 
              Append[row, 1], 
row = 
              Append[row, 
                0]
];
]; 
plist = 
          Append[plist, row];
];

plist

ListPlot3D[plist]

plist[[1]]

z[s_, t_] := s*p[1, x] + t*p[2, x];
v = z[1, 2];
a = Table[p[k, x] . v/n, {k, 1, m}]

w[x_] := E^x;   
plist = {};
For[s = -3, s ≤ 3, 
        s +=  .5, 
Print[s]; 
For[
          t = -3, t ≤ 3, 
          t +=  .5, 
v = z[s, t]; 
a = 
            Table[p[k, x] . v/n, {k, 1, m}]; 
eval = Sort[
                Eigenvalues[
                  hess = h[v, a]]][[1]]; 
If[
            eval > 0, 
ww =  
              Table[(w[x]*v) . p[k,  x]/n,  {k,  1,  m}];   (* 
              coeffs of inhomog term *) 
grad = 
              g[v, a] + 
                ww; 
chi = -Inverse[hess] . 
                  grad; 
c1 = chi[[1]]; 
            c2 = chi[[2]]; 
plist = 
              Append[plist, {{s, t}, {c1, 
                    c2}}]
];
];
];

ppp = 
        Graphics[{PointSize[ .05], 
            Point[{1.6913264205803369`, 0.301203645978688`}]}];

<<Graphics`PlotField`

aaa = ListPlotVectorField[plist, ScaleFactor -> 1/2]

Show[ppp, aaa]

plist

save

perturb newt

\(n\  = \ 
      100; \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \  (*\ 
      grid\ size\ *) 

    dx\  = \ 1/n // N; \ \ \ \ \ \ \ \ \ \ \ \ \ \ \  (*\ delta\ x\ *) 

    f[u_] := u^3; \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \  (*\ 
      nonlinearity\ *) \)
fp[u_] :=  3  u^2;
F[u_] := u^4/4;
\(\(\(w[x_] := ccc*E^x;\)\(
\)
    \)\)
\(kron[b1_, b2_]\  := \ If[b1 == b2, \ 1, \ 0]; \ \ \  (*\ 
      kronecker\ delta\ *) 

    m\  = \ 10; \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \  (*\ 
      number\ of\ Fourier\ Modes\ in\ Expansion\ *)  (*\ Grid, \ 
      table\ of\ x_  0, \ x_  1, \ 
      x_  2, \ \(\(...\)\(\ \)\(x_n\)\(\ \)\(values\)\)\ *) 

    x\  = \ Table[i/n, \ {i, \ 0, \ n}]\  // 
        N;\ \ \ \ \ \ \ \ \ \ \ \ \ \ 

\  (*\ Fourier\ Basis, \ 
      family\ of\ functions, \ eats\ x - grid, \ 
      spits\ out\ y - values\ *) 
 (*\ 
      eigenfunctions\ of\  - \(y'\)' = \(\(\\)\(lambda\)\)\ y, \ 
      normalized\ in\ L^2\ s . t . \ \int_  0^1\ p^2\  = \ 1\ *) \)
\(\(p[k_, \ x_] := \ Sqrt[2] Sin[k\ Pi\ x]\  // N;\)\ 
\)
\(one\  = \ Table[1, \ {i, \ 0, \ n}]; \ \ \ \ \ \  (*\ stupid, \ 
      integrate/Riemann\ sums\ using\ dot\ with\ one\ *) 


    lam\  = \ Table[k^2\ Pi^2, \ {k, \ 1, \ m}] // N; \ \  (*\ 
      eigenvalues\ of\  - \(y'\)' = \(\(\\)\(lambda\)\)\ y\ *) 

 (*\ 
      gradient\ of\ J\ *) 

    g[u_, \ a_] := \ 
      Module[{v\  = \ f[u], \ j\  = \ 0}, 
Table[
          a[\([j]\)]\ lam[\([j]\)]\  - \ \((\((p[j, \ x]\ *\ v + 
                      w[x]*p[j, x])\) . one)\)*dx, \ {j, \ 1, \ m}]];

 (*\ 
      Hessian\ of\ J\ \((derivative\ of\ g)\)\ *) \)
\(\(h[u_, \ a_] := \ 
        Module[{v\  = \ fp[u], \ j\  = \ 0, \ 
            k\  = \ 0}, 
Table[
            kron[j, \ 
                  k]\ lam[\([j]\)]\  - \ \((\((p[k, \ x]*p[j, \ x]\ *\ v)\) . 
                    one)\)*dx, \ {j, \ 1, \ m}, \ {k, \ 1, \ m}]];\)

 (*\ 
      counts\ number\ of\ neg\ evals\ of\ h, \ signature, \ 
      number\ of\ concave\ down\ directions\ of\ critical\ point\ *) \)
sig[elist_] := 
        Module[{cnt = 0, i = 0}, 
For[i = 1, 
            i ≤ 
              Length[elist], i++, 
If[
                elist[[i]] < - .01, cnt++, 
                Null];
]; 
cnt
];

J[u_, a_] := 
        Sum[a[[j]]^2  lam[[j]], {j, 1, m}]/
            2 - (F[u] + w[x]*u) . one/n;
gamma[u_, a_] := 
        Sum[a[[i]]^2  lam[[i]], {i, 1, m}] - (u^4 + w[x]*u) . one/
            n;

a

a  =  Table[0,  {k,  1,  m}];

\(\(\( (*\ 
      positive\ sol\ *) \)\(
\)\(\(a\  = \ Table[0, \ {k, \ 1, \ m}];\)

    \(a[\([1]\)]\  = \ 2.5;\)

    \(a[\([3]\)] = \(-0.11\);\)\)\)\)

ccc = 1.5;

\(\(\( (*\ 
      Main\ Newton\ Loop\ *) \)\(
\)\(
\)\(
\)\(For[
        ccc = 0, ccc ≤ 10, 
        ccc +=  .25, 

 (*\ 
          Fourier\ expansion\ *) 
u\  = \ 
          Sum[a[\([k]\)] p[k, \ x], \ {k, \ 1, \ m}]\  // N; 

        norm\  = \ 1; \ tol\  = \ 10^\((\(-10\))\); \ num\  = \ 0; \ 
        maxits\  = \ 20; 
        delta\  = \ 1.0; 

 (*\ 
          do\ Newton\ iterations\ until\ grad\ is\ small\ or\ too\ many\ iterations\ *) 

        While[\((norm\  > \ tol)\)\  && \ \((num\  < 
                maxits)\), 
\(num++\); 

 (*\(If[norm <  .5, delta =  .5];\)	*) 
 (*\ 
            calc\ hess, \ grad, \ evals\ of\ hess, \ 
            and\ take\ step\ in\ search\ direction\ *) 
hess\  = \ h[u, \ a]; 
evals\  = \ 
            Eigenvalues[hess]; 
grad\  = \ 
            g[u, \ a]; 
a\  = \ 
            a\  - \ delta\ *\ PseudoInverse[hess] . grad; 
	
	 (*\ 
            new\ Fourier\ expan\ with\ new\ coeffs\ *) 
u\  = \ Sum[a[\([k]\)] p[k, \ x], \ {k, \ 1, \ m}]\  // N; 
	
	 (*\ 
            calc\  | 
              g | \(\(^\)\(2\)\)\ and\ print\ with\ its\ and\ signature\ *) 
norm\  = \ 
            grad . grad; 
 (*\(Print[{num, \ norm, \ 
                  sig[evals]}];\)*) 
]; 

Print[{num, \ norm, \ 
            sig[evals]}]; 
Print["J = ", 
          J[u, a], \ ", Max[u] = ", Max[u]]; 
ListPlot[
          u, PlotJoined -> True]; 
Print[ccc]; 
Print[a]; 

        Print[Eigenvalues[hess]]; 

        Plot[{J[tt*u, tt*a]}, {tt, 0, 
            1.1}]

];\)\)\)

a

ListPlot[u, PlotJoined -> True];
Print[ccc];
Print[a];
Print[Eigenvalues[hess]];
Plot[{J[tt*u, tt*a]}, {tt, 0, 1.1}]

gamma[u, a]

Plot[gamma[tt*u, tt*a], {tt, 0, 2}]

Plot[w[x], {x, 0, 1}, PlotRange -> {0, 4}]

ccc

u

a

g[u, a]