successive iteration

Question

I hope you can help with a problem I am having. I want to solve this equation

$$\left( n^2 \Theta^2 + 2\alpha n\Theta -1\right)R_n = \frac{1}{R_0^2}\left( 2 \sum_{m=1}^\infty R_{n+m}R_m + \sum_{m=1}^n R_{n-m}R_m\right)$$,

for different value of $\Theta$ and $\alpha$, here 0.1 and 1. This is done by guessing a trial function $R_n \propto \exp(-n)$ and doing successive iteration. Thereby calculating the sums for a fixed value. Here I choose $m \leq 256$ divided by iteration number. So first one is 256 next is 128 etc. In the end I want to sum over all the number R_n to get a single value. My idea so far is to do the following code, by brute force

    f[n_] = Exp[-n]
    E2 = (Sum[2*f[n + m]*f[m], {m, 256}] + Sum[f[n - m]*f[m], {m, 256}])/(n^2*0.1^2 + 2*n*0.1*1-1)*1/(f[0]^2)
    D1[n_] = E2; 

and then doing it all again

      E3 = (Sum[2*D1[n + m]*D1[m], {m, 256/2}] + Sum[D1[n - m]*D1[m], {m, 256/2}])/(n^2*0.1^2 + 2*n*0.1*1 - 1)*1/(D1[0]^2);
      D2[n_] = E3;

and again

      E4 = (Sum[2*D2[n + m]*D2[m], {m, 256/4}] + Sum[D2[n - m]*D2[m], {m, 256/4}])/(n^2*0.1^2 + 2*n*0.1*1 - 1)*1/(D2[0]^2);
      D3[n_] = E4;

The problem now is that the evaluation of E5 etc. take a long time. Am I doing it the right clever way? I want in the end a list of R_n values. Hope I have provided enough information.

I am trying to replicate the calculation in this paper https://journals.aps.org/prb/abstract/10.1103/PhysRevB.36.1931, from equation 7.1, where they for $\alpha=1$ and $\theta=0.1$ find that

$$ s=2\pi/\Theta \left( 1+1/8 \sum_{n,m=-\infty}^{\infty} R_{n+m} R_n R_m \right)=22.2$$

Answer 1

I am not completely sure but I think, you are looking for the roots of the function F below. For given R, F[R] is the difference of both side of the equation

$$\left( n^2 \Theta^2 + 2\alpha n\Theta -1\right)R_n = \frac{1}{R_0^2}\left( 2 \sum_{m=1}^\infty R_{n+m}R_m + \sum_{m=1}^n R_{n-m}R_m\right)$$

Caution: R[[1]] corresponds to $R_0$, R[[2]] to $R_1$, R[[3]] to $R_2$ and so on.

nn = 100;
R0 = Exp[-Range[0., nn - 1]];
Θ = 0.1;
α = 1.;

F = R \[Function] 
   Join[{(-1) R[[1]] - 2. R[[1 ;;]].R[[1 ;; nn]] / R[[1]]^2}, 
    Table[(n^2 Θ^2 + 2 α n Θ - 1) R[[n + 1]] - (2. R[[n + 1 ;;]].R[[1 ;; nn - n]] + R[[2 ;; n + 1]].R[[n ;; 1 ;; -1]])/R[[1]]^2, {n, 1, nn - 1}]
    ];

We can employ FindRoot to search for roots as follows

R = Array[r, nn];
FindRoot[Evaluate[F[R]], Evaluate@Transpose[{R, Evaluate@R0}],
 MaxIterations -> 1000]

In a nutshell, FindRoot applies Newton's method

$$R^{(k+1)} = R^{(k)} - DF(R^{(k)})^{-1} \, F(R^{(k)})$$

Unfortunately, $R^{(k)}$ seems to converge towards the constant vector $0$ for $k \to \infty$.

Edit

If I got it correctly, the fixed point iteration they use in the article is generated by the following function:

G = R \[Function] 
   Join[{ 2. R[[1 ;;]].R[[1 ;; nn]]/(-1 R[[1]]^2)}, 
    Table[(2. R[[n + 1 ;;]].R[[1 ;; nn - n]] + 
        R[[2 ;; n + 1]].R[[n ;; 1 ;; -1]])/((n^2 Θ^2 + 
          2 α n Θ - 1) R[[1]]^2), {n, 1, 
      nn - 1}]];

This is how you can optain the list of the first 1000 iterates:

sequence = NestList[G, R0, 1000];

And indeed, the iterates seem to converge to a fixed point (the differences converge to the zero vector):

ListPlot[Max /@ Abs[Differences[sequence[[1 ;; 100]]]], PlotRange -> All]

I am surprised.