Computing a Sequence using Recurrence

I posted a similar question earlier today, but unfortunately I realized after several hours that I over-simplified the problem for this forum. There is a very good answer to the Computing a sequence by recurrence question. I just can't do the simplification step that turns the sums into $Q_{N-1}$'s. This is because the form of my problem is not as simple as the question I posed.

I am trying to create a list where each element of the list is calculated using the sum of all previous elements of the list and also uses the $N$th element of another list. The initial value of the list is a known quantity. $\alpha$, $\beta$, $\gamma$, $\delta$, $\phi$, and $\theta$ are all known constants $k_N$ is a known list of $N$ elements.

$\quad \quad Q_N=\alpha\frac{\beta-\sum_{i=0}^{N-1}{Q_i\delta}}{k_N}\sqrt{(\frac{k_N}{\gamma}\frac{\beta+\sum_{i=0}^{N-1}{Q_i\delta}}{\beta-\sum_{i=0}^{N-1}{Q_i\delta}})^{\phi}-(\frac{k_N}{\gamma}\frac{\beta+\sum_{i=0}^{N-1}{Q_i\delta}}{\beta-\sum_{i=0}^{N-1}{Q_i\delta}})^{\theta}},\,N>0$ $\quad \quad Q_{N=0}=Q_0$

I'm not quite sure how to code this. I think the Nest function is somewhat close to what I want, but I need to keep plugging in the sum of all previous values, not just the previous value.

FoldList and Nest seem to be close to what I'm trying to do, but I can't figure it out. It may be relevant that $N$ is going to be rather large. The context is a numeric solution to the time evolution of a system. So I want to have smaller time-steps, which leads to a large number of elements.

Any help would be greatly appreciated.

Note: I'm not asking the forum to figure out this rather complicated form necessarily; I just wanted to show the form to say that the simplification made to use FoldList in the other answer won't work here. I think the computation requires using all previous values, and not just the previous value. Answering the other question just without using the simplification may suffice.

Computing a sequence by recurrence

• How long is your sequence, ie, how large of a value of $N$ are you looking for? Note that defining $R_N=\sum_{i=0}^{N-1}Q_i$, one has $R_{N+1}=R_N+Q_N$, and the expression for $Q_N$ only depends on $R_N$. So as long as $N$ is a couple million or less, you can just compute both $Q_N$ and $R_N$ at each stage, store the results, and then compute $Q_{N+1}$. – DumpsterDoofus Jan 9 '15 at 23:32
• If you're computing as large as $N$ in the tens of millions or billions, then you may have to think more carefully, but the above method is efficient in the sense that each new term is computed in $O(1)$ time, and thus is probably asymptotically optimal. – DumpsterDoofus Jan 9 '15 at 23:33
• I think an upper bound of 1 million for N. – matt Jan 10 '15 at 0:54
• In that case, the method I used in my answer below should work fine. My example take 9 seconds for 1 million terms, so your slightly more complicated example should probably take 1 minute or less to compute. – DumpsterDoofus Jan 10 '15 at 1:05

1 Answer

This is the kind of problem which is most efficiently solved via memoization. As I mentioned in a comment, defining $R_N=\sum_{i=0}^{N-1}Q_i$, we have

\begin{align} Q_N&=f(R_N,k_n)\\ R_N&=R_{N-1}+Q_{N-1} \end{align}

where $f$ is some function which contains the specifics. In your case,

$$f(R_N,k_N)=\alpha\frac{\beta-\delta R_N}{k_N}\sqrt{\left(\frac{k_N}{\gamma}\frac{\beta+\delta R_N}{\beta-\delta R_N}\right)^{\phi}-\left(\frac{k_N}{\gamma}\frac{\beta+\delta R_N}{\beta-\delta R_N}\right)^{\theta}}$$

but I'll use the simpler example

\begin{align} f(R_N,k_N)& =\sin\left(R_N+k_N\right)\\ k_N& =N. \end{align} The following implements the memoization:

Clear[Q, R]
Q[0] = R[0] = 0;
Q[n_] := Q[n] = Sin[1.0 R[n] + n]
R[n_] := R[n] = R[n - 1] + Q[n - 1]


The following computes the first 100,000 values of $Q_N$ and stores their timings in a list, which is then plotted (after smoothing to remove jitter):

times = Table[First@AbsoluteTiming@Q[i], {i, 100000}];
ListLinePlot[GaussianFilter[times, 1000], PlotRange -> All]


As you can see, the performance is $O(1)$ and is thus asymptotically optimal. On my computer, computing 100,000 values takes 0.9 seconds, so you can compute the first 1,000,000,000 values of $Q_N$ in 2.5 hours.

It is straightforward to modify the code to work on your example (you just have to change the definition of Q[n]). If you need more than a billion values of $Q_N$, you may want to consider using a different programming language, or ask a separate question on how to performance-tune memoization.